# MCMC Sampling¶

The `CmdStanModel`

class method `sample`

invokes Stan’s adaptive HMC-NUTS sampler which uses the Hamiltonian Monte Carlo (HMC) algorithm and its adaptive variant the no-U-turn sampler (NUTS) to produce a set of draws from the posterior distribution of the model parameters conditioned on the data. It returns a `CmdStanMCMC`

object which provides properties to retrieve information about the sample, as well as methods to run CmdStan’s summary and diagnostics tools.

In order to evaluate the fit of the model to the data, it is necessary to run several Monte Carlo chains and compare the set of draws returned by each. By default, the `sample`

command runs 4 sampler chains, i.e., CmdStanPy invokes CmdStan 4 times. CmdStanPy uses Python’s `subprocess`

and `multiprocessing`

libraries to run these chains in separate processes. This processing can be done in parallel, up to the number of processor cores available.

## Prerequisites¶

CmdStanPy displays progress bars during sampling via use of package tqdm. In order for these to display properly, you must have the ipywidgets package installed, and depending on your version of Jupyter or JupyterLab, you must enable it via command:

```
[1]:
```

```
!jupyter nbextension enable --py widgetsnbextension
```

```
Enabling notebook extension jupyter-js-widgets/extension...
- Validating: OK
```

For more information, see the the installation instructions, also this tqdm GitHub issue.

## Fitting a model to data¶

In this example we use the CmdStan example model bernoulli.stan and data file bernoulli.data.json.

We instantiate a `CmdStanModel`

from the Stan program file

```
[2]:
```

```
import os
from cmdstanpy import CmdStanModel, cmdstan_path
bernoulli_dir = os.path.join(cmdstan_path(), 'examples', 'bernoulli')
stan_file = os.path.join(bernoulli_dir, 'bernoulli.stan')
data_file = os.path.join(bernoulli_dir, 'bernoulli.data.json')
# instantiate, compile bernoulli model
model = CmdStanModel(stan_file=stan_file)
```

```
INFO:cmdstanpy:compiling stan file /home/docs/checkouts/readthedocs.org/user_builds/cmdstanpy/conda/v1.0.0rc3/bin/cmdstan/examples/bernoulli/bernoulli.stan to exe file /home/docs/checkouts/readthedocs.org/user_builds/cmdstanpy/conda/v1.0.0rc3/bin/cmdstan/examples/bernoulli/bernoulli
INFO:cmdstanpy:compiled model executable: /home/docs/checkouts/readthedocs.org/user_builds/cmdstanpy/conda/v1.0.0rc3/bin/cmdstan/examples/bernoulli/bernoulli
```

By default, the model is compiled during instantiation. The compiled executable is created in the same directory as the program file. If the directory already contains an executable file with a newer timestamp, the model is not recompiled.

We run the sampler on the data using all default settings: 4 chains, each of which runs 1000 warmup and sampling iterations.

```
[3]:
```

```
# run CmdStan's sample method, returns object `CmdStanMCMC`
fit = model.sample(data=data_file)
```

```
INFO:cmdstanpy:CmdStan start procesing
```

```
```

```
INFO:cmdstanpy:CmdStan done processing.
```

```
```

The `sample`

method returns a `CmdStanMCMC`

object, which contains: - metadata - draws - HMC tuning parameters `metric`

, `step_size`

```
[4]:
```

```
print('sampler diagnostic variables:\n{}'.format(fit.metadata.method_vars_cols.keys()))
print('stan model variables:\n{}'.format(fit.metadata.stan_vars_cols.keys()))
```

```
sampler diagnostic variables:
dict_keys(['lp__', 'accept_stat__', 'stepsize__', 'treedepth__', 'n_leapfrog__', 'divergent__', 'energy__'])
stan model variables:
dict_keys(['theta'])
```

```
[5]:
```

```
fit.summary()
```

```
[5]:
```

Mean | MCSE | StdDev | 5% | 50% | 95% | N_Eff | N_Eff/s | R_hat | |
---|---|---|---|---|---|---|---|---|---|

