Hi all stan users and developers,
Does someone know how to find roots for a equation of the form CDF_fun(x) - prob = 0 in Rstan?
Thank you
1 Like
Hi,
you are in luck - Stan has you covered via the algebraic function solver: 9.1 Algebraic equation solver | Stan Functions Reference
Feel free to ask for clarifications if you have trouble using it - it is one of the more tricky parts of Stan language.
Best of luck with your model!
1 Like
Hi @martinmodrak ,
Thank you for the response. I have used algebra_solver_newton function in my weird code. But, I am still struggling to make my model. I am trying to do Bayesian calibration for a Microsimulation model using stan (I could develop and run the code written in R using Metropolis Algorithm, but it’s a very slow converging process). As I have several _rng functions in my code I tried to define them within a user-defined _rng function.
Please see the attached .stan file and then kindly help me to develop it.
Thank you
I looked at the Stan code, but I admit I have no idea on what the model represents and what it should do. Could you be a bit more specific on where you’ve encoutered problems?
Hi @martinmodrak ,
Okay, sorry for my bad explanation. Simply, I am trying to use the simulation algorithm described in the below paper to do individual-level data simulations, and then it gives me a data table with information relevant to disease stages of lung cancer patients. I can easily do it in R.
ijm-00164.pdf (1.4 MB)
But, I need to use Bayesian calibration for parameter estimation of the model. For that, I am trying to run the simulation algorithm 10 times (say) and find the mean number of diagnosed individuals in each subgroup to feed the likelihood of my model (i.e. data_rng function in user-defined functions block and it returns lambda, that is mean number of diagnosed individuals in each subgroup).
My question is how to link this lambda value into the model block? Or do I have to run the simulation code in R separately to find lambda and then use it as data values for cmdstan() or rstan() funtions?
Thank you.
I admit I still don’t really understand what you are trying to do. What exactly do you mean by “Bayesian calibration” - the only sense in which I am familiar with this term is something like Simulation-based calibration, but that appears to be different to what you are trying to do.
Are you trying to estimate the parameters of the simulator? Or are you estimating something else?
Stan needs the likelihood to be completely deterministic. If you only need to generate lambda once and then use the same value in all likelihood evaluations, you can do that in transformed data block (that would be equivalent to running the code in R and passing it as data). However, within the model block, everything needs to be deterministic, so you cannot use a stochastic simulator there.
1 Like
Hi @martinmodrak ,
I am trying to figure out the CDF function of the following function using stan.
h(t) = \dfrac{\nu \mu X[exp(\gamma+2B).t - 1]}{\gamma+B[exp(\gamma+2B).t +1]}
Here, I defined CDF as
CDF = 1 - exp(-\int_0^{t}h(u)du
and then need to find t values when CDF has a uniform distribution.
Can I use algebraic_solver and integrate_1d within a same user defined function in order to pass into the model block in stan?
Can someone help me to write stan code here?
Thank you
Yes you can, but it will likely be very slow.
We did our best general explanations in the Stan documentation for those functions and the relevant chapters of Stan User’s guide. If you feel like you have a more specific issue, feel free to ask for additional advice.
Best of luck with your model!
Hi Martin,
Thank you for the response. I am very keen to see more examples and techniques for Rstan as I am planning to use Rstan in my whole project. How can I collect specific Rstan issues?
I have a question, is it possible to call integrator within another user-defined function instead of in the model block?
Something like that…
functions{
real ht_func(real x, real xc, real[] theta, real[] x_r, int[] x_i){
real NC = theta[1];
real mu = theta[2];
real nu = theta[2];
real alpha1 = theta[3];
real alpha2 = theta[4];
real gamma_ns = theta[5];
real alpha_ns = theta[6];
real q = theta[7];
real ht;
real alpha = alpha_ns * (1 + alpha1* pow(q,alpha2));
real gamma = gamma_ns * (1 + alpha1* pow(q,alpha2));
real B = (-gamma + sqrt(gamma^2 + 4 * alpha * mu))/2;
ht = (nu * mu * NC * (exp(gamma + 2* B) * x -1))/(gamma + (B *(exp(gamma + 2* B) * x + 1)));
return ht;
}
vector HT_func(vector x, vector theta, real[] x1_r, int[] x1_i){
real NC = theta[1];
real mu = theta[2];
real nu = theta[2];
real alpha1 = theta[3];
real alpha2 = theta[4];
real gamma_ns = theta[5];
real alpha_ns = theta[6];
real q = theta[7];
real t;
vector[1] Ht;
Ht = intrgrate_1d(ht_func, 0, t, {NC, mu, nu, alpha1, alpha2, gamma_ns, alpha_ns, q}, {0.0}, {0}, 1e-8);
return Ht;
}
}
Thank you
Yes, that should definitely be possible.
1 Like
Hi @martinmodrak ,
There is another weird thing for me to understand when applying Algebraic equation solver . In my case, algebra_solver_newton doesn’t recognize the function (i.e. function does not exit) but algebra_solver does for the same argumunets.
Do you know why does it happen and any solutions ?
Also, are there any extra references with examples related to algebraic equation solver ?
Thank you.
Most likely, this is because you are using an older version of Stan, algebra_solver_newton was AFAIK introduced around Stan 2.26 (not sure exactly). If you are using rstan, you may need to use Repository for distributing (some) stan-dev R packages | r-packages as the version on CRAN is not updated currently (for stupid reasons, we are trying to resolve).
1 Like
The latest rstan version is 2.21.2, isn’t it? (which is available in CRAN and mine as well). Are there other new versions which compatible with algebra_solver_newton because version 2.21.2 didn’t help me to find roots.
Thank you
Sorry, forgot about this thread.
The latest r-stan version is 2.27 and it is available at Repository for distributing (some) stan-dev R packages | r-packages Unfortunately, the latest CRAN version is still 2.21.2 due to somewhat stupid reasons).
1 Like