# 4 MCMC Sampling

## 4.1 Running the Sampler

To generate a sample from the posterior distribution of
the model conditioned on the data,
we run the executable program with the argument `sample`

or `method=sample`

together with the input data.
The executable can be run from any directory.
Here, we run it in the directory which contains the Stan program and input data,
`<cmdstan-home>/examples/bernoulli`

:

`> cd examples/bernoulli`

To execute sampling of the model under Linux or Mac, use:

`> ./bernoulli sample data file=bernoulli.data.json`

In Windows, the `./`

prefix is not needed:

`> bernoulli.exe sample data file=bernoulli.data.json`

The output is the same across all supported platforms. First, the configuration of the program is echoed to the standard output:

```
method = sample (Default)
sample
num_samples = 1000 (Default)
num_warmup = 1000 (Default)
save_warmup = 0 (Default)
thin = 1 (Default)
adapt
engaged = 1 (Default)
gamma = 0.050000000000000003 (Default)
delta = 0.80000000000000004 (Default)
kappa = 0.75 (Default)
t0 = 10 (Default)
init_buffer = 75 (Default)
term_buffer = 50 (Default)
window = 25 (Default)
algorithm = hmc (Default)
hmc
engine = nuts (Default)
nuts
max_depth = 10 (Default)
metric = diag_e (Default)
metric_file = (Default)
stepsize = 1 (Default)
stepsize_jitter = 0 (Default)
id = 0 (Default)
data
file = bernoulli.data.json
init = 2 (Default)
random
seed = 3252652196 (Default)
output
file = output.csv (Default)
diagnostic_file = (Default)
refresh = 100 (Default)
```

After the configuration has been displayed, a short timing message is given.

```
Gradient evaluation took 1.2e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.12 seconds.
Adjust your expectations accordingly!
```

Next, the sampler reports the iteration number, reporting the percentage complete.

```
Iteration: 1 / 2000 [ 0%] (Warmup)
....
Iteration: 2000 / 2000 [100%] (Sampling)
```

Finally, the sampler reports timing information:

```
Elapsed Time: 0.007 seconds (Warm-up)
0.017 seconds (Sampling)
0.024 seconds (Total)
```

## 4.2 Running Multiple Chains

A Markov chain generates samples from the target distribution only after it has converged to equilibrium. In theory, convergence is only guaranteed asymptotically as the number of draws grows without bound. In practice, diagnostics must be applied to monitor convergence for the finite number of draws actually available. One way to monitor whether a chain has converged to the equilibrium distribution is to compare its behavior to other randomly initialized chains. For robust diagnostics, we recommend running 4 chains.

To run multiple chains given a model and data, either sequentially or in parallel,
we use the Unix or DOS shell `for`

loop to set up index variables needed to identify
each chain and its outputs.

On MacOS or Linux, the for-loop syntax for both the bash and zsh interpreters is:

`for NAME [in LIST]; do COMMANDS; done`

The list can be a simple sequence of numbers, or you can use the shell expansion syntax `{1..N}`

which expands to the sequence from \(1\) to \(N\), e.g. `{1..4}`

expands to `1 2 3 4`

.
Note that the expression `{1..N}`

cannot contain spaces.

To run 4 chains for the example bernoulli model on MacOS or Linux:

```
> for i in {1..4}
do
./bernoulli sample data file=bernoulli.data.json \
output file=output_${i}.csv
done
```

The backslash (`\`

) indicates a line continuation in Unix.
The expression `${i}`

substitutes in the value of loop index variable `i`

.
To run chains in parallel, put an ampersand (`&`

) at the end of the nested sampler command:

```
> for i in {1..4}
do
./bernoulli sample data file=bernoulli.data.json \
output file=output_${i}.csv &
done
```

This pushes each process into the background which allows the loop to continue without waiting for the current chain to finish.

On Windows, the DOS for-loop syntax is one of:

```
for %i in (SET) do COMMAND COMMAND-ARGUMENTS
for /l %i in (START, STEP, END) do COMMAND COMMAND-ARGUMENTS
```

To run 4 chains in parallel on Windows:

```
>for /l %i in (1, 1, 4) do start /b bernoulli.exe sample ^
data file=bernoulli.data.json my_data ^
output file=output_%i.csv
```

The caret (`^`

) indicates a line continuation in DOS.

## 4.3 Stan CSV Output File

Each execution of the model results in draws from a single Markov
chain being written to a file in comma-separated value (CSV) format.
The default name of the output file is `output.csv`

.

The first part of the output file records the version of the
underlying Stan library and the configuration as comments (i.e., lines
beginning with the pound sign (`#`

)).

```
# stan_version_major = 2
# stan_version_minor = 23
# stan_version_patch = 0
# model = bernoulli_model
# method = sample (Default)
# sample
# num_samples = 1000 (Default)
# num_warmup = 1000 (Default)
...
# output
# file = output.csv (Default)
# diagnostic_file = (Default)
# refresh = 100 (Default)
```

This is followed by a CSV header indicating the names of the values sampled.

`lp__,accept_stat__,stepsize__,treedepth__,n_leapfrog__,divergent__,energy__,theta`