name | |||||||||

lp__ | -7.30 | 0.0210 | 0.74 | -8.800 | -7.00 | -6.70 | 1300.0 | 14000.0 | 1.0 |

theta | 0.25 | 0.0031 | 0.12 | 0.074 | 0.23 | 0.46 | 1400.0 | 16000.0 | 1.0 |

The sampling data from the fit can be accessed either as a `numpy`

array or a pandas `DataFrame`

:

```
[6]:
```

```
print(fit.draws().shape)
fit.draws_pd().head()
```

```
(1000, 4, 8)
```

```
[6]:
```

lp__ | accept_stat__ | stepsize__ | treedepth__ | n_leapfrog__ | divergent__ | energy__ | theta | |
---|---|---|---|---|---|---|---|---|

0 | -8.18109 | 0.947147 | 1.08292 | 2.0 | 3.0 | 0.0 | 8.18115 | 0.085541 |

1 | -9.25583 | 0.960750 | 1.08292 | 2.0 | 3.0 | 0.0 | 9.28633 | 0.054002 |

2 | -8.48709 | 1.000000 | 1.08292 | 1.0 | 1.0 | 0.0 | 9.27267 | 0.074519 |

3 | -8.60606 | 0.986149 | 1.08292 | 1.0 | 1.0 | 0.0 | 8.88402 | 0.070754 |

4 | -7.83236 | 1.000000 | 1.08292 | 1.0 | 1.0 | 0.0 | 8.54461 | 0.101193 |

Additionally, if `xarray`

is installed, this data can be accessed another way:

```
[7]:
```

```
fit.draws_xr()
```

```
[7]:
```

<xarray.Dataset> Dimensions: (draw: 1000, chain: 4) Coordinates: * chain (chain) int64 1 2 3 4 * draw (draw) int64 0 1 2 3 4 5 6 7 8 ... 992 993 994 995 996 997 998 999 Data variables: theta (chain, draw) float64 0.08554 0.054 0.07452 ... 0.1688 0.2943 Attributes: stan_version: 2.28.1 model: bernoulli_model num_draws_sampling: 1000

- draw: 1000
- chain: 4

- chain(chain)int641 2 3 4
array([1, 2, 3, 4])

- draw(draw)int640 1 2 3 4 5 ... 995 996 997 998 999
array([ 0, 1, 2, ..., 997, 998, 999])

- theta(chain, draw)float640.08554 0.054 ... 0.1688 0.2943
array([[0.0855405, 0.0540022, 0.0745186, ..., 0.357497 , 0.335188 , 0.229181 ], [0.175374 , 0.297097 , 0.164462 , ..., 0.121707 , 0.253696 , 0.233111 ], [0.154479 , 0.177931 , 0.228719 , ..., 0.171627 , 0.134446 , 0.265674 ], [0.149893 , 0.227655 , 0.237765 , ..., 0.246552 , 0.168819 , 0.294343 ]])

- stan_version :
- 2.28.1
- model :
- bernoulli_model
- num_draws_sampling :
- 1000

The `fit`

object records the command, the return code, and the paths to the sampler output csv and console files. The string representation of this object displays the CmdStan commands and the location of the output files.

Output filenames are composed of the model name, a timestamp in the form YYYYMMDDhhmm and the chain id, plus the corresponding filetype suffix, either ‘.csv’ for the CmdStan output or ‘.txt’ for the console messages, e.g. `bernoulli-201912081451-1.csv`

. Output files written to the temporary directory contain an additional 8-character random string, e.g. `bernoulli-201912081451-1-5nm6as7u.csv`

.

```
[8]:
```

```
fit
```

```
[8]:
```

