cs146-5.2-pre-class-work.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pre-class work\n",
    "Below is the data set from 6 medical trials on the effect of specific allergen immunotherapy (SIT) on eczema patients.\n",
    "\n",
    "| Study          | TG improved      | TG not improved   | CG improved    | CG not improved   |\n",
    "|:-------------- | --------:| ------:| ------:| ------:|\n",
    "| Di Rienzo 2014 | 20       | 3      | 9      | 6      |\n",
    "| Galli 1994     | 10       | 6      | 11     | 7      |\n",
    "| Kaufman 1974   | 13       | 3      | 4      | 6      |\n",
    "| Qin 2014       | 35       | 10     | 21     | 18     |\n",
    "| Sanchez 2012   | 22       | 9      | 12     | 17     |\n",
    "| Silny 2006     | 7        | 3      | 0      | 10     |\n",
    "| **Totals**     | **107**  | **34** | **57** | **64** |\n",
    "\n",
    "* TG = Treatment group\n",
    "* CG = Control group\n",
    "\n",
    "The model we used was that each trial's results were generated from a binomial distribution over the number of improved patients with a common improvement rate parameter shared between all trials.\n",
    "\n",
    "For the treatment group we use a subscript $t$:\n",
    "\n",
    "$$\\begin{align}\n",
    "k_{ti} &\\sim \\text{Binomial}(n_{ti}, p_t) \\qquad i=1,2,\\ldots 6\\\\\n",
    "p_t &\\sim \\text{Beta}(\\alpha=1, \\beta=1)\n",
    "\\end{align}$$\n",
    "\n",
    "For the control group we use a subscript $c$:\n",
    "\n",
    "$$\\begin{align}\n",
    "k_{ci} &\\sim \\text{Binomial}(n_{ci}, p_c) \\qquad i=1,2,\\ldots 6\\\\\n",
    "p_c &\\sim \\text{Beta}(\\alpha=1, \\beta=1)\n",
    "\\end{align}$$\n",
    "\n",
    "So we have the same model structure for the treatment and control groups, just with different data.\n",
    "\n",
    "The code below implements the Stan model for the scenario above.\n",
    "\n",
    "* Carefully **read through the code**, including all comments, to understand how Stan is used to represent the medical trial model.\n",
    "* **Run the code** to see inference results for the treatment group.\n",
    "* **Complete the two tasks** at the end of the notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pystan\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# For Stan we provide all known quantities as data, namely the observed data\n",
    "# and our prior hyperparameters.\n",
    "eczema_data = {\n",
    "    'treatment': {\n",
    "        'alpha': 1,  # fixed prior hyperparameters for the\n",
    "        'beta': 1,   # beta distribution\n",
    "        'num_trials': 6,  # number of trials in the data set\n",
    "        'patients': [23, 16, 16, 45, 31, 10],  # number of patients per trial\n",
    "        'improved': [20, 10, 13, 35, 22, 7]},  # number of improved patients per trial\n",
    "    'control': {\n",
    "        'alpha': 1,\n",
    "        'beta': 1,\n",
    "        'num_trials': 6,\n",
    "        'patients': [15, 18, 10, 39, 29, 10],\n",
    "        'improved': [9, 11, 4, 21, 12, 0]}}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Below is the Stan code for the medical trial data set. Note that the Stan\n",
    "# code is a string that is passed to the StanModel object below.\n",
    "\n",
    "# We have to tell Stan what data to expect, what our parameters are and what\n",
    "# the likelihood and prior are. Since the posterior is just proportional to\n",
    "# the product of the likelihood and the prior, we don't distinguish between\n",
    "# them explicitly in the model below. Every distribution we specify is\n",
    "# automatically incorporated into the product of likelihood * prior.\n",
    "\n",
    "stan_code = \"\"\"\n",
    "\n",
    "// The data block contains all known quantities - typically the observed\n",
    "// data and any constant hyperparameters.\n",
    "data {  \n",
    "    int<lower=1> num_trials;  // number of trials in the data set\n",
    "    int<lower=0> patients[num_trials];  // number of patients per trial\n",
    "    int<lower=0> improved[num_trials];  // number of improved patients per trial\n",
    "    real<lower=0> alpha;  // fixed prior hyperparameter\n",
    "    real<lower=0> beta;   // fixed prior hyperparameter\n",
    "}\n",
    "\n",
    "// The parameters block contains all unknown quantities - typically the\n",
    "// parameters of the model. Stan will generate samples from the posterior\n",
    "// distributions over all parameters.\n",
    "parameters {\n",
    "    real<lower=0,upper=1> p;  // probability of improvement - the\n",
    "                              // parameter of the binomial likelihood\n",
    "}\n",
    "\n",
    "// The model block contains all probability distributions in the model.\n",
    "// This of this as specifying the generative model for the scenario.\n",
    "model {\n",
    "    p ~ beta(alpha, beta);  // prior over p\n",
    "    for(i in 1:num_trials) {\n",
    "        improved[i] ~ binomial(patients[i], p);  // likelihood function\n",
    "    }\n",
    "}\n",
    "\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_4822bea325d0250e03828b3bc1bb8bdd NOW.\n"
     ]
    },
    {
     "ename": "CompileError",
     "evalue": "command 'gcc' failed with exit status 1",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mDistutilsExecError\u001b[0m                        Traceback (most recent call last)",
      "\u001b[0;32m~/anaconda3/lib/python3.7/distutils/unixccompiler.py\u001b[0m in \u001b[0;36m_compile\u001b[0;34m(self, obj, src, ext, cc_args, extra_postargs, pp_opts)\u001b[0m\n\u001b[1;32m    117\u001b[0m             self.spawn(compiler_so + cc_args + [src, '-o', obj] +\n\u001b[0;32m--> 118\u001b[0;31m                        extra_postargs)\n\u001b[0m\u001b[1;32m    119\u001b[0m         \u001b[0;32mexcept\u001b[0m \u001b[0mDistutilsExecError\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mmsg\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.7/distutils/ccompiler.py\u001b[0m in \u001b[0;36mspawn\u001b[0;34m(self, cmd)\u001b[0m\n\u001b[1;32m    908\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mspawn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcmd\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 909\u001b[0;31m         \u001b[0mspawn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcmd\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdry_run\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdry_run\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    910\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.7/distutils/spawn.py\u001b[0m in \u001b[0;36mspawn\u001b[0;34m(cmd, search_path, verbose, dry_run)\u001b[0m\n\u001b[1;32m     35\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mos\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mname\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'posix'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 36\u001b[0;31m         \u001b[0m_spawn_posix\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcmd\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msearch_path\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdry_run\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdry_run\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     37\u001b[0m     \u001b[0;32melif\u001b[0m \u001b[0mos\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mname\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'nt'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.7/distutils/spawn.py\u001b[0m in \u001b[0;36m_spawn_posix\u001b[0;34m(cmd, search_path, verbose, dry_run)\u001b[0m\n\u001b[1;32m    158\u001b[0m                           \u001b[0;34m\"command %r failed with exit status %d\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 159\u001b[0;31m                           % (cmd, exit_status))\n\u001b[0m\u001b[1;32m    160\u001b[0m             \u001b[0;32melif\u001b[0m \u001b[0mos\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mWIFSTOPPED\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mstatus\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mDistutilsExecError\u001b[0m: command 'gcc' failed with exit status 1",
      "\nDuring handling of the above exception, another exception occurred:\n",
      "\u001b[0;31mCompileError\u001b[0m                              Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-5-ea68f1d55dfc>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      6\u001b[0m \u001b[0;31m# to generate samples for different data sets really quickly.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      7\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 8\u001b[0;31m \u001b[0mstan_model\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpystan\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mStanModel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmodel_code\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mstan_code\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m~/anaconda3/lib/python3.7/site-packages/pystan/model.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, file, charset, model_name, model_code, stanc_ret, include_paths, boost_lib, eigen_lib, verbose, obfuscate_model_name, extra_compile_args, allow_undefined, include_dirs, includes)\u001b[0m\n\u001b[1;32m    376\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    377\u001b[0m         \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 378\u001b[0;31m             \u001b[0mbuild_extension\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrun\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    379\u001b[0m         \u001b[0;32mfinally\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    380\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mredirect_stderr\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.7/distutils/command/build_ext.py\u001b[0m in \u001b[0;36mrun\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    338\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    339\u001b[0m         \u001b[0;31m# Now actually compile and link everything.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 340\u001b[0;31m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbuild_extensions\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    341\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    342\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mcheck_extensions_list\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mextensions\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.7/distutils/command/build_ext.py\u001b[0m in \u001b[0;36mbuild_extensions\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    447\u001b[0m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_build_extensions_parallel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    448\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 449\u001b[0;31m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_build_extensions_serial\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    450\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    451\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m_build_extensions_parallel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.7/distutils/command/build_ext.py\u001b[0m in \u001b[0;36m_build_extensions_serial\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    472\u001b[0m         \u001b[0;32mfor\u001b[0m \u001b[0mext\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mextensions\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    473\u001b[0m             \u001b[0;32mwith\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_filter_build_errors\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mext\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 474\u001b[0;31m                 \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbuild_extension\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mext\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    475\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    476\u001b[0m     \u001b[0;34m@\u001b[0m\u001b[0mcontextlib\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcontextmanager\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.7/distutils/command/build_ext.py\u001b[0m in \u001b[0;36mbuild_extension\u001b[0;34m(self, ext)\u001b[0m\n\u001b[1;32m    532\u001b[0m                                          \u001b[0mdebug\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdebug\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    533\u001b[0m                                          \u001b[0mextra_postargs\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mextra_args\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 534\u001b[0;31m                                          depends=ext.depends)\n\u001b[0m\u001b[1;32m    535\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    536\u001b[0m         \u001b[0;31m# XXX outdated variable, kept here in case third-part code\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.7/distutils/ccompiler.py\u001b[0m in \u001b[0;36mcompile\u001b[0;34m(self, sources, output_dir, macros, include_dirs, debug, extra_preargs, extra_postargs, depends)\u001b[0m\n\u001b[1;32m    572\u001b[0m             \u001b[0;32mexcept\u001b[0m \u001b[0mKeyError\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    573\u001b[0m                 \u001b[0;32mcontinue\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 574\u001b[0;31m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_compile\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobj\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msrc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mext\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcc_args\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mextra_postargs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mpp_opts\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    575\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    576\u001b[0m         \u001b[0;31m# Return *all* object filenames, not just the ones we just built.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.7/distutils/unixccompiler.py\u001b[0m in \u001b[0;36m_compile\u001b[0;34m(self, obj, src, ext, cc_args, extra_postargs, pp_opts)\u001b[0m\n\u001b[1;32m    118\u001b[0m                        extra_postargs)\n\u001b[1;32m    119\u001b[0m         \u001b[0;32mexcept\u001b[0m \u001b[0mDistutilsExecError\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mmsg\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 120\u001b[0;31m             \u001b[0;32mraise\u001b[0m \u001b[0mCompileError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmsg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    121\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    122\u001b[0m     def create_static_lib(self, objects, output_libname,\n",
      "\u001b[0;31mCompileError\u001b[0m: command 'gcc' failed with exit status 1"
     ]
    }
   ],
   "source": [
    "# This cell takes a while to run. Compiling a Stan model will feel slow even\n",
    "# on simple models, but it isn't much slower for really complex models. Stan\n",
    "# is translating the model specified above to C++ code and compiling the C++\n",
    "# code to a binary that it can executed. The advantage is that the model needs\n",
    "# to be compiled only once. Once that is done, the same code can be reused\n",
    "# to generate samples for different data sets really quickly.\n",
    "\n",
    "stan_model = pystan.StanModel(model_code=stan_code)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Fit the model to the data. This will generate samples from the posterior over\n",
    "# all parameters of the model. We start by computing posteriors for the treatment\n",
    "# data.\n",
    "\n",
    "stan_results = stan_model.sampling(data=eczema_data['treatment'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Print out the mean, standard deviation and quantiles of all parameters.\n",
    "# These are approximate values derived from the samples generated by Stan.\n",
    "# You can ignore the \"lp__\" row for now. Pay attention to the row for\n",
    "# the \"p\" parameter of the model.\n",
    "#\n",
    "# The columns in the summary are\n",
    "#\n",
    "#  * mean: The expected value of the posterior over the parameter\n",
    "#  * se_mean: The estimated error in the posterior mean\n",
    "#  * sd: The standard deviation of the posterior over the parameter\n",
    "#  * 2.5%, etc.: Percentiles of the posterior over the parameter\n",
    "#  * n_eff: The effective number of samples generated by Stan. The\n",
    "#           larger this value, the better.\n",
    "#  * Rhat: An estimate of the quality of the samples. This should be\n",
    "#          close to 1.0, otherwise there might be a problem with the\n",
    "#          convergence of the sampler.\n",
    "\n",
    "print(stan_results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Specify which parameters you want to see in the summary table using\n",
    "# the \"pars\" keyword argument. Specify which percentiles you want to\n",
    "# see using the \"probs\" keyword argument.\n",
    "#\n",
    "# The statement below shows only the 2.5, 50, 97.5 percentiles for the\n",
    "# parameter p.\n",
    "\n",
    "print(stan_results.stansummary(pars=['p'], probs=[0.025, 0.5, 0.975]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Finally, we can extract the samples generated by Stan so that we\n",
    "# can plot them or calculate any other functions or expected values\n",
    "# we might be interested in.\n",
    "\n",
    "posterior_samples = stan_results.extract()\n",
    "plt.hist(posterior_samples['p'], bins=50, density=True)\n",
    "plt.title('Sampled posterior probability density for p')\n",
    "print(\n",
    "    \"Posterior 95% confidence interval for p:\",\n",
    "    np.percentile(posterior_samples['p'], [2.5, 97.5]))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 1\n",
    "* Reuse the code above to calculate the posterior 95% confidence interval for the probability of improvement in the **control group**.\n",
    "* Plot the posterior histograms of the probability of improvement in the treatment and control groups on the same figure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 2\n",
    "* Using the samples from the treatment and control group posteriors, estimate the probability that treatment is at least 19% (in absolute terms) better than control, $P(p_t > p_c + 0.19)$. We computed this result in Session 3.2 where we solved the same model analytically using the algebra of conjugate distributions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}