source/analyses/fitting/MCMC_for_fitting_models_Workbook.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fitting a model with Markov Chain Monte Carlo\n",
    "\n",
    "Markov Chain Monte Carlo (MCMC) is a way to infer a distribution of model parameters, given that the measurements of the output of the model are influenced by some tractable random process. In this case, performs something akin to the opposite of what a standard Monte Carlo simulation will do. Instead of starting with distributions for the parameters of a model and using them to calculate a distribution (usually related to uncertainty) in the output of the simulation, we start with a distribution of that output and look for input distributions.\n",
    "\n",
    " \n",
    "## Ingredients\n",
    "\n",
    "For this analysis we'll introduce the python package [PyMC](https://pymc-devs.github.io/pymc/README.html) which implements MCMC algorithms for us. Another project which performs similar calculations is [PyStan](https://pystan.readthedocs.org/en/latest/)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%pylab inline\n",
    "import pysd\n",
    "import <<...>>\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this example, we'll revisit the ebola case, only assuming that our data has some noise, and we'll use MCMC to estimate distributions for the parameters for the model. For a more detailed description of this model and the dataset see the recipe [Fitting with Optimization](http://pysd-cookbook.readthedocs.org/en/latest/analyses/fitting/Fitting_with_Optimization.html).\n",
    "\n",
    "We'll assume that the model simulates an underlying process of disease propagation, but that the data is noisy - perhaps it represents admittance rates at a hospital, and so will be missing some cases, and may include some false positives through misdiagnosis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = pysd.read_vensim('../../models/SI_Model/SI_Model.mdl')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.read_csv('../../data/Ebola/Ebola_in_SL_Data.csv', index_col='Weeks')\n",
    "data.plot();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Recipe\n",
    "\n",
    "#### Step 1: Establish prior distributions for the input parameters\n",
    "Our first step is to establish prior distributions for the parameters for which we will be trying to infer posterior distributions. As in the optimization case, we'll modify the effective population and the contact frequency. In real epidemiological modeling, we'd have a more complex model, and we'd have some information about the population, etc. but this makes for a relatively tight example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "population = pymc.<<...>>('total_population', lower=2, upper=50000, value=10000)\n",
    "contact_frequency = pymc.<<..>>('contact_frequency', beta=5, value=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step 2: Perform the simulation\n",
    "We now need to define a function that will run the model given a draw from the input distributions, and return the values that we will take as those of the fundamental underlying phenomena, before measurement noise is included."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@pymc.<<...>>\n",
    "def sim(population=population, contact_frequency=contact_frequency):\n",
    "    result = model.run(params={'total_population':population,\n",
    "                               'contact_frequency':contact_frequency},\n",
    "                       return_columns=['infection_rate'],\n",
    "                       return_timestamps=list(data.index.values))\n",
    "    return result['infection_rate'].values\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step 3: Include noise terms\n",
    "There are several ways we could include noise. If we expected no false positives, we could use a Binomial distribution, such that of `n` possible cases that could be reported, only a fraction `p` would be reported, and other cases missed. If we only want to model false positives, we could assume that there was an average rate of false positives, with the data following a poisson distribution. The full rate would be the sum of these two processes.\n",
    "\n",
    "For now, however, we'll simplify the analysis by only looking at the Poisson noise component. The mean of the poisson process will be the results of our simulation.\n",
    "\n",
    "This is where we include our measured data into the model. PyMC will know how to calculate the log likelihood of seeing the observed data given the assumption that the simulation result represents the underlying process, subject to Poisson noise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "admittances = pymc.Poisson('admittances', mu=<<...>>,\n",
    "                           value=data['New Reported Cases'], observed=<<...>>)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step 4: Perform the MCMC Sampling\n",
    "\n",
    "Now that we have set up the problem for PyMC, we need only to run the MCMC sampler. What this will do, essentially, is take a trial set of points from our prior distribution, simulate the model, and evaluate the likelihood of the data given those input parameters, the simulation model, and the noise distribution. It will then use bayes law to decide whether to keep the trial points or throw them away. It will then choose a new set of points and start over. (There is a lot more cleverness happening than this, of course. If you want to know how it works, I recommend [Bayesian Methods for Hackers](http://camdavidsonpilon.github.io/Probabilistic-Programming-and-Bayesian-Methods-for-Hackers/).\n",
    "\n",
    "First we assemble the various pieces of the data flow that we built up into a model that pymc can recognize, and instantiate a sampler `MCMC` to run the algorithm for us.\n",
    "\n",
    "Then we'll ask the MCMC algorithm to run until it has kept 20000 points. We'll throw out the first 1000 of these, as they are likely to be biased towards the initial values we set up and not representative of the overall distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mcmdl = pymc.Model([population, contact_frequency, sim, admittances])\n",
    "mcmc = pymc.MCMC(<<...>>)\n",
    "mcmc.sample(20000,1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step 5: Look at the distribution\n",
    "\n",
    "We can now evaluate the results by looking at the series of points we 'kept'. These are stored as traces within the population and contact frequency objects we built earlier. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(6,2))\n",
    "plt.hist(population.<<...>>[:], bins=100, histtype='stepfilled', normed=True, alpha=.6);\n",
    "plt.title('Fit distribution for the population parameter');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(6,2))\n",
    "plt.hist(contact_frequency.trace[:], bins=100, histtype='stepfilled', normed=True, alpha=.6);\n",
    "plt.title('Fit distribution for the contact frequency parameter');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(6,1))\n",
    "plt.hist(contact_frequency.trace[:], bins=100, histtype='stepfilled', normed=True, alpha=.6);\n",
    "plt.yticks([])\n",
    "plt.title('Fit distribution for the characteristic time parameter');\n",
    "plt.xlabel('Minutes')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step 6: Understand the correlation between our parameter distributions\n",
    "PyMC has calculated for us the input distributions for the parameters we have based upon our data. This isn't the whole story, however, as these values are not independent of one another. We can see by plotting them against one another that there is correlation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(population.trace[:], contact_frequency.trace[:], '.', alpha=.1)\n",
    "plt.xlabel('Effective Population')\n",
    "plt.ylabel('Contact Frequency');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we use these distributions in the future to propagate uncertainty through the model (say in doing a policy test) we should make sure that we're including this correlation. The simplest way to do this is just to use the traces themselves as inputs into a normal Monte Carlo error propagation calculation, making sure that we use the same index in each trace for any given run. This will automatically include the correlation for us. Clever!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Resources:\n",
    "\n",
    "The analysis in this notebook is inspired by blogt posts by Abraham Flaxman:\n",
    "- http://healthyalgorithms.com/2010/11/29/mcmc-in-python-statistical-model-stuck-on-a-stochastic-system-dynamics-model-in-pymc/\n",
    "\n",
    "- http://healthyalgorithms.com/2010/10/19/mcmc-in-python-how-to-stick-a-statistical-model-on-a-system-dynamics-model-in-pymc/ \n",
    "\n",
    "A great resource for learning about MCMC is [Bayesian Methods for Hackers](http://camdavidsonpilon.github.io/Probabilistic-Programming-and-Bayesian-Methods-for-Hackers/)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.8.10 64-bit",
   "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.10"
  },
  "vscode": {
   "interpreter": {
    "hash": "916dbcbb3f70747c44a77c7bcd40155683ae19c65e1c03b4aa3499c5328201f1"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}