```
CmdStanMCMC: model=bernoulli chains=4['method=sample', 'algorithm=hmc', 'adapt', 'engaged=1']
csv_files:
/tmp/tmpba1_ahip/bernoulli-20211119014442_1.csv
/tmp/tmpba1_ahip/bernoulli-20211119014442_2.csv
/tmp/tmpba1_ahip/bernoulli-20211119014442_3.csv
/tmp/tmpba1_ahip/bernoulli-20211119014442_4.csv
output_files:
/tmp/tmpba1_ahip/bernoulli-20211119014442_0-stdout.txt
/tmp/tmpba1_ahip/bernoulli-20211119014442_1-stdout.txt
/tmp/tmpba1_ahip/bernoulli-20211119014442_2-stdout.txt
/tmp/tmpba1_ahip/bernoulli-20211119014442_3-stdout.txt
```

The sampler output files are written to a temporary directory which is deleted upon session exit unless the `output_dir`

argument is specified. The `save_csvfiles`

function moves the CmdStan CSV output files to a specified directory without having to re-run the sampler. The console output files are not saved. These files are treated as ephemeral; if the sample is valid, all relevant information is recorded in the CSV files.

### Sampler Progress¶

Your model make take a long time to fit. The `sample`

method provides two arguments:

visual progress bar:

`show_progress=True`

stream CmdStan ouput to the console -

`show_console=True`

To illustrate how progress bars work, we will run the bernoulli model. Since the progress bars are only visible while the sampler is running and the bernoulli model takes no time at all to fit, we run this model for 200K iterations, in order to see the progress bars in action.

```
[9]:
```

```
fit = model.sample(data=data_file, iter_warmup=100000, iter_sampling=100000, show_progress=True)
```

```
INFO:cmdstanpy:CmdStan start procesing
```

```
```

```
INFO:cmdstanpy:CmdStan done processing.
```

```
```

The Stan language `print`

statement can be use to monitor the Stan program state. In order to see this information as the sampler is running, use the `show_console=True`

argument. This will stream all CmdStan messages to the terminal while the sampler is running.

```
[10]:
```

```
fit = model.sample(data=data_file, chains=2, parallel_chains=1, show_console=True)
```

```
INFO:cmdstanpy:Chain [1] start processing
INFO:cmdstanpy:Chain [1] done processing
INFO:cmdstanpy:Chain [2] start processing
INFO:cmdstanpy:Chain [2] done processing
```

