cmdstanpy_tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CmdStanPy Tutorial: do inference over a model conditioned on data\n",
    "\n",
    "\n",
    "## Overview\n",
    "\n",
    "The Bayesian workflow provides a framework for model comparison and model expansion.  To facilitate this workflow, CmdStanPy provides a set of objects which manage the model, and the result of doing inference over a model conditioned on data.   The building blocks of a Bayesian data analysis are\n",
    "\n",
    " + the model, written as a Stan program\n",
    " + the data, i.e., the definitions for all data variables in the Stan program\n",
    " + the resulting inferences over the model conditioned on the data:\n",
    "     + the inference engine configuration, i.e., all arguments to the CmdStan command-line interface\n",
    "     + the inference engine outputs:  output and diagnostic files, console message, error messages, and process return codes.\n",
    "\n",
    "#### The model\n",
    "\n",
    "The CmdStanPy class `CmdStanModel` encapsulates the model.  When working in CmdStanPy, all models must be supplied either as a Stan program file or as a compiled model executable file; or both, in which case, the base filenames must match.   While writing a program as a file requires a text editor, as opposed to the convenience of a Python notebook, in the long run, this is a more efficient way to develop complex models. \n",
    "\n",
    "#### The data\n",
    "\n",
    "Definitions of all data variables in the model are supplied either as a Python dictionary or as an on-disk data file which can be in either JSON or Rdump format.  Python dictionaries are converted to JSON inputs to CmdStanPy.  It is up to the user to save and manage the data; the data itself is not a first-class object in CmdStanPy.\n",
    "\n",
    "\n",
    "#### Inference Methods\n",
    "\n",
    "The `CmdStanModel` provides a set of inference methods, each of which returns an inference-specific result class which encapsulates the appropriate information and operations:\n",
    "\n",
    " + `sample` - Run the HMC-NUTS sampler to produce a set of draws from the posterior distribution. Returns a `CmdStanMCMC` object which contains all draws from all sampler chains.\n",
    " + `optimize` - Produce a penalized maximum likelihood estimate (point estimate) of the model parameters.  Returns a `CmdStanMLE` object containing point estimates for all model variables.\n",
    " + `variational` - Run CmdStan’s variational inference algorithm to approximate the posterior distribution. Returns a `CmdStanVB` object containing the approximate estimates.\n",
    " + `generate_quantities` - Run CmdStan’s generate_quantities method to produce additional quantities of interest based on each draw from an existing sample. Returns a `CmdStanGQ` object which contains the draws for the generated quantities variables.\n",
    "\n",
    "\n",
    "#### Output files and logs\n",
    "\n",
    "The CmdStan interface is entirely file-based, therefore CmdStanPy manages the process of converting in-memory data structures to their file-based representation and vice versa.\n",
    "\n",
    "Each run of an inference algorithm produces one or more output files in [Stan CSV format](https://mc-stan.org/docs/cmdstan-guide/stan-csv.html), as well as the messages written to stdout and stderr.  CmdStanPy uses per-run directories to keep all files together; the location of the directory can be specified when setting up the run, or it can be specified and/or changed afterwards via the `save_csvfiles` method, common to all inference method classes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Installation\n",
    "\n",
    "CmdStanPy is a pure-Python3 package, but it relies on CmdStan for all\n",
    "of its functionality. There are several ways to install CmdStan and CmdStanPy, \n",
    "which depend on the kind of user you are. \n",
    "\n",
    "\n",
    "#### Conda users (Recommended)\n",
    "\n",
    "If you use [conda](https://docs.conda.io/en/latest/),\n",
    "installation of both can be done very simply. CmdStanPy\n",
    "and CmdStan are both available via the \n",
    "[conda-forge](https://conda-forge.org/) repository.\n",
    "\n",
    "We recommend creating a new environment for CmdStan[Py]:\n",
    "\n",
    "```\n",
    "conda create -n cmdstan -c conda-forge cmdstanpy\n",
    "```\n",
    "\n",
    "but installation is possible in an existing environment:\n",
    "\n",
    "```\n",
    "conda install -c conda-forge cmdstanpy\n",
    "```\n",
    "\n",
    "These commands will install CmdStanPy, CmdStan, and the\n",
    "required compilers for using CmdStan on your system inside\n",
    "a conda environment. To use them, run `conda activate cmdstan`,\n",
    "or whichever name you used for your environment (following `-n`\n",
    "above).\n",
    "\n",
    "Note that CmdStan is only available on conda for versions\n",
    "2.27.0 and newer. If you require an older version, you must use\n",
    "one of the following methods to install it. If you require a\n",
    "version of CmdStan *newer* than 2.27.0, but not the latest,\n",
    "you can install it in the standard conda way by specifying\n",
    "`cmdstan==VERSION` in the install command.\n",
    "\n",
    "#### Pip (non-Conda) users \n",
    "\n",
    "CmdStan can also be installed from PyPI via URL: https://pypi.org/project/cmdstanpy/ or from the\n",
    "command line using `pip`:\n",
    "\n",
    "```\n",
    "pip install --upgrade cmdstanpy\n",
    "```\n",
    "\n",
    "CmdStanPy will use the following optional packages, if installed:\n",
    "\n",
    "* `xarray`, an n-dimension labeled dataset package which can be used for outputs\n",
    "\n",
    "To install CmdStanPy with all the optional packages:\n",
    "\n",
    "```\n",
    "pip install --upgrade cmdstanpy[all]\n",
    "```\n",
    "\n",
    "To install the current develop branch from GitHub:\n",
    "\n",
    "```\n",
    "pip install -e git+https://github.com/stan-dev/cmdstanpy@/develop#egg=cmdstanpy\n",
    "```\n",
    "\n",
    "If you install CmdStanPy from GitHub,\n",
    "**you must install CmdStan**. The recommended way for Pip users \n",
    "to do so is via the `install_cmdstan` function. \n",
    "\n",
    "*Note for PyStan & RTools users:*  PyStan and CmdStanPy should be installed in \n",
    "separate environments if you are using the RTools toolchain (primarily Windows users).\n",
    "If you already have PyStan installed, you should take care to install CmdStanPy in its own\n",
    "virtual environment.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hello World:  Compile and run CmdStan example model `bernoulli.stan`\n",
    "\n",
    "Every CmdStan release has an `examples/bernoulli` directory which contains the Stan model and test data files.  In this example we compile the model and run the NUTS-HMC sampler to obtain the full posterior density.\n",
    "\n",
    "The `CmdStanModel` class method `sample` invokes Stan's NUTS-HMC sampler which conditions the model on the input data and returns a `CmdStanMCMC` object. The `CmdStanMCMC` object contains the set of draws from the posterior plus metadata."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import CmdStanPy classes and methods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "from cmdstanpy import cmdstan_path, CmdStanModel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compile the model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The CmdStanPy directory `examples/bernoulli` contains the model and data files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bernoulli_dir = os.path.join(cmdstan_path(), 'examples', 'bernoulli')\n",
    "stan_file = os.path.join(bernoulli_dir, 'bernoulli.stan')\n",
    "with open(stan_file, 'r') as f:\n",
    "    print(f.read())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `CmdStanModel` object stores the filepaths of the Stan program and the compiled executable file.  By default, if no executable file is supplied, or if the timestamp on the model file is newer than the timestamp on the executable file, the model will be compiled on instantiation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = CmdStanModel(stan_file=stan_file)\n",
    "print(model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Assemble the data\n",
    "\n",
    "Input data is either a Python `Dict` with entries corresponding to input data values, or it can be a file in JSON or Rdump format.\n",
    "\n",
    "If a `Dict` is specified, CmdStanPy writes it to a temp file in JSON format."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = { \"N\" : 10, \"y\" : [0,1,0,0,0,0,0,0,0,1] }"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the CmdStan `examples/bernoulli` directory, there are data files in both `JSON` and `rdump` formats."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_file = os.path.join(bernoulli_dir, 'bernoulli.data.json')\n",
    "with open(data_file, 'r') as f:\n",
    "    print(f.read())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Do full Baysian inference using the HMC-NUTS sampler\n",
    "\n",
    "The `CmdStanModel` class `sample` method runs the HMC-NUTS sampler to produce a set of draws from the posterior distribution.  By default, `sample` runs 4 sampler chains."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fit = model.sample(data=data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `sample` method returns a `CmdStanMCMC` object.  This object records all sampler configuration information and the location of all outputs produced."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(fit)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The CmdStanMCMC class has properties and methods to access the sampler outputs and all metadata.  Accessor properties include:\n",
    "\n",
    "- `draws`  - all draws from all chains, stored as a 3-D numpy.ndarray.\n",
    "- `num_draws` - total number of warmup and sampling draws\n",
    "- `chains` - number of chains run by sampler\n",
    "- `metric`, `step_size` - per-chain HMC tuning parameters\n",
    "- `csv_files` - list of Stan csv output files which comprise the sample\n",
    "    \n",
    "The `draws` array is created only as needed; therefore the first time that this property is accessed, CmdStanPy parses the set of csv output files, at the same time it parses the contain step size and metric information.  Depending on the size of the sample and your computing environment, this may take a few seconds or more.\n",
    "\n",
    "The numpy ndarray is stored column major format so that values for each parameter are stored contiguously in memory, likewise all draws from a chain are contiguous.  Thus the dimensions of the ndarray are arranged as follows:  (draws, chains, columns).  The draws array contains both the sampler variables and the model variables.  Sampler variables report the sampler state.  For the NUTS-HMC sampler, there are 7 reported variables, all of which end in `__`.  The example model contains a single variable `theta`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fit.draws().shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Viewing the sample as a set of variables:  sampler variables and program variables\n",
    "\n",
    "The `CmdStanMCMC` property `draws` is a 3-D numpy ndarray which contains all draws across all chains.  \n",
    "    \n",
    "- `stan_variables` - dictionary mapping Stan program variables names to the draws for that variable.\n",
    "- `sampler_variables` - dictionary mapping the sampler variables names, (e.g., `lp__`) to the draws for that variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vars = fit.stan_variables()\n",
    "for (k,v) in vars.items():\n",
    "    print(k, v.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Get HMC sampler tuning parameters\n",
    "\n",
    "##### step_size\n",
    "\n",
    "The `CmdStanMCMC` property `step_size` property is a 1-D numpy ndarray which contains the step size used by the sampler for each chain.  This array is created at the same time as the `sample` and `metric` arrays are created.\n",
    "\n",
    "At the end of adaptation, the step size for the 4 chains in this example is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fit.step_size"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### metric_type, metric\n",
    "\n",
    "The `metric` property is an numpy ndarray which contains the metric used by the sampler for each chain.  This array is created at the same time as the `sample` and `step_size` arrays are created.\n",
    "\n",
    "At the end of adaptation, the metric for the 4 chains in this example is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fit.metric_type,  fit.metric"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Summarize the results\n",
    "\n",
    "The `CmdStanMCMC` method `summary` runs the CmdStan `bin/stansummary` over the fitted sample csv files.  The output of the command is streamed to the console.  The method returns a pandas.DataFrame which contains\n",
    "estimates of posterior means, standard deviations, Monte-Carlo standard error, effective sample size, and convergence diagnostic statistic for all parameters, transformed parameters, and generated quantities variables:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(fit.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Run sampler diagnostics\n",
    "\n",
    "The `CmdStanMCMC` method `diagnose` prints the output of CmdStan bin/stansummary to the console.  This is useful when the summary report shows `R_hat` values outside of the range of .99 to 1.01, or the number of effective samples (`N_eff` value) is below 5% of the total draws in the sample.   It returns a string which contains the console output.\n",
    "\n",
    "\n",
    "   + The `CmdStanMCMC` class method `diagnose` invokes CmdStan's `diagnose` utility which checks for the following problems:\n",
    "    + transitions that hit the maximum treedepth\n",
    "    + divergent transitions\n",
    "    + low E-BFMI values (sampler transitions HMC potential energy)\n",
    "    + low effective sample sizes\n",
    "    + high R-hat values\n",
    "   + See the Stan reference manual section on [posterior analysis](https://mc-stan.org/docs/reference-manual/analysis-chapter.html) for further details.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(fit.diagnose())"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "d31ce8e45781476cfd394e192e0962028add96ff436d4fd4e560a347d206b9cb"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}