Fit models for use in examples

Source: R/example.R

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 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.
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__     -65.9   -65.6   1.41   1.18   -68.8    -64.2     1.00    1006.
#>  2 alpha      0.385   0.382 0.220  0.217    0.0236   0.761   1.00    2337.
#>  3 beta[1]   -0.662  -0.654 0.242  0.241   -1.06    -0.276   1.00    1987.
#>  4 beta[2]   -0.279  -0.278 0.221  0.216   -0.648    0.0864  1.00    1949.
#>  5 beta[3]    0.680   0.673 0.266  0.271    0.270    1.12    1.00    1976.
#>  6 log_lik…  -0.514  -0.509 0.0984 0.0965  -0.687   -0.361   1.00    2280.
#>  7 log_lik…  -0.406  -0.389 0.146  0.137   -0.666   -0.204   1.00    2144.
#>  8 log_lik…  -0.499  -0.463 0.213  0.202   -0.903   -0.207   1.00    2098.
#>  9 log_lik…  -0.449  -0.438 0.148  0.146   -0.717   -0.233   1.00    1869.
#> 10 log_lik…  -1.18   -1.16  0.280  0.273   -1.68    -0.763   1.00    2465.
#> # … 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.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 
#> # … 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.2   -65.9   1.56   1.46   -69.2    -64.4  
#>  2 lp_approx__  -2.03   -1.72  1.44   1.27    -4.83    -0.311
#>  3 alpha         0.410   0.406 0.230  0.221    0.0365   0.806
#>  4 beta[1]      -0.790  -0.795 0.222  0.220   -1.16    -0.401
#>  5 beta[2]      -0.219  -0.228 0.248  0.262   -0.618    0.196
#>  6 beta[3]       0.761   0.753 0.252  0.246    0.333    1.16 
#>  7 log_lik[1]   -0.488  -0.482 0.0980 0.0954  -0.660   -0.338
#>  8 log_lik[2]   -0.350  -0.328 0.143  0.126   -0.615   -0.160
#>  9 log_lik[3]   -0.411  -0.373 0.204  0.192   -0.807   -0.151
#> 10 log_lik[4]   -0.450  -0.431 0.148  0.144   -0.725   -0.239
#> # … 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: 142 of 4000 (4.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.2  -58.6   5.26  5.39 -66.4    -48.8  1.03     141.     80.0
#>  2 mu         6.35   6.09  4.26  4.22  -0.616   13.3  1.01     761.   1722. 
#>  3 tau        5.38   4.60  3.56  3.34   1.18    12.3  1.04     129.     47.1
#>  4 theta[1]   9.13   8.26  7.12  6.34  -0.994   21.9  1.01     672.   1870. 
#>  5 theta[2]   6.64   6.42  5.72  5.18  -2.68    16.1  1.01    1598.   2092. 
#>  6 theta[3]   5.27   5.27  6.73  5.94  -6.14    16.0  1.01    1434.   2085. 
#>  7 theta[4]   6.43   6.19  5.94  5.14  -3.29    16.3  1.02    1352.   1663. 
#>  8 theta[5]   4.46   4.86  5.62  5.31  -5.27    13.1  1.01    1219.   2009. 
#>  9 theta[6]   5.18   5.40  5.97  5.37  -5.28    14.4  1.02    1243.   1882. 
#> 10 theta[7]   9.02   8.53  6.21  5.90  -0.0847  20.0  1.01     736.   1886. 
#> 11 theta[8]   7.03   6.49  6.92  5.67  -3.58    19.0  1.01    1352.   1806. 

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...
#> 
#> Warning: 3 of 4000 (0.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_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__     -46.9    -46.6    2.47  2.38  -51.3   -43.4   1.00    1468.    2252.
#>  2 mu         6.40     6.43   4.32  4.29   -0.516  13.4   1.00    2844.    2415.
#>  3 tau        4.80     4.05   3.62  3.42    0.431  11.9   1.00    1955.    1843.
#>  4 theta_r…   0.369    0.355  0.953 0.951  -1.23    1.95  1.00    3487.    2595.
#>  5 theta_r…   0.0449   0.0320 0.941 0.891  -1.46    1.60  1.00    3600.    2546.
#>  6 theta_r…  -0.123   -0.116  0.954 0.948  -1.68    1.45  1.00    4329.    2805.
#>  7 theta_r…   0.0187   0.0280 0.921 0.898  -1.52    1.50  1.00    3642.    2631.
#>  8 theta_r…  -0.285   -0.287  0.910 0.884  -1.79    1.22  1.00    3780.    2694.
#>  9 theta_r…  -0.158   -0.151  0.932 0.911  -1.69    1.37  1.00    3950.    2791.
#> 10 theta_r…   0.342    0.358  0.915 0.915  -1.16    1.79  1.00    3871.    3323.
#> 11 theta_r…   0.0502   0.0393 0.966 0.968  -1.51    1.63  1.00    4149.    2930.
#> 12 theta[1]   8.89     8.26   6.65  5.76   -0.834  20.7   1.00    3839.    3069.
#> 13 theta[2]   6.62     6.62   5.63  5.12   -2.71   15.7   1.00    3882.    2963.
#> 14 theta[3]   5.49     5.91   6.69  5.87   -6.12   15.7   1.00    3271.    2519.
#> 15 theta[4]   6.49     6.56   5.85  5.36   -2.93   16.0   1.00    3688.    3014.
#> 16 theta[5]   4.69     5.10   5.71  5.17   -5.21   13.4   1.00    3578.    2767.
#> 17 theta[6]   5.44     5.81   5.99  5.44   -4.86   14.7   1.00    3768.    2909.
#> 18 theta[7]   8.61     8.12   5.84  5.41   -0.231  18.8   1.00    3933.    3144.
#> 19 theta[8]   6.83     6.79   6.70  5.71   -3.72   17.8   1.00    3796.    2917.

# optimization fails for hierarchical model
cmdstanr_example("schools", "optimize", quiet = FALSE)
#> Model executable is up to date!
#> Initial log joint probability = -57.9499 
#>     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
#>       99       113.325      0.319253   2.69628e+09      0.5294      0.5294      169    
#>     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
#>      183       240.446    0.00518547   3.03565e+15       1e-12       0.001      411  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__       240.45
#>  mu           1.32
#>  tau          0.00
#>  theta[1]     1.32
#>  theta[2]     1.32
#>  theta[3]     1.32
#>  theta[4]     1.32
#>  theta[5]     1.32
#>  theta[6]     1.32
#>  theta[7]     1.32
#> 
#>  # showing 10 of 11 rows (change via 'max_rows' argument)
# }