```
Chain [1] method = sample (Default)
Chain [1] sample
Chain [1] num_samples = 1000 (Default)
Chain [1] num_warmup = 1000 (Default)
Chain [1] save_warmup = 0 (Default)
Chain [1] thin = 1 (Default)
Chain [1] adapt
Chain [1] engaged = 1 (Default)
Chain [1] gamma = 0.050000000000000003 (Default)
Chain [1] delta = 0.80000000000000004 (Default)
Chain [1] kappa = 0.75 (Default)
Chain [1] t0 = 10 (Default)
Chain [1] init_buffer = 75 (Default)
Chain [1] term_buffer = 50 (Default)
Chain [1] window = 25 (Default)
Chain [1] algorithm = hmc (Default)
Chain [1] hmc
Chain [1] engine = nuts (Default)
Chain [1] nuts
Chain [1] max_depth = 10 (Default)
Chain [1] metric = diag_e (Default)
Chain [1] metric_file = (Default)
Chain [1] stepsize = 1 (Default)
Chain [1] stepsize_jitter = 0 (Default)
Chain [1] num_chains = 1 (Default)
Chain [1] id = 1 (Default)
Chain [1] data
Chain [1] file = /home/docs/checkouts/readthedocs.org/user_builds/cmdstanpy/conda/v1.0.0rc3/bin/cmdstan/examples/bernoulli/bernoulli.data.json
Chain [1] init = 2 (Default)
Chain [1] random
Chain [1] seed = 50018
Chain [1] output
Chain [1] file = /tmp/tmpba1_ahip/bernoulli-20211119014451_1.csv
Chain [1] diagnostic_file = (Default)
Chain [1] refresh = 100 (Default)
Chain [1] sig_figs = -1 (Default)
Chain [1] profile_file = profile.csv (Default)
Chain [1] num_threads = 1 (Default)
Chain [1]
Chain [1]
Chain [1] Gradient evaluation took 5e-06 seconds
Chain [1] 1000 transitions using 10 leapfrog steps per transition would take 0.05 seconds.
Chain [1] Adjust your expectations accordingly!
Chain [1]
Chain [1]
Chain [1] Iteration: 1 / 2000 [ 0%] (Warmup)
Chain [1] Iteration: 100 / 2000 [ 5%] (Warmup)
Chain [1] Iteration: 200 / 2000 [ 10%] (Warmup)
Chain [1] Iteration: 300 / 2000 [ 15%] (Warmup)
Chain [1] Iteration: 400 / 2000 [ 20%] (Warmup)
Chain [1] Iteration: 500 / 2000 [ 25%] (Warmup)
Chain [1] Iteration: 600 / 2000 [ 30%] (Warmup)
Chain [1] Iteration: 700 / 2000 [ 35%] (Warmup)
Chain [1] Iteration: 800 / 2000 [ 40%] (Warmup)
Chain [1] Iteration: 900 / 2000 [ 45%] (Warmup)
Chain [1] Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain [1] Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain [1] Iteration: 1100 / 2000 [ 55%] (Sampling)
Chain [1] Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain [1] Iteration: 1300 / 2000 [ 65%] (Sampling)
Chain [1] Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain [1] Iteration: 1500 / 2000 [ 75%] (Sampling)
Chain [1] Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain [1] Iteration: 1700 / 2000 [ 85%] (Sampling)
Chain [1] Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain [1] Iteration: 1900 / 2000 [ 95%] (Sampling)
Chain [1] Iteration: 2000 / 2000 [100%] (Sampling)
Chain [1]
Chain [1] Elapsed Time: 0.008 seconds (Warm-up)
Chain [1] 0.013 seconds (Sampling)
Chain [1] 0.021 seconds (Total)
Chain [1]
Chain [1]
Chain [2] method = sample (Default)
Chain [2] sample
Chain [2] num_samples = 1000 (Default)
Chain [2] num_warmup = 1000 (Default)
Chain [2] save_warmup = 0 (Default)
Chain [2] thin = 1 (Default)
Chain [2] adapt
Chain [2] engaged = 1 (Default)
Chain [2] gamma = 0.050000000000000003 (Default)
Chain [2] delta = 0.80000000000000004 (Default)
Chain [2] kappa = 0.75 (Default)
Chain [2] t0 = 10 (Default)
Chain [2] init_buffer = 75 (Default)
Chain [2] term_buffer = 50 (Default)
Chain [2] window = 25 (Default)
Chain [2] algorithm = hmc (Default)
Chain [2] hmc
Chain [2] engine = nuts (Default)
Chain [2] nuts
Chain [2] max_depth = 10 (Default)
Chain [2] metric = diag_e (Default)
Chain [2] metric_file = (Default)
Chain [2] stepsize = 1 (Default)
Chain [2] stepsize_jitter = 0 (Default)
Chain [2] num_chains = 1 (Default)
Chain [2] id = 2
Chain [2] data
Chain [2] file = /home/docs/checkouts/readthedocs.org/user_builds/cmdstanpy/conda/v1.0.0rc3/bin/cmdstan/examples/bernoulli/bernoulli.data.json
Chain [2] init = 2 (Default)
Chain [2] random
Chain [2] seed = 50018
Chain [2] output
Chain [2] file = /tmp/tmpba1_ahip/bernoulli-20211119014451_2.csv
Chain [2] diagnostic_file = (Default)
Chain [2] refresh = 100 (Default)
Chain [2] sig_figs = -1 (Default)
Chain [2] profile_file = profile.csv (Default)
Chain [2] num_threads = 1 (Default)
Chain [2]
Chain [2]
Chain [2] Gradient evaluation took 4e-06 seconds
Chain [2] 1000 transitions using 10 leapfrog steps per transition would take 0.04 seconds.
Chain [2] Adjust your expectations accordingly!
Chain [2]
Chain [2]
Chain [2] Iteration: 1 / 2000 [ 0%] (Warmup)
Chain [2] Iteration: 100 / 2000 [ 5%] (Warmup)
Chain [2] Iteration: 200 / 2000 [ 10%] (Warmup)
Chain [2] Iteration: 300 / 2000 [ 15%] (Warmup)
Chain [2] Iteration: 400 / 2000 [ 20%] (Warmup)
Chain [2] Iteration: 500 / 2000 [ 25%] (Warmup)
Chain [2] Iteration: 600 / 2000 [ 30%] (Warmup)
Chain [2] Iteration: 700 / 2000 [ 35%] (Warmup)
Chain [2] Iteration: 800 / 2000 [ 40%] (Warmup)
Chain [2] Iteration: 900 / 2000 [ 45%] (Warmup)
Chain [2] Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain [2] Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain [2] Iteration: 1100 / 2000 [ 55%] (Sampling)
Chain [2] Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain [2] Iteration: 1300 / 2000 [ 65%] (Sampling)
Chain [2] Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain [2] Iteration: 1500 / 2000 [ 75%] (Sampling)
Chain [2] Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain [2] Iteration: 1700 / 2000 [ 85%] (Sampling)
Chain [2] Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain [2] Iteration: 1900 / 2000 [ 95%] (Sampling)
Chain [2] Iteration: 2000 / 2000 [100%] (Sampling)
Chain [2]
Chain [2] Elapsed Time: 0.007 seconds (Warm-up)
Chain [2] 0.014 seconds (Sampling)
Chain [2] 0.021 seconds (Total)
Chain [2]
Chain [2]
```

