Solving one algebraic equation repeatedly with different data input

Modeling
1

I’m trying to use the algebraic solver in Stan to solve an implicit equation about y,

\frac{d_1}{\theta_3 \left( \frac{y}{|\theta_1 - y|} \right) ^{1/\theta_5}} + \frac{d_2}{\theta_4 \left( \frac{y}{|\theta_2 - y|} \right) ^{1/\theta_6}} - \theta_8 = 0,

where \theta's are the parameters, and d_1 and d_2 are data points.

The data d_1 and d_2 look like the following,
d1 d2
0.1 0.1
0.5 0.1
1.0 0.1
2.0 0.1
3.0 0.1
10.0 0.1
0.1 0.5
0.5 0.5
…
I want to solve for y for each row of this data with the same set of parameter values.

Can anybody give me some idea how to realize it in Stan?
Any help will be truly appreciated.

Angel

2

If you have N data points, just do N different solves. Check “20.3. Calling the Algebraic Solver” in the 2.17 manual.

If the code for one solve is:

data {
  real x_r[2];
}
...
transformed parameters {
  vector[1] y;
  y = algebra_solver(system, y_guess, theta, x_r, x_i);
}

You’ll just want to write something like (leaving out a lot of the model):

data {
  int N;
  real x_r[N, 2];
}
...
transformed parameters {
  vector[1] y[N];
  for(n in 1:N) {
    y[n] = algebra_solver(system, y_guess, theta, x_r[n], x_i);
  }
}
3

Thanks so much for your reply, Ben! Okay, I understand that part. I’ve been building code for one solve, and am having trouble getting it right… Here below is my stan code. I assume x_r is my d1 and d2 (i.e. data) that I pass to the system function. Would you mind taking a quick look at it? The error message is “PARSER EXPECTED: whitespace to end of file. FOUND AT line 2:”, but I suspect it is probably something else that causes the error. Thanks a lot!!

laModel <- "
vector system(vector y, vector theta, real[] x_r, int[] x_i){
  real z;
  z = x_r[1]/(theta[3]*(y/fabs(theta[1] - y))^(1/theta[5])) + x_r[2]/(theta[4]*(y/fabs(theta[2] - y))^(1/theta[6])) - theta[7];
  return z;
}

data{
  real x_r[2];
}

transformed data {
  vector y_guess = 0.6;
  //real x_r[0];
  int x_i[0];
}

transformed parameters {
  vector[7] theta = {1, 1, 7.3891, 3.4903, 1, 2, 0.4966};
  vector y;
  y = algebra_solver(system, y_guess, theta, x_r, x_i);
}

"

stan_laModel <- stan_model(model_code = laModel)
4

There were a few syntax errors. I changed it so it’d compile, but I’m not sure it’s correct from a modeling perspective.

functions { // You need to define functions inside a functions block
  vector system(vector y, vector theta, real[] x_r, int[] x_i){
    vector[1] z; // Vector sizes must be declared
    z[1] = x_r[1]/(theta[3]*(y[1]/fabs(theta[1] - y[1]))^(1/theta[5])) + x_r[2]/(theta[4]*(y[1]/fabs(theta[2] - y[1]))^(1/theta[6])) - theta[7];
    return z; // On the previous line, there is a vectorized version of ^, so I added [1] indices to the y variables. Not sure if that is right
  }
}

data {
  real x_r[2];
}

transformed data {
  vector[1] y_guess = [ 0.6 ]'; // Vector sizes must be declared
  //real x_r[0];
  int x_i[0];
}

transformed parameters {
  vector[7] theta = [ 1, 1, 7.3891, 3.4903, 1, 2, 0.4966 ]'; // I changed the initialization here
  vector[1] y; // Vector sizes must be declared
  y = algebra_solver(system, y_guess, theta, x_r, x_i);
}
5

Thanks a lot! It is such a great help. I updated the stan code and it compiled. Then I proceeded with the R code to invoke it, like following

stan_laModel <- stan_model(model_code = laModel)

stanFit <- function(model_obj){
  dataList <- list(x_r = c(0.1, 0.1))
  stanFit <- sampling(object = model_obj, 
                      data   = dataList,
                      #init = initf1,
                      algorithm="Fixed_param",
                      seed = 42,
                      chains = 1,
                      iter =1000) #include/pars opposite
  
  chains <- data.frame(rstan::extract(stanFit))
  return(list("stanFit" = stanFit, "df" = chains))
}

la_modelFit <- stanFit(stan_laModel)

The error message suggests to increase the relative tolerance and max_num_steps. Do you have any advice about the relative tolenrance and ma_num_steps? Or the actual problem is something else? Thanks again!

6

If the optimization is failing to converge, then the first thing to do is verify that there are no bugs in the function you’re optimizing and that an optimum of some sort actually exists.

Probly how I’d go about debugging this is by first using expose_stan_functions to expose the function you’re trying to optimize and play around with it in R (using the optim function, for instance) to see if you can reliably find optima.