Fit a Stan model with MCMC

Fit the eight schools hierarchical model with four MCMC chains, then print a posterior summary, the sampler’s warnings, and draws of tau.

Compiling a model needs CmdStan and a C++ toolchain, as described in Installation. The data are the estimated effects y of a coaching program on SAT-V scores in eight high schools and their standard errors sigma (Rubin 1981), read from schools.csv in the PHPstats datasets repository.

use function CmdStan\cmdstan_model;use function CmdStan\write_stan_file;$code = <<<'STAN'data {  int<lower=1> J;  vector[J] y;  vector<lower=0>[J] sigma;}parameters {  real mu;  real<lower=0> tau;  vector[J] theta_raw;}transformed parameters {  vector[J] theta = mu + tau * theta_raw;}model {  mu ~ normal(0, 10);  tau ~ normal(0, 10);  theta_raw ~ std_normal();  y ~ normal(theta, sigma);}STAN;$file = write_stan_file(    $code,    dir: __DIR__,    basename: 'schools',);$model = cmdstan_model($file);$url = 'https://raw.githubusercontent.com/'    . 'phpstats/datasets/main/csv/schools.csv';$lines = file($url, FILE_IGNORE_NEW_LINES);$rows = array_map(    fn($line) => str_getcsv($line, escape: ''),    array_slice($lines, 1),);$data = [    'J' => count($rows),    'y' => array_map(fn($r) => (float)$r[1], $rows),    'sigma' => array_map(fn($r) => (float)$r[2], $rows),];$fit = $model->sample(    data: $data,    seed: 123,    chains: 4,    parallel_chains: 4,    adapt_delta: 0.95,    show_messages: false,);$rows = $fit->summary(['mu', 'tau', 'theta']);echo "variable   mean    sd     q5    q95  rhat ess_bulk\n";foreach ($rows as $r) {    printf(        "%-8s %6.2f %5.2f %6.2f %6.2f %5.3f %8.0f\n",        $r['variable'],        $r['mean'],        $r['sd'],        $r['q5'],        $r['q95'],        $r['rhat'],        $r['ess_bulk'],    );}$warnings = $fit->diagnosticWarnings();echo $warnings ? implode("\n", $warnings) : 'No warnings';echo "\n";$draws = $fit->draws('tau', 'draws_matrix');$tau = array_column($draws, 'tau');echo count($tau), " draws of tau, the first five:\n";$first = array_slice($tau, 0, 5);printf("%.2f, %.2f, %.2f, %.2f, %.2f\n", ...$first);

Output:

variable   mean    sd     q5    q95  rhat ess_bulk
mu         6.40  4.06  -0.31  13.12 1.001     3582
tau        4.78  3.75   0.37  12.00 1.000     2176
theta[1]   8.87  6.90  -0.84  21.11 1.000     4184
theta[2]   6.72  5.47  -2.11  15.66 1.001     5706
theta[3]   5.37  6.35  -5.21  14.78 1.001     4311
theta[4]   6.64  5.63  -2.76  15.59 1.001     4613
theta[5]   4.83  5.50  -4.85  13.22 1.001     4452
theta[6]   5.46  5.79  -4.55  14.29 1.000     4957
theta[7]   8.75  5.91   0.12  19.22 1.002     4966
theta[8]   6.90  6.39  -3.13  17.47 1.001     4643
No warnings
4000 draws of tau, the first five:
1.41, 1.36, 1.63, 0.63, 0.06

• write_stan_file() rewrites schools.stan only when the program changes, so the model compiles on the first run and later runs reuse the executable.
• Four chains of 1,000 sampling iterations, after 1,000 warmup iterations, give 4,000 draws per variable.
• theta[1] to theta[8] are the effects in schools A to H, and mu and tau are the mean and standard deviation of the effects across schools, all in SAT-V points.
• Partial pooling pulls each school’s theta toward mu. The pull is stronger when tau is small or a school’s sigma is large.
• The non-centered form, theta = mu + tau * theta_raw, usually produces far fewer divergent transitions than sampling theta directly. It does not rule them out.
• Before using the estimates, check that rhat is below 1.01, ess_bulk is at least 100 per chain, and there are no warnings. Passing these checks does not prove convergence.
• A fixed seed reproduces the draws only with the same CmdStan version, compiler, and machine.

Compiling models · MCMC sampling · Draws, summaries, and saving fits · Functions and methods · Interactive version