## Running a data-generating model¶

In this example we use the CmdStan example model data_filegen.stan to generate a simulated dataset given fixed data values.

```
[11]:
```

```
model_datagen = CmdStanModel(stan_file='bernoulli_datagen.stan')
datagen_data = {'N':300, 'theta':0.3}
fit_sim = model_datagen.sample(data=datagen_data, fixed_param=True)
fit_sim.summary()
```

```
INFO:cmdstanpy:compiling stan file /home/docs/checkouts/readthedocs.org/user_builds/cmdstanpy/checkouts/v1.0.0rc3/docsrc/examples/bernoulli_datagen.stan to exe file /home/docs/checkouts/readthedocs.org/user_builds/cmdstanpy/checkouts/v1.0.0rc3/docsrc/examples/bernoulli_datagen
INFO:cmdstanpy:compiled model executable: /home/docs/checkouts/readthedocs.org/user_builds/cmdstanpy/checkouts/v1.0.0rc3/docsrc/examples/bernoulli_datagen
INFO:cmdstanpy:CmdStan start procesing
```

```
```

```
INFO:cmdstanpy:CmdStan done processing.
```

```
```

```
[11]:
```

Mean | MCSE | StdDev | 5% | 50% | 95% | N_Eff | N_Eff/s | R_hat | |
---|---|---|---|---|---|---|---|---|---|

name | |||||||||

lp__ | 0 | NaN | 0.0 | 0 | 0 | 0.0 | NaN | NaN | NaN |

theta_rep | 90 | 0.26 | 8.1 | 76 | 90 | 100.0 | 970.0 | 140000.0 | 1.0 |

Compute, plot histogram of total successes for `N`

Bernoulli trials with chance of success `theta`

:

```
[12]:
```

```
drawset_pd = fit_sim.draws_pd()
drawset_pd.columns
# restrict to columns over new outcomes of N Bernoulli trials
y_sims = drawset_pd.drop(columns=['lp__', 'accept_stat__'])
# plot total number of successes per draw
y_sums = y_sims.sum(axis=1)
y_sums.astype('int32').plot.hist(range(0,datagen_data['N']+1))
```

```
[12]:
```

```
<AxesSubplot:ylabel='Frequency'>
```