bbbales2,
I see, that makes sense to me, thank you.
I coded a toy example to test the viability of the secant solver; the program did compile and run, but came back with several warnings-- Which doesn’t seem to me to be of too much concern, but as a new user to Stan (and MCMC in general), I am unsure if they indicate a more deeply rooted issue or just bad code on my part.
Some context: I intend to use the output B as the input to something else, i.e. I am not actually interested in how B is distributed, but I wanted to see if this small segment of root-finding code in Stan is well-behaved before nesting it in another function. However, I am unsure whether this toy example is sufficient to conclude that autodiff is working as intended.
[EDIT] I just read the documentation on the first warning regarding the divergent transition, I’ll try raising the adept_delta next time, but I am not sure if it’s something I’d need to worry about moving forward.
[EDIT 2] Thinking about the problem, it makes sense that many samples were divergent (if I understand the sampling process correctly), especially since the solutions to B is exact (numerical errors aside of course) because B isn’t sourced from data i.e. there isn’t any random noise to contend with. So in that sense, the values of lo_B and up_B will be well defined, and the sampler will “step off the edge” more often than not because there isn’t any scatter in the “data” for B, i.e. lo_B and up_B’s distributions will be a sharp peak centered at their respective values. Does that sound right?
Thank you for your time,
Chiot
Code in R:
library("rstan")
set.seed(1234)
N = 50
Ci <- runif(N, min = 12000, max = 12920)
C <- runif(N, min = 1.4, max = 1.6)
E <- runif(N, min = 0.7, max = 0.9)
init_val <- function () list(lo_B = 0.05, hi_B = 0.2)
# Example scenario, want to generate values of B
fit_B <- stan(
file = "Z01_algebra_solver_example.stan", # Stan program
data=c("N","Ci","C","E"), #algorithm="Fixed_param" when using monte carlo simulation,
chains = 4, # number of Markov chains
warmup = 1000, # number of warmup iterations per chain
iter = 5000, # total number of iterations per chain
cores = 2, # number of cores (could use one per chain)
refresh = 0, # no progress shown
init = init_val
)
options(max.print=100000)
print(fit_B, pars = c("B"), include = FALSE)
The code in Stan:
functions {
real B_root(real Ci, real B, real C, real E, real K, real f, real Fz){
real root_eqx;
root_eqx = Ci - (B*C*K*f*Fz)/(sin(C*atan(B*((1-E)*K + (E/B)*atan(B*K)))));
return root_eqx;
}
// Secant method root solver
real secant_root(real Ci, real C, real E, real K, real f, real Fz){
real x0 = 0.01;
real x1 = 0.5;
real x2;
real f0;
real f1;
for (i in 1:10){
f0 = B_root(Ci, x0, C, E, K, f, Fz);
f1 = B_root(Ci, x1, C, E, K, f, Fz);
x2 = x1 - ((x0 - x1)*f1)/(f0 - f1);
if ( fabs(x2 - x1) < 1E-10 ){
break;
}
x0 = x1;
x1 = x2;
}
return x2;
}
real shell_func(real Ci, real C, real E, real K, real f, real Fz){
real B_out = secant_root(Ci, C, E, K, f, Fz);
return B_out;
}
}
data {
int<lower = 0> N;
vector[N] Ci;
vector[N] C;
vector[N] E;
}
parameters {
real lo_B;
real up_B;
}
transformed parameters {
real K = 100;
real Fz = 677.9820;
real<lower=0> f = 0.9;
vector[N] B;
for (i in 1:N){
B[i] = shell_func(Ci[i], C[i], E[i], K, f, Fz);
}
}
model {
lo_B ~ gamma(0.001,0.001);
up_B ~ gamma(0.001,0.001);
B ~ uniform(lo_B, up_B);
}
generated quantities {
//real Ci = uniform_rng(12000,12920);
//real C = uniform_rng(1.4,1.6);
//real E = uniform_rng(0.7,0.9);
//real B = shell_func(Ci, C, E, K, f, Fz);
}
The output and warnings in R:
DIAGNOSTIC(S) FROM PARSER:
Info:
Left-hand side of sampling statement (~) may contain a non-linear transform of a parameter or local variable.
If it does, you need to include a target += statement with the log absolute determinant of the Jacobian of the transform.
Left-hand-side of sampling statement:
B ~ uniform(...)
Warning messages:
1: There were 13984 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
2: Examine the pairs() plot to diagnose sampling problems
>
> options(max.print=100000)
> print(fit_B, pars = c("B"), include = FALSE)
Inference for Stan model: Z01_algebra_solver_example.
4 chains, each with iter=5000; warmup=1000; thin=1;
post-warmup draws per chain=4000, total post-warmup draws=16000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
lo_B 0.11 0.00 0.00 0.11 0.11 0.11 0.11 0.11 648 1.00
up_B 0.14 0.00 0.00 0.14 0.14 0.14 0.14 0.14 699 1.01
K 100.00 NaN 0.00 100.00 100.00 100.00 100.00 100.00 NaN NaN
Fz 677.98 0.00 0.00 677.98 677.98 677.98 677.98 677.98 2 1.00
f 0.90 0.00 0.00 0.90 0.90 0.90 0.90 0.90 2 1.00
lp__ 177.26 0.06 1.48 173.60 176.53 177.59 178.35 179.14 577 1.01
Samples were drawn using NUTS(diag_e) at Sun Oct 27 16:31:41 2019.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).