Fit models for use in examples

cmdstanr_example(
  example = c("logistic", "schools", "schools_ncp"),
  method = c("sample", "optimize", "variational"),
  ...,
  quiet = TRUE
)

print_example_program(example = c("logistic", "schools", "schools_ncp"))

Arguments

example	Name of the example. The currently available examples are: `"logistic"`: logistic regression with intercept and 3 predictors. `"schools"`: the so-called "eight schools" model, a hierarchical meta-analysis. Fitting this model will result in warnings about divergences. `"schools_ncp"`: non-centered parameterization of the "eight schools" model that fixes the problem with divergences. To print the Stan code for a given example use `print_example_program()`.
method	Which fitting method should be used? The default is the `"sample"` method (MCMC).
...	Arguments passed to the chosen `method`. See the help pages for the individual methods for details.
quiet	If `TRUE` (the default) then fitting the model is wrapped in `utils::capture.output()`.

Value

The fitted model object returned by the selected method.

Examples

# \dontrun{
print_example_program("logistic")
#> data {
#>   int<lower=0> N;
#>   int<lower=0> K;
#>   int<lower=0,upper=1> y[N];
#>   matrix[N, K] X;
#> }
#> parameters {
#>   real alpha;
#>   vector[K] beta;
#> }
#> model {
#>   target += normal_lpdf(alpha | 0, 1);
#>   target += normal_lpdf(beta | 0, 1);
#>   target += bernoulli_logit_glm_lpmf(y | X, alpha, beta);
#> }
#> generated quantities {
#>   vector[N] log_lik;
#>   for (n in 1:N) log_lik[n] = bernoulli_logit_lpmf(y[n] | alpha + X[n] * beta);
#> }
fit_logistic_mcmc <- cmdstanr_example("logistic", chains = 2)
#> Compiling Stan program...
#> 
-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

 
fit_logistic_mcmc$summary()
#> # A tibble: 105 x 10
#>    variable    mean  median     sd    mad       q5      q95  rhat ess_bulk
#>    <chr>      <dbl>   <dbl>  <dbl>  <dbl>    <dbl>    <dbl> <dbl>    <dbl>
#>  1 lp__     -66.1   -65.7   1.50   1.30   -68.9    -64.3    1.00     1034.
#>  2 alpha      0.375   0.376 0.219  0.215    0.0119   0.727  1.00     1869.
#>  3 beta[1]   -0.664  -0.664 0.251  0.244   -1.09    -0.247  0.999    1653.
#>  4 beta[2]   -0.278  -0.271 0.235  0.247   -0.672    0.0868 1.00     1539.
#>  5 beta[3]    0.672   0.669 0.279  0.275    0.221    1.16   1.00     1884.
#>  6 log_lik…  -0.518  -0.510 0.0992 0.0997  -0.691   -0.365  1.00     1844.
#>  7 log_lik…  -0.407  -0.387 0.154  0.144   -0.694   -0.194  1.00     2025.
#>  8 log_lik…  -0.506  -0.471 0.224  0.207   -0.912   -0.209  1.00     1754.
#>  9 log_lik…  -0.449  -0.430 0.152  0.149   -0.722   -0.236  1.00     1755.
#> 10 log_lik…  -1.17   -1.14  0.286  0.281   -1.69    -0.744  1.00     2105.
#> # … with 95 more rows, and 1 more variable: ess_tail <dbl>

fit_logistic_optim <- cmdstanr_example("logistic", method = "optimize")
#> Model executable is up to date!
fit_logistic_optim$summary()
#> # A tibble: 105 x 2
#>    variable   estimate
#>    <chr>         <dbl>
#>  1 lp__        -63.9  
#>  2 alpha         0.364
#>  3 beta[1]      -0.632
#>  4 beta[2]      -0.259
#>  5 beta[3]       0.649
#>  6 log_lik[1]   -0.515
#>  7 log_lik[2]   -0.394
#>  8 log_lik[3]   -0.469
#>  9 log_lik[4]   -0.442
#> 10 log_lik[5]   -1.14 
#> # … with 95 more rows

