Reproducing the parameter transformation example from the statmodeling blog

Modeling
1

I am trying to run the stan example from the statmodeling blog post on Transforming parameters in a time series model

parameters {
  real xi;
}
transformed parameters {
  real rho;
  rho = sqrt(1 - 1/xi);
}
model {
  target += - log(rho) + 2*log(1 - square(rho));
}

I am running this stan model through the following R code.

library(rstan)
print(stan("transformation_test.stan"))

But I get the following warnings and no matter how much I change the control parameters the model does not converge. At the bottom you can see that the mean of rho is not close to 0.5 and the number of samples is low at 48. The following result is on Mac.

Warning message:
“There were 2018 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup”Warning message:
“There were 278 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded”Warning message:
“There were 3 chains where the estimated Bayesian Fraction of Missing Information was low. See
http://mc-stan.org/misc/warnings.html#bfmi-low”Warning message:
“Examine the pairs() plot to diagnose sampling problems
”
Inference for Stan model: transformation_test.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

      mean se_mean   sd  2.5%   25%   50%  75% 97.5% n_eff Rhat
xi    4.59    0.79 8.61  1.01  1.21  1.75 3.73 26.87   118 1.04
rho   0.62    0.04 0.26  0.11  0.42  0.65 0.86  0.98    48 1.11
lp__ -1.20    0.31 2.29 -6.56 -2.48 -0.69 0.49  2.15    54 1.10

Samples were drawn using NUTS(diag_e) at Wed Jan 30 18:14:19 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).

I also tried to run the above on Rstudio cloud based on the instructions on the stan web page but the code failed to compile. The public R workspace is here

https://rstudio.cloud/project/193726

2

Okay figured out the problem by reading the comments, it seems that some constraints on the parameters are missing in the blog currently. Adding a lower value constraint fixes the convergence.

parameters {
  real<lower = 1> xi;
}
transformed parameters {
  real rho;
  rho = sqrt(1 - 1/xi);
}
model {
  target += - log(rho) + 2*log(1 - square(rho));
}
2 Likes