Basic Algebra Solver Example Code

Hi, I’ve been playing around with Stan in R with brms and eventually moved on to PyStan and then CmdStanPy for some simple linear regression models. I’ve gotten to the point that I want to incorporate some more advanced features of Stan, namely the Newton solver to find the root of an equation. I’ve read the Stan manual sections on the solver, but I have not been able to get the example discussed (12 Solving Algebraic Equations | Stan User’s Guide) to work.

I’m convinced I’m missing something really basic. For example, to get results in the algebra solver section, how should I specify the data, model, parameters block? Or do you even need to for the example given in the text. It’d help me a lot to understand if I could see a more detailed (full?) example of the Algebra Solver example in the manual. Thanks.

I had some time to experiment more, so I’ll add an update to this post. Looks like the answer to my question is that indeed the model and data blocks are necessary. I was able get the example from the Stan User Guide 12.3 to work, and I coded it up in a Jupyter notebook, see below and attached html.
solver_example.html (579.3 KB)

Solver Example with CmdStanPy
import os
import cmdstanpy
from cmdstanpy import cmdstan_path, CmdStanModel
import numpy as np

Specify the Stan model file

stan_file = './example_12_3_stan.stan'
with open(stan_file,'w') as f:
    f.write("""functions {

vector system(vector y,              // unknowns
              vector theta,          // parameters
              data array[] real x_r, // data (real)
              array[] int x_i) {     // data (integer)
  vector[2] z;
  z[1] = y[1] - theta[1];
  z[2] = y[1] * y[2] - theta[2];
  return z;
}

}

data {
   int N;
   vector[N] Y;
}

transformed data {
  vector[2] y_guess;
  y_guess[1] = 1;
  y_guess[2] = 1;
  array[1] real x_r;
  array[1] int x_i;
}

parameters {
  real alpha;
}

transformed parameters {
  vector[2] theta;
  theta[1] = 3;
  theta[2] = 6;
  vector[2] y;
  y = algebra_solver_newton(system, y_guess, theta, x_r, x_i);
}

model {

  Y ~ normal(alpha,0.1);
 

}

""")

Specify the input data dictionary

y = np.random.normal(loc=0,scale=0.1,size=10)
data_dict = {'N':10, 'Y':y}

Instantiate the model; compiles the Stan program as needed.

model = CmdStanModel(stan_file=stan_file)

Fit the model

fit = model.sample(data=data_dict)

Observe the summary

fit.summary()

Solver found solution to system of equation and values are stored in vector y. y[1] and y[2] are the solutions 3.0 and 2.0 as described in the example text.

If you’re just trying to run the algebra_solver, but not for evaluating the target density in Markov chain Monte Carlo, then you can use the fixed_param sampler which doesn’t require a parameters or model block. See for example the “*_tune.stan” programs in Some Ruminations on Containment Prior Modeling. The configuration of the fixed_param sampler will be similar although slightly different for CmdStanPy.