fit_logistic_vb <- cmdstanr_example("logistic", method = "variational")
#> Model executable is up to date!
fit_logistic_vb$summary()
#> # A tibble: 106 x 7
#>    variable       mean  median    sd   mad        q5     q95
#>    <chr>         <dbl>   <dbl> <dbl> <dbl>     <dbl>   <dbl>
#>  1 lp__        -66.3   -65.9   1.85  1.44  -70.0     -64.3  
#>  2 lp_approx__  -2.00   -1.58  1.53  1.14   -4.95     -0.334
#>  3 alpha         0.357   0.355 0.225 0.227  -0.00322   0.727
#>  4 beta[1]      -0.560  -0.566 0.206 0.203  -0.900    -0.199
#>  5 beta[2]      -0.283  -0.291 0.247 0.246  -0.671     0.118
#>  6 beta[3]       0.682   0.684 0.314 0.300   0.157     1.19 
#>  7 log_lik[1]   -0.540  -0.529 0.102 0.103  -0.720    -0.392
#>  8 log_lik[2]   -0.424  -0.395 0.186 0.164  -0.773    -0.187
#>  9 log_lik[3]   -0.548  -0.501 0.257 0.236  -1.03     -0.215
#> 10 log_lik[4]   -0.447  -0.426 0.156 0.151  -0.728    -0.229
#> # … with 96 more rows

print_example_program("schools")
#> data {
#>   int<lower=1> J;
#>   vector<lower=0>[J] sigma;
#>   vector[J] y;
#> }
#> parameters {
#>   real mu;
#>   real<lower=0> tau;
#>   vector[J] theta;
#> }
#> model {
#>   target += normal_lpdf(tau | 0, 10);
#>   target += normal_lpdf(mu | 0, 10);
#>   target += normal_lpdf(theta | mu, tau);
#>   target += normal_lpdf(y | theta, sigma);
#> }
fit_schools_mcmc <- cmdstanr_example("schools")
#> Compiling Stan program...
#> 
#> Warning: 116 of 4000 (3.0%) transitions ended with a divergence.
#> This may indicate insufficient exploration of the posterior distribution.
#> Possible remedies include: 
#>   * Increasing adapt_delta closer to 1 (default is 0.8) 
#>   * Reparameterizing the model (e.g. using a non-centered parameterization)
#>   * Using informative or weakly informative prior distributions 
#> 
-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

 
fit_schools_mcmc$summary()
#> # A tibble: 11 x 10
#>    variable   mean median    sd   mad      q5   q95  rhat ess_bulk ess_tail
#>    <chr>     <dbl>  <dbl> <dbl> <dbl>   <dbl> <dbl> <dbl>    <dbl>    <dbl>
#>  1 lp__     -58.1  -58.4   5.01  5.00 -66.0   -49.4  1.02     151.     146.
#>  2 mu         6.46   6.48  4.17  4.20  -0.430  13.3  1.01     583.     892.
#>  3 tau        5.31   4.57  3.46  3.16   1.29   12.0  1.02     149.     117.
#>  4 theta[1]   9.03   8.31  6.55  5.78  -0.315  20.9  1.00     896.    1634.
#>  5 theta[2]   6.85   6.77  5.78  5.41  -2.59   16.5  1.00    1153.    2092.
#>  6 theta[3]   5.33   5.72  6.49  5.78  -5.50   15.2  1.01     909.    1843.
#>  7 theta[4]   6.68   6.69  5.93  5.62  -2.88   16.2  1.01    1018.    1835.
#>  8 theta[5]   4.58   4.91  5.69  5.30  -5.62   13.5  1.00     821.    1447.
#>  9 theta[6]   5.47   5.61  5.97  5.43  -4.73   14.9  1.00    1010.    2114.
#> 10 theta[7]   9.11   8.84  5.92  5.64   0.279  19.4  1.01     647.    1611.
#> 11 theta[8]   7.11   6.93  6.50  5.67  -3.15   17.9  1.00    1103.    1891.

