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.45 1.20 -68.7 -64.3 1.00 934.
#> 2 alpha 0.384 0.386 0.218 0.218 0.0341 0.748 1.00 1959.
#> 3 beta[1] -0.667 -0.661 0.247 0.248 -1.08 -0.273 1.00 1840.
#> 4 beta[2] -0.275 -0.273 0.231 0.231 -0.662 0.114 1.00 1954.
#> 5 beta[3] 0.687 0.685 0.271 0.269 0.254 1.13 1.00 2212.
#> 6 log_lik[1] -0.515 -0.509 0.101 0.0967 -0.698 -0.360 1.00 2038.
#> 7 log_lik[2] -0.402 -0.380 0.150 0.140 -0.673 -0.195 1.00 2565.
#> 8 log_lik[3] -0.498 -0.467 0.223 0.212 -0.922 -0.201 1.00 2174.
#> 9 log_lik[4] -0.449 -0.430 0.154 0.148 -0.731 -0.232 0.999 1785.
#> 10 log_lik[5] -1.19 -1.17 0.283 0.281 -1.69 -0.772 1.00 2119.
#> # ℹ 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.5 -66.0 1.82 1.62 -69.9 -64.3
#> 2 lp_approx__ -1.96 -1.62 1.41 1.19 -4.62 -0.318
#> 3 alpha 0.503 0.492 0.219 0.222 0.144 0.876
#> 4 beta[1] -0.582 -0.568 0.301 0.288 -1.10 -0.0916
#> 5 beta[2] -0.278 -0.275 0.185 0.183 -0.592 0.0118
#> 6 beta[3] 0.688 0.696 0.301 0.307 0.211 1.18
#> 7 log_lik[1] -0.480 -0.475 0.0981 0.0994 -0.648 -0.332
#> 8 log_lik[2] -0.466 -0.443 0.190 0.185 -0.836 -0.197
#> 9 log_lik[3] -0.470 -0.443 0.195 0.188 -0.838 -0.209
#> 10 log_lik[4] -0.495 -0.480 0.152 0.148 -0.776 -0.277
#> # ℹ 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: 146 of 4000 (4.0%) transitions ended with a divergence.
#> 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.4 -58.5 5.06 5.42 -66.4 -50.2 1.02 139. 71.0
#> 2 mu 6.36 6.39 4.24 4.27 -0.684 13.0 1.01 487. 345.
#> 3 tau 5.42 4.61 3.53 3.33 1.25 12.3 1.02 136. 67.5
#> 4 theta[1] 9.11 8.41 7.20 6.29 -1.35 21.8 1.01 774. 601.
#> 5 theta[2] 6.65 6.69 5.88 5.59 -2.86 16.2 1.01 1011. 2093.
#> 6 theta[3] 5.38 5.61 6.74 6.18 -5.83 15.9 1.01 863. 1953.
#> 7 theta[4] 6.59 6.59 6.01 5.73 -2.94 16.6 1.01 946. 2109.
#> 8 theta[5] 4.37 4.72 5.66 5.51 -5.63 12.8 1.01 745. 2276.
#> 9 theta[6] 5.28 5.49 5.96 5.60 -4.99 14.6 1.01 879. 1751.
#> 10 theta[7] 9.06 8.65 6.17 5.72 -0.151 20.1 1.01 678. 1072.
#> 11 theta[8] 6.91 6.84 6.86 6.09 -4.10 17.8 1.01 939. 2286.
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.49 2.36 -51.6 -43.6 1.00 1229. 1844.
#> 2 mu 6.43 6.46e+0 4.08 3.93 -0.238 13.1 1.00 2985. 2197.
#> 3 tau 4.74 4.03e+0 3.70 3.41 0.334 11.9 1.00 1724. 1496.
#> 4 theta_r… 0.351 3.53e-1 0.979 0.955 -1.26 1.96 1.00 3404. 2385.
#> 5 theta_r… 0.0484 6.67e-2 0.924 0.913 -1.49 1.53 1.00 3680. 2435.
#> 6 theta_r… -0.131 -1.50e-1 0.962 0.938 -1.71 1.50 1.00 3067. 2557.
#> 7 theta_r… 0.0221 5.18e-3 0.919 0.905 -1.49 1.59 1.00 3748. 2583.
#> 8 theta_r… -0.273 -3.14e-1 0.931 0.892 -1.77 1.34 1.00 4139. 2891.
#> 9 theta_r… -0.157 -1.57e-1 0.930 0.929 -1.69 1.36 1.00 3560. 2797.
#> 10 theta_r… 0.363 4.01e-1 0.924 0.911 -1.15 1.87 1.00 3687. 2941.
#> 11 theta_r… 0.0558 7.25e-2 0.951 0.924 -1.55 1.65 1.00 4077. 2555.
#> 12 theta[1] 8.85 8.09e+0 6.85 5.77 -0.638 21.2 1.00 2994. 2531.
#> 13 theta[2] 6.79 6.75e+0 5.59 5.01 -2.46 15.9 1.00 4593. 3367.
#> 14 theta[3] 5.40 5.88e+0 6.64 5.63 -5.94 15.5 1.00 3448. 3181.
#> 15 theta[4] 6.65 6.62e+0 5.69 5.09 -2.25 16.0 1.00 4537. 3155.
#> 16 theta[5] 4.78 5.06e+0 5.54 5.26 -4.69 13.2 1.00 3996. 3191.
#> 17 theta[6] 5.44 5.63e+0 5.97 5.29 -4.66 14.8 1.00 3934. 2877.
#> 18 theta[7] 8.83 8.37e+0 5.98 5.56 0.141 19.4 1.00 4104. 3336.
#> 19 theta[8] 6.89 6.73e+0 6.34 5.33 -3.03 17.1 1.00 3988. 3025.
# optimization fails for hierarchical model
cmdstanr_example("schools", "optimize", quiet = FALSE)
#> Initial log joint probability = -56.2449
#> Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
#> 99 121.706 0.109668 2.48002e+09 1 1 174
#> Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
#> 187 244.785 0.253817 8.00296 1e-12 0.001 400 LS failed, Hessian reset
#> Line search failed to achieve a sufficient decrease, no more progress can be made
#> Chain 1 Optimization terminated with error:
#> Warning: Fitting finished unexpectedly! Use the $output() method for more information.
#> Error: Fitting failed. Unable to print.
# }