CS146 Session 5.2 Pre-Class Work

CS146 Session 5.2 Pre-Class Work.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "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.6.3"
    },
    "colab": {
      "name": "CS146 Session 5.2 Pre-Class Work.ipynb",
      "provenance": []
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "FfdZ_rKWcqrU",
        "colab_type": "text"
      },
      "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",
      "metadata": {
        "id": "-Ef5NgaqcqrX",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import pystan\n",
        "\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "iCvfmVX2cqra",
        "colab_type": "code",
        "colab": {}
      },
      "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]}}"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "6lVl_Dvfcqrd",
        "colab_type": "code",
        "colab": {}
      },
      "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",
        "\"\"\""
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "_lxOIvTLcqrf",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 35
        },
        "outputId": "63f913e0-e566-41d0-a1c4-59994454bfb6"
      },
      "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)"
      ],
      "execution_count": 4,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_4822bea325d0250e03828b3bc1bb8bdd NOW.\n"
          ],
          "name": "stderr"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "CI73N2P-cqri",
        "colab_type": "code",
        "colab": {}
      },
      "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'])"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "4Bm_5rV6cqrk",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 233
        },
        "outputId": "664b1197-cf7a-430e-ea56-74378d0d79ea"
      },
      "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)"
      ],
      "execution_count": 6,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Inference for Stan model: anon_model_4822bea325d0250e03828b3bc1bb8bdd.\n",
            "4 chains, each with iter=2000; warmup=1000; thin=1; \n",
            "post-warmup draws per chain=1000, total post-warmup draws=4000.\n",
            "\n",
            "       mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat\n",
            "p      0.76  8.7e-4   0.04   0.68   0.73   0.76   0.78   0.82   1615    1.0\n",
            "lp__ -80.06    0.02   0.68 -82.04 -80.22  -79.8 -79.63 -79.58   1912    1.0\n",
            "\n",
            "Samples were drawn using NUTS at Thu Oct 10 11:35:44 2019.\n",
            "For each parameter, n_eff is a crude measure of effective sample size,\n",
            "and Rhat is the potential scale reduction factor on split chains (at \n",
            "convergence, Rhat=1).\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "iP_GwtPscqrp",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 215
        },
        "outputId": "be199c35-b054-4169-9e18-e95e09554fd9"
      },
      "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]))"
      ],
      "execution_count": 7,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Inference for Stan model: anon_model_4822bea325d0250e03828b3bc1bb8bdd.\n",
            "4 chains, each with iter=2000; warmup=1000; thin=1; \n",
            "post-warmup draws per chain=1000, total post-warmup draws=4000.\n",
            "\n",
            "    mean se_mean     sd   2.5%    50%  97.5%  n_eff   Rhat\n",
            "p   0.76  8.7e-4   0.04   0.68   0.76   0.82   1615    1.0\n",
            "\n",
            "Samples were drawn using NUTS at Thu Oct 10 11:35:44 2019.\n",
            "For each parameter, n_eff is a crude measure of effective sample size,\n",
            "and Rhat is the potential scale reduction factor on split chains (at \n",
            "convergence, Rhat=1).\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "khmDn0Iacqrr",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 299
        },
        "outputId": "9803d52c-a602-4608-e8fc-c90066b65973"
      },
      "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()"
      ],
      "execution_count": 8,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Posterior 95% confidence interval for p: [0.68420083 0.82184982]\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEICAYAAABPgw/pAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAFrpJREFUeJzt3Xm0ZWV55/HvT0ZFEAhXW4aioMEJ\ntButjhnUuCK2CCq20QiKDUis1rViNI4YTUtsbcmgbXol6iLBgGJQorbSJnbEAWkV1GJSAQfEUopi\nKJkEtMOQp//Yu6jj9Y7nnHtP3be+n7Xuuvvs/Z69n/3ufZ77nvfde99UFZKkle8Bkw5AkjQeJnRJ\naoQJXZIaYUKXpEaY0CWpESZ0SWqECX0rk+SUJGct93snKckVSZ466ThmMuLxOCHJl+dY/pkkx89U\nNsmdSQ4cZrsLiOuMJG9finUPbGNsxzTJw5JckOSOJO8axzpbtf2kA9haJHkS8GfAIcB9wFXAq6vq\nGxMNbCuXpICDq+rqYddRVYeMMaQVo6qeOceyB2+eTnIGsKGq3rIccY3D4DFNcgpwUFUdN+Tq1gI/\nAXYrb5yZkwkdSLIb8GngFcA5wI7Ak4F/mWRcrUuyfVXdO8H3B0hV/euw69Cy2B+4cphkPuo5stLY\n5dJ5BEBVnV1V91XVz6vqs1X1TYAk/zbJF5LcnOQnST6cZPfNb06yPsnrk3wzyV1JTu+/Jn6m/5r4\nuSR79GVXJ6kka5NsTHJ9ktfNFliSX0vy1SS3Jbl88GtskgOSfKnfxnnAXnOs56lJNiT5o34f1id5\n8cDyhyT5YJJNSX6U5C1JHtAvO6jfzu39ez/az7+gf/vlfRfBC/v5z0pyWR/zV5M8blpdvTHJN4G7\nkmzfzzu8X75Tkvf0dbOxn95p2j68MckNwN/NsJ8nJPlKkr/q4/1OkqcNLD8/yTuSfAX4GXBgkr2T\nnJvkliRXJ3nZtNXunOSjfT1fkuTfDazv5CQ/6JddmeQ//XJIc8bye7Mcr+rrfS3wYuANfR3/7/5c\n+/i08v8zyV/Osq7D+rjv6I/dztOWz3e8Xtef27f39bBzv2yvJJ/u33dLkv87cM6sT3J4kiOAPwJe\n2Md/eZIXJLl4WgyvSfKpGWI/Azh+YP8PX+pzZEWrqm3+B9gNuBk4E3gmsMe05QcBTwd2AqaAC4D3\nDCxfD1wEPAzYB7gJuAQ4jO7D8wXgrX3Z1UABZwO7AI8FNgGH98tPAc7qp/fp4zqS7o/v0/vXU/3y\nC4F393E9Bbhj83tn2MenAvcOlP8t4C7gkf3yDwKfAnbtY/wecFK/7GzgzX0MOwNPGlhv0X2d3vz6\nsH7/nwhsR/dhXA/sNFBXlwH7AQ8cmLd5/9/W1+VD+7r+KvDfpu3Dn/b78MAZ9vOEvswfAjsALwRu\nB/bsl58P/Jiua237vswFwHv7ffv3/fH47YHjcQ/w/L7s64AfAjv0y18A7N3XzQv7On34ImL5vYGy\nX56pXoEzgLcPLHt4v53d+9fb93X+hBnqY0fgRwMxPL/fn7cv4nh9vd/HPem6Il/eL3sn8P5+vTvQ\nfavNDMf0FAbOy/7Y3QI8emDepcDvzHLuTt//JT1HVvLPxAPYWn6AR/cnzob+YJ8LPGyWss8FLh14\nvR548cDrjwPvG3j9SuCT/fTq/sP6qIHlfwac3k/ff/IDbwQ+NG3b/9x/6Fb1ce4ysOzvmT+hD5Y/\nB/jj/oN8N/CYgWX/BTi/n/4gcBqw7wzrnZ7Q37f5wzUw77vAbw3U1UunLR/88P8AOHJg2TOA9QP7\ncDew8xzH8QRgI31i6ed9HXhJP30+8LaBZfvRjZnsOjDvncAZA8fjooFlDwCuB548y/YvA45eRCyL\nTuj9vM8AL+unn0XXJTFTPE+ZIYavsiWhL+R4HTftXH1/P/02ukbAQTNsd/CYnsK087Lf7jv66UOA\nW+n/iMywrl/Y/6U+R1byj10uvaq6qqpOqKp9gUPpWiTvgftH2T+S5LokPwXO4pe7N24cmP75DK8f\n/IvFuXZg+kf99qbbH3hB/5X2tiS3AU+ia6HtDdxaVXdNW89cZiq/d78vO0x7/4/oviEAvAEI8PV0\nVy+8dI5t7A+8dlrM+03bv2tnfiv05abHMfjeTVX1/+Z4P8B11X9KZ1nH4Pb3Bm6pqjumld9npvLV\n9bdv2Ly+JP95oLviNrpzZ/DcmC+WYZ0JbB5kPA740Czl9p4lhs0WcrxuGJj+GVvO5T8HrgY+m+Sa\nJCcvMv4XJQnwEuCcqlromNVynCMrkgl9BlX1HbpWwaH9rP9O12J6bFXtRvcByoib2W9gehVdi2G6\na+la6LsP/OxSVafStRL3SLLLtPXMZabyG+muILiH7sM9uOw6gKq6oapeVlV707Xc35vkoFm2cS1d\ny2sw5gdV1dkDZeYa3No4QxyDdbOQgbF9+kSxkHVsBPZMsuu08tcNvL7/WPV9xPsCG5PsD/wN8PvA\nr1TV7sC3+cVzY75YFmKmff4k8Lgkh9K10D88y3uvnyWGzRZyvGYOquqOqnptVR0IPAd4zSx90b8U\nf1VdRNeSfjLwImb/gzST5ThHViQTOpDkUUlem2Tf/vV+wLF0/XTQ9SvfCdyeZB/g9WPY7B8neVCS\nQ4ATgY/OUOYs4NlJnpFkuyQ794M++1bVj4B1wJ8k2THdZZfPXsB2N5d/Ml0i+Iequo+u++UdSXbt\nE9Vr+u3TD2Lt27//VroPzOYrQ24EBq+X/hvg5UmemM4uSY6aljDncjbwliRTSfYC/uvmOBbhocAf\nJNkhyQvoutP+aaaCVXUtXRfEO/v6fRxw0rRtPiHJ85JsD7ya7uqni+jGQIquz50kJ7KlEbDoWOYw\nvY7pW6Afo+tm+3pV/XiW915I19W2OYbnAb86sHzo45VuMPWgPjHeTtd1NdMVQzcCqzcPmA74IPBX\nwD1VNev1+jNY1nNkJTGhd+6gGxT6WpK76D6s3wZe2y//E+DxdCftPwKfGMM2v0T3dfXzwF9U1Wen\nF+iTzdF0VwlsomtNvZ4tx+1Ffdy3AG+l+4DM5Qa6hLyRrkX38v7bCHT9/HcB1wBfpksUH+iX/Qe6\nurmTbmzhVVV1Tb/sFODM/uv671bVOuBldB/UW/t9PGH+6rjf2+n+UH0T+Bbd4PJib4L5GnAw3TeP\ndwDPr6qb5yh/LN3Yxkbgf9ENYH9uYPmn6AbObqXrHnheVd1TVVcC76JLmjfSDXB/ZcRYZnI68Ji+\njj85MP/Mfpuztm6r6m7geXTH4JZ+Pz4xsHyU43Uw8Dm6xs6FwHur6oszlPuH/vfNSS4ZmP8huj+A\ni03GkzhHVoTNI9JaJklWs+UqiWW7Pjbd5Y5n9WMEzUpyAt1A45MmHctSS7IK+A7wb6rqp5OOZ7GS\nPJDuCpvHV9X3l3G7J9DoOWILXVqB+u6L1wAfWYnJvPcK4BvLmcxb552i0grTD2zfSHdlxhETDmco\nSdbTDR4/d8KhNGXeLpckH6AbPLupqg7t5/053QDc3XTXhJ5YVbctcaySpDkspMvlDH65FXAecGhV\nPY7ujsI3jTkuSdIizdvlUlUX9AN5g/MGr8i4iO524nnttddetXr16nnLSZK2uPjii39SVVPzlRtH\nH/pLmfkaagDSPVxoLcCqVatYt27dGDYpSduOJPPdBQ6MeJVLkjfT3bQw211qVNVpVbWmqtZMTc37\nB0aSNKShW+j9tZzPAp5WXswuSRM3VEJP94zjN9A9ke1n4w1JkjSMebtckpxNd1vvI9M9OP4kutuE\ndwXO65809/4ljlOSNI+FXOVy7AyzT1+CWCRJI/DWf0lqhAldkhphQpekRpjQJakRPm1R25TVJ//j\njPPXn3rUMkcijZ8tdElqhAldkhphQpekRtiHLs1htj53sN9dWx8Tupo0VyKWWmWXiyQ1wha6NCQv\ngdTWxha6JDXChC5JjbDLRSvauAY/HURVC2yhS1IjTOiS1AgTuiQ1woQuSY0woUtSI0zoktQIE7ok\nNcKELkmNMKFLUiO8U1QrgndySvOzhS5JjbCFrq2KLXFpePO20JN8IMlNSb49MG/PJOcl+X7/e4+l\nDVOSNJ+FdLmcARwxbd7JwOer6mDg8/1rSdIEzZvQq+oC4JZps48GzuynzwSeO+a4JEmLNGwf+sOq\n6vp++gbgYbMVTLIWWAuwatWqITcnrXz+yzottZGvcqmqAmqO5adV1ZqqWjM1NTXq5iRJsxg2od+Y\n5OEA/e+bxheSJGkYwyb0c4Hj++njgU+NJxxJ0rAWctni2cCFwCOTbEhyEnAq8PQk3wcO719LkiZo\n3kHRqjp2lkVPG3MskqQReOu/JDXChC5JjTChS1IjTOiS1AgTuiQ1woQuSY0woUtSI0zoktQI/2OR\nNGb+1yVNii10SWqECV2SGmGXiybCbglp/GyhS1IjTOiS1AgTuiQ1woQuSY0woUtSI0zoktQIE7ok\nNcKELkmNMKFLUiNM6JLUCBO6JDXChC5JjTChS1IjTOiS1AgTuiQ1YqSEnuQPk1yR5NtJzk6y87gC\nkyQtztD/4CLJPsAfAI+pqp8nOQc4BjhjTLGpAf4jC2n5jNrlsj3wwCTbAw8CNo4ekiRpGEMn9Kq6\nDvgL4MfA9cDtVfXZ6eWSrE2yLsm6TZs2DR+pJGlOQyf0JHsARwMHAHsDuyQ5bnq5qjqtqtZU1Zqp\nqanhI5UkzWmULpfDgR9W1aaqugf4BPAb4wlLkrRYoyT0HwO/luRBSQI8DbhqPGFJkhZr6Ktcqupr\nST4GXALcC1wKnDauwLR1mu2qlfWnHrXMkUiabuiEDlBVbwXeOqZYJEkj8E5RSWqECV2SGmFCl6RG\nmNAlqREmdElqxEhXuUib+RAuafJM6NIK470Amo1dLpLUCBO6JDXChC5JjTChS1IjTOiS1AgTuiQ1\nwssWpQnzMkSNiy10SWqECV2SGmFCl6RGmNAlqREmdElqhAldkhphQpekRpjQJakRJnRJaoQJXZIa\nYUKXpEb4LBdpK+X/adVi2UKXpEaMlNCT7J7kY0m+k+SqJL8+rsAkSYszapfLXwL/p6qen2RH4EFj\niEkT5ld9aWUaOqEneQjwFOAEgKq6G7h7PGFJkhZrlC6XA4BNwN8luTTJ3ybZZXqhJGuTrEuybtOm\nTSNsTpI0l1ES+vbA44H3VdVhwF3AydMLVdVpVbWmqtZMTU2NsDlJ0lxGSegbgA1V9bX+9cfoErwk\naQKGTuhVdQNwbZJH9rOeBlw5lqgkSYs26lUurwQ+3F/hcg1w4ughSZKGMVJCr6rLgDVjikXSEpjt\nMtT1px61zJFoqXmnqCQ1woQuSY0woUtSI0zoktQIH58rNcJn8MgWuiQ1woQuSY0woUtSI0zoktQI\nE7okNcKELkmNMKFLUiNM6JLUCBO6JDXChC5JjTChS1IjTOiS1AgfziVto/xPRu2xhS5JjbCFvg3z\ncatSW2yhS1IjTOiS1AgTuiQ1woQuSY0woUtSI0zoktQIE7okNWLkhJ5kuySXJvn0OAKSJA1nHC30\nVwFXjWE9kqQRjJTQk+wLHAX87XjCkSQNa9QW+nuANwD/OluBJGuTrEuybtOmTSNuTpI0m6ETepJn\nATdV1cVzlauq06pqTVWtmZqaGnZzkqR5jPJwrt8EnpPkSGBnYLckZ1XVceMJTYvl41ClbdvQLfSq\nelNV7VtVq4FjgC+YzCVpcrwOXZIaMZbnoVfV+cD541iXxs/nnkvbBlvoktQIE7okNcKELkmNMKFL\nUiNM6JLUiLFc5SKpHd6gtnLZQpekRpjQJakRJnRJaoQJXZIaYUKXpEaY0CWpESZ0SWqECV2SGuGN\nRZIWxBuOtn620CWpESZ0SWqECV2SGmFCl6RGmNAlqREmdElqhAldkhphQpekRpjQJakRJnRJaoQJ\nXZIa4bNcVqDZnqkhads2dAs9yX5JvpjkyiRXJHnVOAOTJC3OKC30e4HXVtUlSXYFLk5yXlVdOabY\nJEmLMHRCr6rrgev76TuSXAXsA5jQx8SuFUmLMZY+9CSrgcOAr41jfZJWjrkaHj4rfXmNfJVLkgcD\nHwdeXVU/nWH52iTrkqzbtGnTqJuTJM1ipISeZAe6ZP7hqvrETGWq6rSqWlNVa6ampkbZnCRpDqNc\n5RLgdOCqqnr3+EKSJA1jlBb6bwIvAX47yWX9z5FjikuStEijXOXyZSBjjEXSNsJ/OL00vPVfkhph\nQpekRpjQJakRPpxrGXnnp7Y1nvPLyxa6JDXChC5JjTChS1IjTOiS1AgTuiQ1woQuSY0woUtSI0zo\nktQIbywagQ8YkrQ1sYUuSY0woUtSI0zoktQIE7okNcJB0SXgE+YkTYItdElqhAldkhrRbJeL14hL\n2tbYQpekRjTbQp/NXAOWs7XeHeSUtBJscwl9LiZuabLG9RncVrtWTeiSthmtj63Zhy5JjTChS1Ij\nVkyXS+tflSSNz7Y6HjZSCz3JEUm+m+TqJCePKyhJ0uIN3UJPsh3w18DTgQ3AN5KcW1VXjiu4hdhW\n/xJLWnorrWdglC6XXwWurqprAJJ8BDgaWNaELknLbZiG5HL8ERgloe8DXDvwegPwxOmFkqwF1vYv\n70zy3RG2OZe9gJ8s0bpXGutiC+uiYz1s8Ut1kT9d+o2OuI39F1JoyQdFq+o04LSl3k6SdVW1Zqm3\nsxJYF1tYFx3rYYuW62KUQdHrgP0GXu/bz5MkTcAoCf0bwMFJDkiyI3AMcO54wpIkLdbQXS5VdW+S\n3wf+GdgO+EBVXTG2yBZvybt1VhDrYgvromM9bNFsXaSqJh2DJGkMvPVfkhphQpekRqyIhL6QRwwk\n+d0kVya5IsnfD8y/L8ll/c+KHrSdrx6S/I+Bff1ektsGlh2f5Pv9z/HLG/n4jVgXzZwTsKC6WJXk\ni0kuTfLNJEcOLHtT/77vJnnG8kY+XsPWQ5LVSX4+cE68f/mjH5Oq2qp/6AZcfwAcCOwIXA48ZlqZ\ng4FLgT361w8dWHbnpPdhuephWvlX0g1UA+wJXNP/3qOf3mPS+zSJumjpnFhoXdANAr6in34MsH5g\n+nJgJ+CAfj3bTXqfJlAPq4FvT3ofxvGzElro9z9ioKruBjY/YmDQy4C/rqpbAarqpmWOcTkspB4G\nHQuc3U8/Azivqm7p6+g84IgljXZpjVIXrVlIXRSwWz/9EGBjP3008JGq+peq+iFwdb++lWiUemjG\nSkjoMz1iYJ9pZR4BPCLJV5JclGQwWe2cZF0//7lLHewSWkg9AJBkf7oW1xcW+94VYpS6gHbOCVhY\nXZwCHJdkA/BPdN9YFvrelWKUegA4oO+K+VKSJy9ppEtoxTwPfR7b03W7PJXujtULkjy2qm4D9q+q\n65IcCHwhybeq6gcTjHU5HAN8rKrum3QgW4GZ6mJbOyeOBc6oqncl+XXgQ0kOnXRQEzBbPVwPrKqq\nm5M8AfhkkkOq6qcTjXYIK6GFvpBHDGwAzq2qe/qvjt+jS/BU1XX972uA84HDljrgJbKYRy0cwy92\nMbT2mIZR6qKlcwIWVhcnAecAVNWFwM50D6hq6bwYuh76Lqeb+/kX0/XFP2LJI14Kk+7En++HrvV9\nDd3X5s2DHYdMK3MEcGY/vRfdV69foRsA3Glg/veZY/Bsa/5ZSD305R4FrKe/aayftyfww74+9uin\n95z0Pk2oLpo5JxZaF8BngBP66UfT9R0HOIRfHBS9hpU7KDpKPUxt3m+6QdXrVurnY+IBLPBgHUnX\n6v4B8OZ+3tuA5/TTAd5N9yz2bwHH9PN/o399ef/7pEnvy1LWQ//6FODUGd77UrpBr6uBEye9L5Oq\ni9bOiYXUBd0VHV/p9/ky4D8OvPfN/fu+Czxz0vsyiXoAfge4op93CfDsSe/LsD/e+i9JjVgJfeiS\npAUwoUtSI0zoktQIE7okNcKELkmNMKFLUiNM6JLUiP8P3F80jgX1vmwAAAAASUVORK5CYII=\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "tWPB5p1lcqrt",
        "colab_type": "text"
      },
      "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",
      "metadata": {
        "id": "0_nG0el6cqru",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 479
        },
        "outputId": "79511c08-e6b5-4c1a-e210-c5816fb83fd1"
      },
      "source": [
        "stan_results_control = stan_model.sampling(data=eczema_data['control']) # Change 'treatment' to 'control' for control group data\n",
        "print(stan_results_control.stansummary(pars=['p'], probs=[0.025, 0.5, 0.975]))\n",
        "\n",
        "posterior_samples_treatment = stan_results.extract()\n",
        "posterior_samples_control = stan_results_control\n",
        "plt.hist(posterior_samples_treatment['p'], bins=50, density=True)\n",
        "plt.hist(posterior_samples_control['p'], bins=50, density=True)\n",
        "plt.title('Sampled posterior probability density for p for treatment and control')\n",
        "\n",
        "plt.show()"
      ],
      "execution_count": 11,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Inference for Stan model: anon_model_4822bea325d0250e03828b3bc1bb8bdd.\n",
            "4 chains, each with iter=2000; warmup=1000; thin=1; \n",
            "post-warmup draws per chain=1000, total post-warmup draws=4000.\n",
            "\n",
            "    mean se_mean     sd   2.5%    50%  97.5%  n_eff   Rhat\n",
            "p   0.47  1.2e-3   0.04   0.38   0.47   0.56   1413    1.0\n",
            "\n",
            "Samples were drawn using NUTS at Thu Oct 10 12:20:47 2019.\n",
            "For each parameter, n_eff is a crude measure of effective sample size,\n",
            "and Rhat is the potential scale reduction factor on split chains (at \n",
            "convergence, Rhat=1).\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAa4AAAEICAYAAAAHsBBpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAHD9JREFUeJzt3Xm0JVV96PHvj3mQUVoiNNCOKBLj\n0EafUWEJviAQMDiBwdCKIq4oGtGIT32iQsREjeQ5pR0CogERJ+IQJSoxyGQjaAQcsaWZWyYBB0B/\n74+9L10czr333KFPVd3+ftbq1eecqnPqt2vvU7/au/apG5mJJEl9sV7bAUiSNBMmLklSr5i4JEm9\nYuKSJPWKiUuS1CsmLklSr/QmcUXEsRHxiXG/t00RcWlE7Nl2HMPMsT6WRcQ5Uyz/SkQcNmzdiLg9\nIh48m+2OENdJEXHc2vjsxjbmrU4jYvuI+FZE3BYR756Pz5zh9neNiEvq9o8a9/Y1tem+Z10VEWdH\nxEumWmfaxBURT4mIcyPi1oi4KSK+HRFPmL8wF6aIyIh46Fw+IzMflZlnz1NIvZGZz8zMkydZdr/M\nvALGk2jmW7NO5+GE6gjgl8CWmXn0fMQ3Q38HfDMzt8jMf57rh43jBDMiVkbE3mtzG1Nsu5eJZCZG\nSTrzYcrEFRFbAl8E/h+wLbAj8Fbgd2s7sHVZRGzQ8vsjInrTG1+H7QJclrO4i8Bc20hj+5fO5o2z\n2b7tsv/mqd1BZk76D1gK3DLF8ocA3wBupJz5fRLYurF8JfA64PvAHcBHge2BrwC3Af8JbFPXXQIk\n5SzyGuBa4LWNzzoW+ETj+ZOAc4FbgO8BezaWPQj4r7qNs4D3Nd87UIY9gauA/1PLsBL4q8byrYCP\nA6uBXwBvAtaryx5at3Nrfe+n6uvfqmW5A7gdeH59fX/gkhrzucCjB/bV6+u++h2wQX1t77p8Y+C9\ndd9cUx9vPFCG1wPXAacMKecy4Nt1X9wK/BDYq7H8bOD4us5vatl2AM4EbgJ+Crx0oD7OAD5V9/N3\ngT9pLD8G+FlddhnwlzOM5SWNdc9pLMsa2xHAXcCddR//O6WtfWag3P8MnDhJ3T+2xn1bLcdpwHGN\n5dPV12trfd1a379JXbYd5YTvlrrv/ps1bWYlsDewT439rhr/94DnAhcNxPga4AtDYj9poPx7r+02\nMrDuN4DfA7+t2384U39XJj77nyjHi+MGPu8++2OKdrkV5VhyLXA1cByw/nTHJOAU4A/1c26n9BiX\nUNrUi4BVwM3AkcATat3eArxvINYXA5fXdb8K7DLQPo8EflLf+34ggEfWffX7uu2hx9Uax+WUNnkF\n8LIhx6qjgRtq+V/UWH5/yvf1V8CFwNtpfHeGbOsprDmGrgKWjXDMWwacA7yrlv/nwDPrsuMH2sT7\nGvvkb+o++Xl97cnAdyjt7DvAk4d9/yeNfZrEtWVtACcDz6QmmcbyhwLPoHxhFlEO2O8d+HKfT0lW\nO9ad/V3KAWOT2sDeMpC4TgU2B/647riJA/ex1ORTP+tGYF9Kr/EZ9fmiuvw84D01rqfVRjBV4rq7\nsf4elISza13+ceALwBY1xh8Dh9dlpwJvrDFsAjxl8AA7cJC8AXgisD5wWN0/Gzf21SXATsCmzYNc\nffy2ui8fUPf1ucDbB8rwzlqGTSc5KN0N/C2wIfD82mi2bTSWK4FHUZLmhrU+P1DL9phaH09v1Mdd\nwHPquq+lNOIN6/LnUhLfenVbdwAPnEEsUyauxsG7mWgeWLczcaDaoO7zxw/ZHxtRvpQTMTynlue4\nGdTXhbWM21IONkfWZe8APlQ/d0PgqUAMqdNjuffJ2MaURPfIxmsXA8+epO0Oln+ttpEh699TTyN8\nVyY++5W1XoZt/177Y4p2+TngXyjHiQfUenjZDI5JezeeL6G0qQ9R2vn/phx4P18/e+K4tUdd/0DK\nSdwjazxvAs4daJ9fBLYGdqZ8Z/YZ1pYn2af7UZJvUI5FvwYeN1CHb6v7Yd+6fOLk/zTg9Lpfdqck\n9aHbo/SWbwMOqZ91f+AxI9bjXcBLKd+Ll1NOkqJRXy8Z2FZSOhDbApvW/28GXlj34SH1+f3nJXHV\nD3kk5QtyVd1pZwLbT7Lus4CLBxpJs/fyGeCDjeevBD4/0IAe0Vj+D8BHBxs15azxlIFtf5VycNm5\nxrl5Y9m/MX3iaq5/OvDmWjF3Ars1lr0MOLtRwcuBxUM+dzBxfZB6EGm89iPWfCFWAi8eWL6SNQe5\nnwH7Npb9ObCyUYY7qWf8k5RzWbOB1dcuBF7YaCxvayzbiXL2tEXjtXcAJzXq4/zGsvUoZ4BPnWT7\nlwAHziCWGSeu+tpXqD1DSo/pskniedqQGM5lTeIapb4OHWirH6qP30b54j90yHabdXos9z1QfxA4\nvj5+FOULvfEkZbhX+dd2GxmyfrOepvuuLAOunGzbU+yPs7l3u9yeMiKxaeO1QyjX2kY9Jg1LXDs2\nXruROkpSn38GeHWjfR0+0O5/Te111c9qnsCeDhwzrC2P8o+SQF/VqMPfABs0lt9AGX1an5JQmsfP\nv59se8AbgM8NeX2UevxpY9lmtcx/NNgmBr6zT288fyFw4cA657Gmx3efzxj8N+14cWZenpnLMnMx\nJYvvQBmCmJjVdFpEXB0RvwI+QRkmabq+8fg3Q57fb2D9VY3Hv6jbG7QL8NyIuGXiH6Xb+8C6/s2Z\necfA50xl2Po71LJsOPD+X1DOwqAMNQRwYZ0t9uIptrELcPRAzDsNlG/V8LdCXW8wjuZ7V2fmb6d4\nP8DVWVvGJJ/R3P4OwE2ZedvA+jsOWz8z/0A5udkBICL+us44myjr7ty7bUwXy2ydDBxaHx9KGR4a\nZodJYpgwSn1d13j8a9a05X+knJV/LSKuiIhjZhj/CyIiKF/w0zNz1GvK42gjk5nuuwJTt++pNN+3\nS93OtY16+RdK72jUY9Iwox6ndgFObGz7JsoxoFnOydrFtCLimRFxfp0IdwulV9WM/8bMvHvI5y+i\n9F4Gj5+T2YlyojNolHq8p3yZ+ev6cLoyDh5bBmMb3MaUZnShMzN/SDnL272+9PeUbPrHmbkl5UAR\nM/nMIXZqPN6ZcgY4aBWlx7V149/mmXkC5ax/m4jYfOBzpjJs/WsoY+R3URprc9nVAJl5XWa+NDN3\noJyVfGCKmYSrKGfSzZg3y8xTG+vkJO+lxjMYR3PfTPXeCTvWA+Ion3ENsG1EbDGw/tWN5/fUVb1o\nvhi4JiJ2AT4MvILS/d8a+AH3bhvTxTKKYWX+PPDoiNid0uP65CTvvXaSGCaMUl/Dg8q8LTOPzswH\nAwcAr4mIvUaJPzPPp5zxPhV4AZMn3mHG0UYmM+V3ZcTtT7a8+foqSo9ru0a9bJmZj6rLpzsmjbIP\nprKKMizZbBebZua5I7x3ym1HxMaU3t27KKNaWwNfZrRj6mrKyNHg8XMyqyhDkoNGqcepjFKHg+10\nptuYdlbhIyLi6IhYXJ/vROmWn19X2YJyEe7WiNiRcnF8rt4cEZtFxKMoFyo/NWSdTwB/ERF/HhHr\nR8QmEbFnRCzOzF8AK4C3RsRGEfEU4C9G2O7E+k+lHPA+nZm/p3T1j4+ILeoB+TV1+0TEcyf2DWVI\nJykXf6GcsTV/b/Rh4MiIeGKdHbV5ROw3kBimcirwpohYFBHbAf93Io4ZeABwVERsGBHPpQwDf3nY\nipm5ijJ09o66fx8NHD6wzcdHxEF1ptCrKQeU8ylj7En5MhERL2LNyc6MY5nC4D6m9ijOoAwPX5iZ\nV07y3vMoX/SJGA4C/rSxfNb1FRH7R8RDawK4lTLk+ochq14PLBkyU+7jlAkSd2XmTKZPj7WNNE33\nXRnRZPujuZ1rga8B746ILSNivYh4SETsUVeZ7ph0nzYzQx8C3lCPT0TEVnU/jeJ6YHFEbDTJ8o0o\n1+ZWA3dHxDMp19ymVff/Z4Fj6/FzN8qlk8l8Etg7Ip4XERtExP0j4jHzUI+j7N8vAw+PiBfUbT8f\n2I1ybXAk0/W4bqNcnL4gIu6gHJR+QJnVAmVq/OMoX84vUXbcXP0XZZjl68C7MvNrgyvUg+qBlJmA\nqylnD69jTXleUOO+CXgL5UAwlesoiecaSoUeWXuXUK7D3UGZ4XMO5YD4sbrsCZR9czvl2t+rsv7G\niDJef3IdUnheZq6gXNB8X93WTynjxaM6jpKQvw/8D2WSy0x/w3QB8DDKWdXxwHMy88Yp1j+Ecg3g\nGsoF8bdk5n82ln+BcgF/4kLrQZl5V2ZeBrybkhyup0y0+fYcYxnmo8BudR9/vvH6yXWbk/ZWMvNO\n4CBKHdxUy/HZxvK51NfDKDNmb6fsgw9k5jeHrPfp+v+NEfHdxuunUBL9TJNOG22kaarvyigm2x+D\n/ppykL+MUjdnUC4TwPTHpHdQkvstEfHaGcQGQGZ+jjLB5bQ6FPkDysS1UXyD8vOB6yLil0M++zbg\nKEriuJlyHDtzBuG9gjJkdx1lZOxfJ1uxntDtSzmW30S5Bv0ndfFc6vFE4DkRcXNEDP1tX21P+9dt\n30i55LJ/Zt5nn0xmYiZI6yJiCWtmpd099drzut09KReEF0+3bp9FxDLKBc+ntB3L2hYRO1Omcv9R\nZv6q7XhmKiI2pVx0f1xm/mSM213GOtJG1G/+mE8LSh1meg1wWh+TVvVy4DvjTFpSn8zPr5ilDqgT\nbK6nzFDap+VwZiUiVlIuxj+r5VCkzmplqDAiPkYZ47whM3evr/0jZRLFnZRpmi/KzFvGHpwkqdPa\nGio8ifueEZ8F7J6Zj6b8UvsN4w5KktR9rQwVZua36mSM5mvN2YPnU27BM63tttsulyxZMu16kqQ1\nLrrool9m5qK245iNrl7jejHDf78FQEQcQbnJKjvvvDMrVqwYV1yStCBExHR3FOqszs0qjIg3Un4Y\nOtkdD8jM5Zm5NDOXLlrUyxMGSdIsdarHVX9Hsj/lTyl04wdmkqRO6Uziioh9KL+g3qNx40ZJku6l\nlaHCiDiVciucXSPiqog4nHJrnS2As6LcVfxDbcQmSeq2tmYVHjLk5Y+OPRBJUu90bnKGJElTMXFJ\nknrFxCVJ6hUTlySpVzozHV6SZmLJMV+65/HKE/ZrMRKNmz0uSVKvmLgkSb1i4pIk9YqJS9KCseSY\nL93r2pcWJidnSOoVE5PscUmSesUel6QFx6nyC5s9LklSr5i4JEm94lChpM6bbkKGEzbWLfa4JEm9\nYuKSJPWKiUuS1CsmLklSr5i4JEm9YuKSJPWKiUuS1CsmLklSr5i4JEm94p0zJHWWd8TQMPa4JEm9\nYo9LUqfYy9J0WutxRcTHIuKGiPhB47VtI+KsiPhJ/X+btuKTJHVTm0OFJwH7DLx2DPD1zHwY8PX6\nXJKke7SWuDLzW8BNAy8fCJxcH58MPGusQUmSOq9r17i2z8xr6+PrgO2HrRQRRwBHAOy8885jCk1S\nnzWvna08Yb8WI9FcdXZWYWYmkJMsW56ZSzNz6aJFi8YcmSSpTV1LXNdHxAMB6v83tByPJKljupa4\nzgQOq48PA77QYiySpA5qczr8qcB5wK4RcVVEHA6cADwjIn4C7F2fS5J0j9YmZ2TmIZMs2musgUiS\neqVrQ4WSJE3JxCVJ6hUTlySpV0xckqReMXFJknrFxCVJ6hUTlySpV7p2k131ybFbDTy/tZ04JK1T\nTFySFjT/ovLC41ChJKlXTFySpF5xqFBSJzikp1HZ45Ik9YqJS5LUKyYuSVKvmLgkSb1i4pIk9Yqz\nCjUezbtseIcNSXNgj0uS1CsmLklSrzhUqLVj8Aa8kjRPTFyaPyYrSWNg4tL4+edQJM2BiUvtc8ah\npBlwcoYkqVdMXJKkXjFxSZJ6pXOJKyL+NiIujYgfRMSpEbFJ2zFJkrqjU5MzImJH4Chgt8z8TUSc\nDhwMnNRqYJLWCv94pGajcz0uSjLdNCI2ADYDrmk5HklSh3QqcWXm1cC7gCuBa4FbM/Nrg+tFxBER\nsSIiVqxevXrcYUqSWtSpxBUR2wAHAg8CdgA2j4hDB9fLzOWZuTQzly5atGjcYUqSWtSpxAXsDfw8\nM1dn5l3AZ4EntxyTJKlDupa4rgSeFBGbRUQAewGXtxyTJKlDOjWrMDMviIgzgO8CdwMXA8vbjUr3\n4o10NUvNGYQrT9ivxUjUd51KXACZ+RbgLW3HsU7z3oGSOqxrQ4WSJE2pcz0udYxDg5I6xh6XJKlX\nTFySpF4xcUmSesVrXJLGzpvrai7scUkSJZmaUPvBxCVJ6hUTlySpV0xckqReMXFJknrFWYXqlsE7\ndXivREkD7HFJknrFxCVpnePU934zcUmSesXEJUnqFROXJKlXnFUo/+aWpF6xxyVJ6hUTlySpV0xc\nkqReMXFJknrFxCVJ6hUTlySpV0xckqRe8Xdc6rbmb8y8U7zmmfcr7Cd7XJKkXulc4oqIrSPijIj4\nYURcHhH/q+2YJEnd0cWhwhOB/8jM50TERsBmbQckafYcjtN861TiioitgKcBywAy807gzjZjkiR1\nS9eGCh8ErAb+NSIujoiPRMTmgytFxBERsSIiVqxevXr8UUqSWtOpHhclnscBr8zMCyLiROAY4M3N\nlTJzObAcYOnSpTn2KPvKu8BLWgC61uO6CrgqMy+oz8+gJDJJkoCOJa7MvA5YFRG71pf2Ai5rMSRJ\nUsd0bagQ4JXAJ+uMwiuAF7UcjySpQzqXuDLzEmBp23Gogwav0XknDa1lE1P5V56wX8uRqKlTQ4WS\nJE3HxCVJ6hUTlySpV0xckqRe6dzkDElqk/dW7D57XJKkXjFxSZJ6xcQlSeoVE5ckqVdMXJKkXjFx\nSZJ6xcQlSeoVE5ckqVdMXJKkXjFxSZJ6xcQlSeoVE5ckqVe8ya4kTaN5413/GnL77HFJknrFHpek\neeefBtHaZI9LktQr9rjUX8du1Xh8a3txSBore1ySpF6xx7XQ2AuRtMDZ45Ik9YqJS5LUKyYuSVKv\ndDJxRcT6EXFxRHyx7VgkSd3SycQFvAq4vO0gJEnd07nEFRGLgf2Aj7QdiySpezqXuID3An8H/GGy\nFSLiiIhYERErVq9ePb7IJEmt69TvuCJif+CGzLwoIvacbL3MXA4sB1i6dGmOKTx1WfP3a+Bv2KQF\nrFOJC/gz4ICI2BfYBNgyIj6RmYe2HJekSfgnPzRunRoqzMw3ZObizFwCHAx8w6QlSWrqVOKSJGk6\nXRsqvEdmng2c3XIYkmbAv8Olcehs4tI8GJywIEkLgEOFkqReMXFJknrFxCVJ6hUTlySpV0xckjQD\nS475krMnW2bikiT1itPh+84p75LWMfa4JEm9YuKSJPWKiUuS1CsmLklSr5i4JEm94qxCLUzN2Zb+\nNWRpQTFxSdIs+Jef2+NQoSSpV0xckqReMXFJknrFxCVJ6hUTlySpV0xckqRecTq8Fr7BO+j7uy6p\n1+xxSZJ6xcQlSeoVE5ckqVdMXJKkXnFyhqRZad6rTxqnTiWuiNgJ+DiwPZDA8sw8sd2oOmhwlpwk\nrUM6lbiAu4GjM/O7EbEFcFFEnJWZl7UdmCSpGzqVuDLzWuDa+vi2iLgc2BFYtxOXPSx1hMOD6oLO\nTs6IiCXAY4EL2o1Ekqa25JgvmdTHqJOJKyLuB3wGeHVm/mrI8iMiYkVErFi9evX4A5QktaZziSsi\nNqQkrU9m5meHrZOZyzNzaWYuXbRo0XgDlCS1qlPXuCIigI8Cl2fme9qORwtU85qh9y2UeqdrPa4/\nA14IPD0iLqn/9m07KElSd3Sqx5WZ5wDRdhySNFcTkzVWnrBfy5EsPF3rcUmSNCUTlySpV0xckqRe\n6dQ1Lql1zji8hz+onTn32XiYuLrK2zxJ0lAmLq3bPEGQesdrXJKkXjFxSZJ6xcQlSeoVE5ckqVec\nnNEVThKQpJHY45Ik9Yo9rjbZy+q2wfpZwD9Ibv5w1pvCquvscUmSesXEJUnqFROXJKlXTFySpF5x\ncoake/EO5+o6e1ySpF4xcUmSesWhwnHyd1saM3+fpYXIxCWNyr+OLHWCiWtts5eljpjofTV7Xk7E\nUB+ZuKR1jMlqvKbb3w7hzpyJS5oPDiNqjrweOTpnFUqSesUe13zwOta6Z6o6X4fuKi+1wcQlrW1j\nGkZ0qKmfvOY4c51LXBGxD3AisD7wkcw8oeWQ7sseliS1plOJKyLWB94PPAO4CvhORJyZmZetlQ3O\nJAE53KP5MOUQ40Abm0NPzbP4hcFe9HCdSlzAnwI/zcwrACLiNOBAYO0krpmwl6W1zetmmsKwk5F1\nNZl1LXHtCKxqPL8KeOLgShFxBHBEfXp7RPxois/cDvjlvEXYHQu1XLBwyzZ/5XprzMvHzCPrbB7F\nO+d3vUnsOqd3t6hriWskmbkcWD7KuhGxIjOXruWQxm6hlgsWbtkWarlg4ZZtoZYLStnajmG2uvY7\nrquBnRrPF9fXJEkCupe4vgM8LCIeFBEbAQcDZ7YckySpQzo1VJiZd0fEK4CvUqbDfywzL53jx440\npNhDC7VcsHDLtlDLBQu3bAu1XNDjskVmth2DJEkj69pQoSRJUzJxSZJ6ZcEkrojYJyJ+FBE/jYhj\nhiw/MiL+JyIuiYhzImK3NuKcqenK1Vjv2RGREdGbqbsj1NmyiFhd6+ySiHhJG3HO1Ch1FhHPi4jL\nIuLSiPi3ccc4WyPU2T816uvHEXFLG3HO1Ajl2jkivhkRF0fE9yNi3zbinKkRyrVLRHy9lunsiFjc\nRpwzlpm9/0eZyPEz4MHARsD3gN0G1tmy8fgA4D/ajns+ylXX2wL4FnA+sLTtuOexzpYB72s71rVQ\nrocBFwPb1OcPaDvu+SrbwPqvpEywaj32eaiz5cDL6+PdgJVtxz1P5fo0cFh9/HTglLbjHuXfQulx\n3XOrqMy8E5i4VdQ9MvNXjaebA32YlTJtuaq3A+8EfjvO4OZo1LL1zSjleinw/sy8GSAzbxhzjLM1\n0zo7BDh1LJHNzSjlSmDL+ngr4Joxxjdbo5RrN+Ab9fE3hyzvpIWSuIbdKmrHwZUi4m8i4mfAPwBH\njSm2uZi2XBHxOGCnzOzbXVVHqjPg2XUY44yI2GnI8q4ZpVwPBx4eEd+OiPPrX0Tog1HrjIjYBXgQ\naw6KXTZKuY4FDo2Iq4AvU3qTXTdKub4HHFQf/yWwRUTcfwyxzclCSVwjycz3Z+ZDgNcDb2o7nrmK\niPWA9wBHtx3LWvLvwJLMfDRwFnByy/HMlw0ow4V7UnolH46IrVuNaP4dDJyRmb9vO5B5cghwUmYu\nBvYFTqnfv757LbBHRFwM7EG5U1Hn62wh7HiY+a2iTgOetVYjmh/TlWsLYHfg7IhYCTwJOLMnEzSm\nrbPMvDEzf1effgR4/Jhim4tR2uJVwJmZeVdm/hz4MSWRdd1MvmcH049hQhitXIcDpwNk5nnAJpQb\n8HbZKN+xazLzoMx8LPDG+lrnJ9QslMQ17a2iIqJ5YNgP+MkY45utKcuVmbdm5naZuSQzl1AmZxyQ\nmX24eeYodfbAxtMDgMvHGN9sjXLbss9TeltExHaUocMrxhnkLI10S7aIeASwDXDemOObrVHKdSWw\nF0BEPJKSuFaPNcqZG+U7tl2j5/gG4GNjjnFWFkTiysy7gYlbRV0OnJ6Zl0bE2yLigLraK+rU40uA\n1wCHtRTuyEYsVy+NWLajap19j3JNclk70Y5uxHJ9FbgxIi6jXBB/XWbe2E7Eo5tBezwYOC3rVLWu\nG7FcRwMvrW3xVGBZ18s3Yrn2BH4UET8GtgeObyXYGfKWT5KkXlkQPS5J0rrDxCVJ6hUTlySpV0xc\nkqReMXFJknrFxCVJ6hUTlySpV/4/viFhCoPlBw0AAAAASUVORK5CYII=\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Onc-7QJUcqrw",
        "colab_type": "text"
      },
      "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",
      "metadata": {
        "id": "ZqeOJahwcqrx",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 35
        },
        "outputId": "f6fc6d5b-6456-45b4-e156-c0e2fb890c1e"
      },
      "source": [
        "prob_treatment_better = np.mean(posterior_samples_treatment['p'] > posterior_samples_control['p'] + 0.19)\n",
        "print(prob_treatment_better)"
      ],
      "execution_count": 15,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "0.94975\n"
          ],
          "name": "stdout"
        }
      ]
    }
  ]
}