print_example_program("schools_ncp")
#> data {
#>   int<lower=1> J;
#>   vector<lower=0>[J] sigma;
#>   vector[J] y;
#> }
#> parameters {
#>   real mu;
#>   real<lower=0> tau;
#>   vector[J] theta_raw;
#> }
#> transformed parameters {
#>   vector[J] theta = mu + tau * theta_raw;
#> }
#> model {
#>   target += normal_lpdf(tau | 0, 10);
#>   target += normal_lpdf(mu | 0, 10);
#>   target += normal_lpdf(theta_raw | 0, 1);
#>   target += normal_lpdf(y | theta, sigma);
#> }
fit_schools_ncp_mcmc <- cmdstanr_example("schools_ncp")
#> Compiling Stan program...
#> 
-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

 
fit_schools_ncp_mcmc$summary()
#> # A tibble: 19 x 10
#>    variable     mean   median    sd   mad      q5    q95  rhat ess_bulk ess_tail
#>    <chr>       <dbl>    <dbl> <dbl> <dbl>   <dbl>  <dbl> <dbl>    <dbl>    <dbl>
#>  1 lp__     -4.69e+1 -4.65e+1 2.51  2.41  -51.5   -43.5   1.00    1398.    1954.
#>  2 mu        6.59e+0  6.59e+0 4.17  4.18   -0.509  13.5   1.00    3431.    2695.
#>  3 tau       4.84e+0  4.02e+0 3.75  3.61    0.369  12.1   1.00    2089.    1665.
#>  4 theta_r…  3.65e-1  3.82e-1 0.947 0.936  -1.21    1.90  1.00    4626.    2859.
#>  5 theta_r…  3.27e-2  3.44e-2 0.902 0.889  -1.44    1.54  1.00    4825.    3037.
#>  6 theta_r… -1.53e-1 -1.37e-1 0.967 0.977  -1.75    1.44  1.00    4805.    2670.
#>  7 theta_r… -3.93e-3 -1.73e-3 0.920 0.957  -1.51    1.50  1.00    4508.    2788.
#>  8 theta_r… -2.67e-1 -2.85e-1 0.909 0.906  -1.72    1.27  1.00    4716.    3210.
#>  9 theta_r… -1.58e-1 -1.74e-1 0.921 0.923  -1.61    1.39  1.00    4615.    2763.
#> 10 theta_r…  3.74e-1  3.90e-1 0.939 0.932  -1.18    1.90  1.00    4511.    2778.
#> 11 theta_r…  7.66e-2  8.40e-2 0.940 0.895  -1.51    1.65  1.00    4680.    2747.
#> 12 theta[1]  9.11e+0  8.56e+0 6.55  5.72   -0.396  20.4   1.00    4355.    3407.
#> 13 theta[2]  6.78e+0  6.76e+0 5.48  4.97   -2.36   15.8   1.00    4868.    3188.
#> 14 theta[3]  5.67e+0  5.87e+0 6.55  5.90   -5.20   15.7   1.00    4349.    3277.
#> 15 theta[4]  6.65e+0  6.67e+0 5.90  5.30   -3.04   16.3   1.00    4622.    3311.
#> 16 theta[5]  4.96e+0  5.27e+0 5.60  5.22   -4.77   13.6   1.00    4428.    3107.
#> 17 theta[6]  5.65e+0  5.88e+0 5.75  5.32   -4.10   14.5   1.00    4341.    3406.
#> 18 theta[7]  8.99e+0  8.47e+0 6.00  5.58    0.206  19.4   1.00    4601.    3579.
#> 19 theta[8]  7.18e+0  7.04e+0 6.43  5.56   -2.92   18.0   1.00    4291.    3486.

# optimization fails for hierarchical model
cmdstanr_example("schools", "optimize", quiet = FALSE)
#> Model executable is up to date!
#> Initial log joint probability = -55.0632 
#>     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
#>       99       108.449     0.0452483   8.52143e+08     0.09533     0.09533      150    
#>     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
#>      189       242.335      0.011797   5.10505e+15       1e-12       0.001      378  LS failed, Hessian reset  
#> Optimization terminated with error: 
#>   Line search failed to achieve a sufficient decrease, no more progress can be made
#> Finished in  0.1 seconds.
#>  variable estimate
#>  lp__       242.34
#>  mu           1.24
#>  tau          0.00
#>  theta[1]     1.24
#>  theta[2]     1.24
#>  theta[3]     1.24
#>  theta[4]     1.24
#>  theta[5]     1.24
#>  theta[6]     1.24
#>  theta[7]     1.24
#> 
#>  # showing 10 of 11 rows (change via 'max_rows' argument or 'cmdstanr_max_rows' option)
# }