Fit models for use in examples
cmdstanr_example(
example = c("logistic", "schools", "schools_ncp"),
method = c("sample", "optimize", "laplace", "variational", "pathfinder", "diagnose"),
...,
quiet = TRUE,
force_recompile = getOption("cmdstanr_force_recompile", default = FALSE)
)
print_example_program(example = c("logistic", "schools", "schools_ncp"))
(string) The 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(example)
.
(string) 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.
(logical) If TRUE
(the default) then fitting the model is
wrapped in utils::capture.output()
.
Passed to the $compile() method.
The fitted model object returned by the selected method
.
# \dontrun{
print_example_program("logistic")
#> data {
#> int<lower=0> N;
#> int<lower=0> K;
#> array[N] int<lower=0, upper=1> y;
#> 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)
fit_logistic_mcmc$summary()
#> # A tibble: 105 × 10
#> variable mean median sd mad q5 q95 rhat ess_bulk
#> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 lp__ -66.0 -65.6 1.43 1.22 -68.8 -64.3 1.00 963.
#> 2 alpha 0.376 0.373 0.220 0.221 0.0226 0.745 1.00 2237.
#> 3 beta[1] -0.671 -0.663 0.250 0.252 -1.10 -0.258 1.00 1930.
#> 4 beta[2] -0.262 -0.261 0.222 0.229 -0.629 0.0934 1.00 1801.
#> 5 beta[3] 0.674 0.681 0.264 0.263 0.242 1.10 1.00 1979.
#> 6 log_lik[1] -0.515 -0.509 0.0991 0.0964 -0.681 -0.363 1.00 2100.
#> 7 log_lik[2] -0.401 -0.381 0.147 0.138 -0.671 -0.197 1.00 2010.
#> 8 log_lik[3] -0.489 -0.456 0.215 0.202 -0.891 -0.207 1.00 1914.
#> 9 log_lik[4] -0.455 -0.435 0.154 0.152 -0.726 -0.242 1.00 1918.
#> 10 log_lik[5] -1.18 -1.16 0.283 0.282 -1.69 -0.749 1.00 2470.
#> # ℹ 95 more rows
#> # ℹ 1 more variable: ess_tail <dbl>
fit_logistic_optim <- cmdstanr_example("logistic", method = "optimize")
fit_logistic_optim$summary()
#> # A tibble: 105 × 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.648
#> 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
#> # ℹ 95 more rows
fit_logistic_vb <- cmdstanr_example("logistic", method = "variational")
fit_logistic_vb$summary()
#> # A tibble: 106 × 7
#> variable mean median sd mad q5 q95
#> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 lp__ -66.4 -65.9 1.84 1.52 -69.9 -64.3
#> 2 lp_approx__ -1.98 -1.65 1.42 1.25 -4.64 -0.331
#> 3 alpha 0.377 0.370 0.296 0.307 -0.116 0.869
#> 4 beta[1] -0.646 -0.648 0.241 0.234 -1.04 -0.241
#> 5 beta[2] -0.252 -0.257 0.201 0.191 -0.579 0.0845
#> 6 beta[3] 0.702 0.695 0.280 0.269 0.236 1.16
#> 7 log_lik[1] -0.523 -0.518 0.128 0.130 -0.747 -0.331
#> 8 log_lik[2] -0.398 -0.369 0.168 0.155 -0.716 -0.174
#> 9 log_lik[3] -0.480 -0.450 0.204 0.193 -0.859 -0.211
#> 10 log_lik[4] -0.455 -0.431 0.159 0.160 -0.739 -0.235
#> # ℹ 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")
#> Warning: 260 of 4000 (6.0%) transitions ended with a divergence.
#> See https://mc-stan.org/misc/warnings for details.
#> Warning: 1 of 4 chains had an E-BFMI less than 0.3.
#> See https://mc-stan.org/misc/warnings for details.
fit_schools_mcmc$summary()
#> # A tibble: 11 × 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.0 -58.4 5.27 5.30 -66.3 -48.5 1.07 43.0 27.8
#> 2 mu 6.38 6.13 4.15 3.80 -0.237 13.5 1.01 627. 964.
#> 3 tau 5.31 4.49 3.57 3.28 1.10 12.2 1.07 37.6 21.2
#> 4 theta[1] 9.14 8.16 7.00 6.00 -0.685 21.7 1.01 966. 1427.
#> 5 theta[2] 6.73 6.51 5.64 4.88 -2.30 16.4 1.02 1127. 1732.
#> 6 theta[3] 5.08 5.36 6.58 5.64 -6.24 15.5 1.02 792. 1549.
#> 7 theta[4] 6.57 6.30 5.92 5.19 -2.79 16.1 1.01 1328. 1869.
#> 8 theta[5] 4.56 4.87 5.47 5.08 -5.00 13.1 1.02 707. 1473.
#> 9 theta[6] 5.31 5.53 6.01 5.11 -5.18 14.7 1.01 834. 1685.
#> 10 theta[7] 9.04 8.41 5.90 5.30 0.408 19.3 1.01 701. 1390.
#> 11 theta[8] 6.88 6.54 6.69 5.37 -4.06 18.2 1.02 1238. 2106.
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")
fit_schools_ncp_mcmc$summary()
#> # A tibble: 19 × 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__ -46.9 -4.66e+1 2.40 2.30 -51.2 -43.5 1.00 1569. 2358.
#> 2 mu 6.51 6.55e+0 4.22 4.14 -0.355 13.3 1.00 3009. 2465.
#> 3 tau 4.86 4.06e+0 3.67 3.43 0.444 12.0 1.00 1792. 1620.
#> 4 theta_r… 0.348 3.50e-1 0.956 0.960 -1.26 1.93 1.00 3527. 2677.
#> 5 theta_r… 0.0410 5.49e-2 0.889 0.880 -1.47 1.46 1.00 4110. 2991.
#> 6 theta_r… -0.142 -1.47e-1 0.953 0.945 -1.69 1.43 1.00 4686. 2441.
#> 7 theta_r… -0.0107 -9.40e-3 0.937 0.935 -1.53 1.53 1.00 4262. 2854.
#> 8 theta_r… -0.296 -2.93e-1 0.931 0.915 -1.83 1.24 1.00 3286. 2516.
#> 9 theta_r… -0.172 -1.84e-1 0.911 0.881 -1.69 1.36 1.00 3801. 2532.
#> 10 theta_r… 0.360 3.84e-1 0.933 0.904 -1.21 1.83 1.00 3900. 2825.
#> 11 theta_r… 0.0773 7.71e-2 0.982 1.00 -1.54 1.67 1.00 4239. 2766.
#> 12 theta[1] 9.01 8.06e+0 6.94 5.65 -0.546 22.1 1.00 3731. 3102.
#> 13 theta[2] 6.85 6.80e+0 5.61 5.13 -2.12 16.3 1.00 4538. 3046.
#> 14 theta[3] 5.47 5.82e+0 6.48 5.46 -6.06 15.3 1.00 4254. 3124.
#> 15 theta[4] 6.52 6.47e+0 5.74 5.18 -2.72 15.9 1.00 4625. 3465.
#> 16 theta[5] 4.79 5.10e+0 5.61 5.15 -4.87 13.3 1.00 4852. 3150.
#> 17 theta[6] 5.54 5.73e+0 5.60 5.11 -3.94 14.1 1.00 4021. 3062.
#> 18 theta[7] 8.95 8.27e+0 6.02 5.41 0.317 19.7 1.00 4087. 3495.
#> 19 theta[8] 7.00 7.00e+0 6.58 5.65 -3.39 17.8 1.00 3992. 3246.
# optimization fails for hierarchical model
cmdstanr_example("schools", "optimize", quiet = FALSE)
#> Initial log joint probability = -57.1999
#> Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
#> 99 137.364 0.389882 2.12196e+10 0.1758 0.3216 199
#> Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
#> 175 252.319 0.0285374 7.72538e+16 1e-12 0.001 386 LS failed, Hessian reset
#> Chain 1 Optimization terminated with error:
#> Chain 1 Line search failed to achieve a sufficient decrease, no more progress can be made
#> Warning: Fitting finished unexpectedly! Use the $output() method for more information.
#> Finished in 0.2 seconds.
#> Error: Fitting failed. Unable to print.
# }