The first output columns report the HMC sampler information:

`lp__`

- the total log probability density (up to an additive constant) at each sample`accept_stat__`

- the average Metropolis acceptance probability over each simulated Hamiltonian trajectory`stepsize__`

- integrator step size`treedepth__`

- depth of tree used by NUTS (NUTS sampler)`n_leapfrog__`

- number of leapfrog calculations (NUTS sampler)`divergent__`

- has value`1`

if trajectory diverged, otherwise`0`

. (NUTS sampler)`energy__`

- value of the Hamiltonian`int_time__`

- total integration time (static HMC sampler)

Because the above header is from the NUTS sampler, it has columns `treedepth__`

, `n_leapfrog__`

, and `divergent__`

and doesn’t have column `int_time__`

.
The remaining columns correspond to model parameters. For the
Bernoulli model, it is just the final column, `theta`

.

The header line is written to the output file before warmup begins.
If option `save_warmup`

is set to `1`

, the warmup draws are output directly after the header.
The total number of warmup draws saved is `num_warmup`

divided by `thin`

, rounded up (i.e., `ceiling`

).

Following the warmup draws (if any), are comments which record the results of adaptation: the stepsize, and inverse mass metric used during sampling:

```
# Adaptation terminated
# Step size = 0.884484
# Diagonal elements of inverse mass matrix:
# 0.535006
```

The default sampler is NUTS with an adapted step size and a diagonal
inverse mass matrix. For this example, the step size is 0.884484, and
the inverse mass contains the single entry 0.535006 corresponding to
the parameter `theta`

.

Draws from the posterior distribution are printed out next, each line containing a single draw with the columns corresponding to the header.

```
-6.84097,0.974135,0.884484,1,3,0,6.89299,0.198853
-6.91767,0.985167,0.884484,1,1,0,6.92236,0.182295
-7.04879,0.976609,0.884484,1,1,0,7.05641,0.162299
-6.88712,1,0.884484,1,1,0,7.02101,0.188229
-7.22917,0.899446,0.884484,1,3,0,7.73663,0.383596
...
```

The output ends with timing details:

```
# Elapsed Time: 0.007 seconds (Warm-up)
# 0.017 seconds (Sampling)
# 0.024 seconds (Total)
```

## 4.4 Summarizing Sampler Output(s) with `stansummary`

The `stansummary`

utility processes one or more output files from a run
or set of runs of Stan’s HMC sampler given a model and data.
For all columns in the Stan CSV output file `stansummary`

reports a set of statistics
including mean, standard deviation, percentiles, effective number of samples, and \(\hat{R}\) values.

To run `stansummary`

on the output files generated by the for loop above,
by the above run of the `bernoulli`

model on Mac or Linux:

`<cmdstan-home>/bin/stansummary output_*.csv`

On Windows, use backslashes to call the `stansummary.exe`

.

`<cmdstan-home>\bin\stansummary.exe output_*.csv`

The stansummary output consists of one row of statistics per column
in the Stan CSV output file. Therefore, the first rows in the
stansummary report statistics over the sampler state.
The final row of output summarizes the estimates of the model variable `theta`

:

```
Inference for Stan model: bernoulli_model
4 chains: each with iter=(1000,1000,1000,1000); warmup=(0,0,0,0); thin=(1,1,1,1); 4000 iterations saved.
Warmup took (0.0070, 0.0070, 0.0070, 0.0070) seconds, 0.028 seconds total
Sampling took (0.020, 0.017, 0.021, 0.019) seconds, 0.077 seconds total
Mean MCSE StdDev 5% 50% 95% N_Eff N_Eff/s R_hat
lp__ -7.3 1.8e-02 0.75 -8.8 -7.0 -6.8 1.8e+03 2.4e+04 1.0e+00
accept_stat__ 0.89 2.7e-03 0.17 0.52 0.96 1.0 3.9e+03 5.1e+04 1.0e+00
stepsize__ 1.1 7.5e-02 0.11 0.93 1.2 1.2 2.0e+00 2.6e+01 2.5e+13
treedepth__ 1.4 8.1e-03 0.49 1.0 1.0 2.0 3.6e+03 4.7e+04 1.0e+00
n_leapfrog__ 2.3 1.7e-02 0.98 1.0 3.0 3.0 3.3e+03 4.3e+04 1.0e+00
divergent__ 0.00 nan 0.00 0.00 0.00 0.00 nan nan nan
energy__ 7.8 2.6e-02 1.0 6.8 7.5 9.9 1.7e+03 2.2e+04 1.0e+00
theta 0.25 2.9e-03 0.12 0.079 0.23 0.46 1.7e+03 2.1e+04 1.0e+00
Samples were drawn using hmc with nuts.
For each parameter, N_Eff is a crude measure of effective sample size,
and R_hat is the potential scale reduction factor on split chains (at
convergence, R_hat=1).
```

In this example, we conditioned the model on a dataset consisting of the outcomes of
10 bernoulli trials, where only 2 trials reported success. The 5%, 50%, and 95%
percentile values for `theta`

reflect the uncertainty in our estimate, due to the
small amount of data, given the prior of `beta(1, 1)`