{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bayesian estimation of vaccine effectiveness\n",
    "\n",
    "Copyright 2020 Allen Downey\n",
    "\n",
    "[MIT License](https://opensource.org/licenses/MIT)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# If we're running on Colab, install empiricaldist\n",
    "# https://pypi.org/project/empiricaldist/\n",
    "\n",
    "import sys\n",
    "IN_COLAB = 'google.colab' in sys.modules\n",
    "\n",
    "if IN_COLAB:\n",
    "    !pip install empiricaldist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from empiricaldist import Pmf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "According to [this press release](https://www.pfizer.com/news/press-release/press-release-detail/pfizer-and-biontech-announce-vaccine-candidate-against)\n",
    "\n",
    "> The first set of results from our Phase 3 COVID-19 vaccine trial provides the initial evidence of our vaccine’s ability to prevent COVID-19.\n",
    "\n",
    "The press release includes the following details about the results\n",
    "\n",
    "> The ... trial ... has enrolled 43,538 participants to date, 38,955 of whom have received a second dose of the vaccine candidate as of November 8, 2020.\n",
    ">\n",
    ">  ... the evaluable case count reached 94 and ... the case split between vaccinated individuals and those who received the placebo indicates a vaccine efficacy rate above 90%, at 7 days after the second dose.\n",
    "\n",
    "The press release provides only a point estimate for the effectiveness of the vaccine, and it does not provide enough information to make a better estimate.\n",
    "\n",
    "But with some guesswork, we can compute the posterior distribution of effectiveness, and use it to estimate the lower bound of the credible interval."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we don't know how many people are in each branch of the trial, I'll assume that it is approximately equal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_control = 38955 / 2\n",
    "n_treatment = 38955 / 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We know there were a total of 94 infections in the two branches.  Since the estimated effectiveness is 90%, I'll guess that there were 86 infections in the control branch and 8 in the treatment branch.\n",
    "\n",
    "We can make a beta distribution that represents the posterior distribution of the infection rate in the control branch, starting with a uniform distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.4446602437964785"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.stats import beta\n",
    "\n",
    "dist_control = beta(86+1, n_control+1)\n",
    "dist_control.mean() * 100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here's the posterior distribution for the treatment branch."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.046183450930083386"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dist_treatment = beta(8+1, n_treatment+1)\n",
    "dist_treatment.mean() * 100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The risk ratio is about 10:1, which is consistent with 90% effectiveness.\n",
    "\n",
    "To compute the distribution of risk ratios, I'll make a discrete approximation to the two posterior distributions, using the `Pmf` object from `empiricaldist`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_beta(dist):\n",
    "    \"\"\"PMF to approximate a beta distribution.\n",
    "    \n",
    "    dist: `beta` object\n",
    "    \n",
    "    returns: Pmf\n",
    "    \"\"\"\n",
    "    qs = np.linspace(8e-6, 0.008, 1000)\n",
    "    ps = dist.pdf(qs)\n",
    "    pmf = Pmf(ps, qs)\n",
    "    pmf.normalize()\n",
    "    return pmf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here are the `Pmf` objects:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "pmf_control = make_beta(dist_control)\n",
    "pmf_treatment = make_beta(dist_treatment)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here's what they look like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEGCAYAAABy53LJAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAvNUlEQVR4nO3deZwU1bn/8c8zPRvMMOw7KCAgOwgIGEFB1AjGLQmJxGiIUaJxuZqrV3OTGJObm2uM+/3lShSNeqOGmGjEG43GBeOKgCICiiyyjCD7NjPMTPfM+f1RNTgOs3V31fT0zPf9evWruqvqVD3dszx9zqk6x5xziIiINFZGqgMQEZH0osQhIiJxUeIQEZG4KHGIiEhclDhERCQumakOoCl06dLF9evXL9VhiIiklWXLlu1yznWtub5VJI5+/fqxdOnSVIchIpJWzGxTbevVVCUiInFR4hARkbgocYiISFxaRR+HiLQM0WiUwsJCSktLUx1Ki5Kbm0ufPn3Iyspq1P5KHCKSNgoLC2nXrh39+vXDzFIdTovgnGP37t0UFhbSv3//RpVRU5WIpI3S0lI6d+6spBEgM6Nz585x1eKUOEQkrShpBC/ez1SJIyyH9sLb90LpgVRHIiISKCWOsLz8S/j7jfDiz1IdiYgE7LPPPuP888/nmGOOYdiwYcycOZOPP/447uPcddddlJSUxF0uPz8/7jJBUuIIy5q/+8vnQJNlibQYzjnOO+88pk6dyvr161m9ejW/+tWv2L59e9zHqi9xVFRUJBtqaJQ4wlC6Hw4UQkEfOLjNe4hIi/DKK6+QlZXFZZdddnjdmDFjmDx5Mtdffz0jRoxg5MiRLFiwAIBFixYxdepUvv71rzNkyBAuuOACnHPcc889bN26lWnTpjFt2jTAq0ncdNNNTJw4kbfeeos77riDESNGMGLECO66665UvN1a6XLcMOxe7y2HzIR37oNda6GgV2pjEmlhfv7MKlZvDbYPcVivAn521vB691m5ciXjxo07Yv2TTz7J8uXLef/999m1axfHH388J510EgDvvfceq1atolevXpx44om88cYbXH311dxxxx288sordOnSBYDi4mJGjBjBL37xC5YtW8bvf/97Fi9ejHOOiRMncvLJJ3PccccF+p4ToRpHGKoSx+Av+6/XpS4WEWkSr7/+OrNnzyYSidC9e3dOPvlklixZAsCECRPo06cPGRkZjBkzho0bN9Z6jEgkwte+9rXDxzvvvPPIy8sjPz+fr371q7z22mtN9XbqpRpHGPasBwyOPhGy2ipxiISgoZpBWIYPH86f//znI9a7evoyc3JyDj+PRCLEYrFa98vNzSUSiTR4vFRTjSMM+7dAfnfIagMd+8OeT1IdkYgE5JRTTqGsrIz777//8LolS5bQsWNHFixYQEVFBTt37uSf//wnEyZMqPdY7dq14+DBg7VuO+mkk/jrX/9KSUkJxcXFPPXUU0yZMiXQ95Io1TjCULQT8v25Twp6wcGtqY1HRAJjZjz11FNcc8013HLLLeTm5tKvXz/uuusuioqKGD16NGbGrbfeSo8ePfjoo4/qPNbcuXOZMWMGPXv25JVXXvnCtrFjxzJnzpzDyeeSSy5pFv0bANacq0NBGT9+vGvSiZzumwptOsGFT8LCq7xLc69f23TnF2mhPvzwQ4YOHZrqMFqk2j5bM1vmnBtfc181VYWhaCfkd/Oet+sFxTuhIpramEREAqLEETTnoHgH5PlNVe16AA6K4r85SESkOVLiCFrpfqgo/7zGUXX/xgHdBCgiLYMSR9CKd3rLPD9x5Hf3lqpxiEgLEWriMLMzzGyNma0zsxtr2W5mdo+/fYWZjfXX9zWzV8zsQzNbZWb/Uq1MJzP7h5mt9Zcdw3wPcTu011u27ewt87w7QinZnZp4REQCFlriMLMI8FtgBjAMmG1mw2rsNgMY5D/mAvf662PAvzrnhgKTgCuqlb0ReMk5Nwh4yX/dfBza5y1z23vLNp28pRKHiLQQYd7HMQFY55zbAGBmfwTOAVZX2+cc4BHnXRP8tpl1MLOezrltwDYA59xBM/sQ6O2XPQeY6pd/GFgE3BDi+4hP6X5v2aaDt8xu6909rsQhkvZ2797N9OnTAW9o9UgkQteu3oUw77zzDtnZ2XEfc9GiRWRnZ/OlL30p0FjDPE+YiaM3sKXa60JgYiP26Y2fNADMrB9wHLDYX9XdTyw457aZWbfaTm5mc/FqMRx11FEJv4m4le7zllU1DvCarUr2NF0MIhKKzp07s3z5cgBuvvlm8vPzue666w5vj8ViZGbG92910aJF5OfnN0niCOo8YfZx1DYXYc27Devdx8zygb8A1zjn4hoG0zl3n3NuvHNufNU3giZxuKmqw+fr2naCkl1NF4OINJk5c+bwwx/+kGnTpnHDDTewfv16zjjjDMaNG8eUKVMO3zn+zDPPMHHiRI477jhOPfVUtm/fzsaNG5k3bx533nknY8aM4bXXXmPOnDlcfvnlTJs2jQEDBvDqq69y8cUXM3ToUObMmXP4vC+88AInnHACY8eOZdasWRQVFQHQr18/fvaznzF27FhGjhzJRx99VOt5khFmjaMQ6FvtdR+g5tgbde5jZll4SeNR59yT1fbZXtWcZWY9gR2BR56M0n1e01RmtSpr2y5qqhIJ2nM3wmcfBHvMHiNhxi1xF/v444958cUXiUQiTJ8+nXnz5jFo0CAWL17MD37wA15++WUmT57M22+/jZkxf/58br31Vm6//XYuu+yyL9RcHnjgAfbu3cvLL7/MwoULOeuss3jjjTeYP38+xx9/PMuXL6dPnz788pe/5MUXXyQvL49f//rX3HHHHdx0000AdOnShXfffZf/+Z//4bbbbmP+/PlHnCcZYSaOJcAgM+sPfAqcD3yrxj4LgSv9/o+JwH4/IRjwAPChc+6OWsp8B7jFXz4d4nuIX+m+LzZTgddUtWd9SsIRkfDNmjWLSCRCUVERb775JrNmzTq8raysDIDCwkK++c1vsm3bNsrLy+nfv3+dxzvrrLMwM0aOHEn37t0ZOXIk4I3Mu3HjRgoLC1m9ejUnnngiAOXl5ZxwwgmHy3/1q18FYNy4cTz55JNHniBJoSUO51zMzK4EngciwIPOuVVmdpm/fR7wLDATWAeUAN/1i58IXAh8YGbL/XX/7px7Fi9h/MnMvgdsBj7/CTUHh/Z9sZkK1MchEoYEagZhycvLA6CyspIOHToc7gep7qqrruKHP/whZ599NosWLeLmm2+u83hVw7BnZGR8YUj2jIwMYrEYkUiE0047jccff7ze8vUN4Z6MUEfH9f/RP1tj3bxqzx1wRS3lXqf2/g+cc7uB6cFGGqDS/bXXOMoOQKwMMnNqLyciaa+goID+/fvzxBNPMGvWLJxzrFixgtGjR7N//3569+4NwMMPP3y4TLt27ThwIL6ZDCdNmsQVV1zBunXrGDhwICUlJRQWFjJ48OA6yyRynrrozvGg1Zo4qu7lUK1DpKV79NFHeeCBBxg9ejTDhw/n6ae91vSbb76ZWbNmMWXKlMNTxYLXLPXUU0/F1WndtWtXHnroIWbPns2oUaOYNGlSvcO3J3qeumhY9aDdMxZ6jYGvP/j5utVPw58ugste9zrfRCQhGlY9PBpWPZXKiyA774vrqvo8qm4OFBFJY0ocQSsrgux2X1xX1XRVdY+HiEgaU+IIUmUlRIshJ/+L66uGH1GNQyRpraF5vanF+5kqcQQpWuwtj2iq8mscShwiScnNzWX37t1KHgFyzrF7925yc3MbXSbUy3FbnTLvln+ya9Q4ctoD9vk4ViKSkD59+lBYWMjOnTtTHUqLkpubS58+fRq9vxJHkMr9GkdOjT6OjAzIKVCNQyRJWVlZ9d5xLU1DTVVBKj/oLWs2VQG0aa/OcRFpEZQ4glRXUxV4/RyqcYhIC6DEEaTDTVW1JY4O6uMQkRZBiSNI5apxiEjLp8QRpLKqPo66ahxKHCKS/pQ4glRfU1WbDuocF5EWQYkjSA01VUWLoSLatDGJiARMiSNIZQchsw1kRI7cpoEORaSFUOIIUnlR7c1UoGFHRKTFUOIIUnlx7c1UUC1x7GuycEREwqDEEaSyoroTR9UIueogF5E0p8QRpGgJZLetfZuaqkSkhVDiCFKsFDLrGJr4cOf4vqaKRkQkFEocQYqWQFZdNY4Cb1l6oOniEREJgRJHkKKlkFVHjSOrLWRkqqlKRNKeEkeQooe8+zhqY+bNyVGmGoeIpDcljiDFDkFWHYkDNNChiLQIShxBijaUOArUxyEiaU+JIyjONSJxqMYhIulPiSMoFeWAq/tyXFAfh4i0CEocQYmWeMu6LscFzckhIi2CEkdQoqXesq7LcUF9HCLSIihxBKVRNY72UH4QKiuaJiYRkRAocQQl5tc4GurjAPVziEhaU+IISvSQt6y3xqFhR0Qk/SlxBOVw4qivj0Mj5IpI+lPiCEpjahxqqhKRFkCJIygxP3HU18ehGoeItABKHEE5XONoYMgRUB+HiKS1UBOHmZ1hZmvMbJ2Z3VjLdjOze/ztK8xsbLVtD5rZDjNbWaPMzWb2qZkt9x8zw3wPjdaoxNHBW6rGISJpLLTEYWYR4LfADGAYMNvMhtXYbQYwyH/MBe6ttu0h4Iw6Dn+nc26M/3g20MATFW1EU1VOO2+pPg4RSWNh1jgmAOuccxucc+XAH4FzauxzDvCI87wNdDCzngDOuX8Ce0KML1ixRnSOR7IgK081DhFJa2Emjt7AlmqvC/118e5Tmyv9pq0HzaxjbTuY2VwzW2pmS3fu3BlP3ImJlgIGmTn175dboMQhImktzMRhtaxzCexT073AMcAYYBtwe207Oefuc86Nd86N79q1awOHDEC0xOvfsNreUjUaWl1E0lyYiaMQ6FvtdR9gawL7fIFzbrtzrsI5Vwncj9cklnqx0vr7N6poaHURSXNhJo4lwCAz629m2cD5wMIa+ywELvKvrpoE7HfObavvoFV9IL7zgJV17dukoofq79+oohqHiKS5zLAO7JyLmdmVwPNABHjQObfKzC7zt88DngVmAuuAEuC7VeXN7HFgKtDFzAqBnznnHgBuNbMxeE1aG4Hvh/Ue4hI9VP9wI1VyC2DPhvDjEREJSWiJA8C/VPbZGuvmVXvugCvqKDu7jvUXBhljYBqaNraKahwikuZ053hQYocgsxGJQ30cIpLmlDiC0ugaR4E3P3nVjIEiImlGiSMo8TRVgZqrRCRtKXEEpbGJI8dPHGquEpE0pcQRlFhp4/o4VOMQkTSnxBGUaEnjL8cFJQ4RSVtKHEGJljb+BkBQ4hCRtKXEEQTn/MtxGznkCKiPQ0TSlhJHECrKwVXqqioRaRWUOILQmNn/qmTngUU0fayIpC0ljiDEkzjMNCeHiKQ1JY4gVM3+15jLcUHDjohIWlPiCEI8NQ7QQIciktaUOIJQNe5UXIlDNQ4RSU9KHEGIlnjLxiYONVWJSBpT4ghCzK9xNLaPQ01VIpLGlDiCEG+NI7dATVUikraUOIKQSB9H2QGorAwvJhGRkChxBKGqxtGYIUfAH3bEQfnB0EISEQmLEkcQYgnUOED9HCKSlpQ4gpBIHweon0NE0pISRxCq+jga21SlGoeIpDEljiDEDnmX4po1bn8NrS4iaUyJIwiNnW+8imocIpLG6k0cZvZQteffCT2adBUtTTBxqMYhIumnoRrH6GrP/yXMQNJatCS+xJGjecdFJH01lDhck0SR7mKljR9uBCAz29u/TIlDRNJPZgPb+5jZPYBVe36Yc+7q0CJLJ/HWOECTOYlI2moocVxf7fnSMANJa9FSyGrkpbhVNLS6iKSpehOHc+7hpgokrUVLoE2H+MpoaHURSVP1Jg4zW1jfdufc2cGGk6ZipY2/+a9Kbnso3RdKOCIiYWqoqeoEYAvwOLAYr69Daor3Pg7w+jj2bQ4nHhGREDWUOHoApwGzgW8BfwMed86tCjuwtJJQ4tBkTiKSnuq9HNc5V+Gc+7tz7jvAJGAdsMjMrmqS6NJF9FB8l+OC+jhEJG01VOPAzHKAM/FqHf2Ae4Anww0rjTjnjVWVSI0jVgqxMsjMCSc2EZEQNNQ5/jAwAngO+LlzbmWTRJVOKqLgKhO7HBe8S3LzuwYfl4hISBq6c/xCYDDecCNvmdkB/3HQzBpsZzGzM8xsjZmtM7Mba9luZnaPv32FmY2ttu1BM9thZitrlOlkZv8ws7X+smPj3mpIDs/F0Ta+chroUETSVEN9HBnOuXbVHgX+o51zrqC+smYWAX4LzACGAbPNbFiN3WYAg/zHXODeatseAs6o5dA3Ai855wYBL/mvUycW51wcVQ4Pra7EISLppaHRcXPN7Boz+39mNtfMGuwTqWYCsM45t8E5Vw78ETinxj7nAI84z9tABzPrCeCc+yewp5bjngNU3Zj4MHBuHDEFTzUOEWllGmqqehgYD3wAzARuj+PYvfHuAalS6K+Ld5+aujvntgH4y2617eQnuqVmtnTnzp1xhB2nqtn/4u7j0PSxIpKeGqpBDHPOjQQwsweAd+I4dm03C9Ycbbcx+yTEOXcfcB/A+PHjwxvlN3rIW8Zb49DQ6iKSphqqcUSrnjjnYnEeuxDoW+11H2BrAvvUtL2qOctf7ogzrmDF/MSRyJAjoHs5RCTtNDiRU/UrqYBRcVxVtQQYZGb9zSwbOB+oOfbVQuAi/+qqScD+qmaoeiwEqmYj/A7wdAP7hyvRGkd2PmBqqhKRtNPQ6LiRRA/snIuZ2ZXA80AEeNA5t8rMLvO3zwOexes7WQeUAN+tKm9mjwNTgS5mVgj8zDn3AHAL8Ccz+x6wGZiVaIyBOJw44qxxZGRoTg4RSUvxXCUVN+fcs3jJofq6edWeO+CKOsrOrmP9bmB6gGEmJ9EaB0BOezVViUjaaaipShqSaB8HaKBDEUlLShzJOlzjiHOsKvCbqlTjEJH0osSRrKQSh2ocIpJ+lDiSleiQI+APra7EISLpRYkjWdESby4OS2ByRNU4RCQNKXEkK1oa/6W4VXILoOwgVFYGG5OISIiUOJIVPZTYpbjg1ThcJZQXBRuTiEiIlDiSFTuUWP8GVBtaXVdWiUj6UOJIVlI1Dg10KCLpR4kjWdFDSfRxVJs+VkQkTShxJCuaRFPV4cSxL7BwRETCpsSRrGgJZOclVrZNJ29ZUttEhyIizZMSR7KS6eNo29lbHlLiEJH0ocSRrGhJ4okjpx1kZEHJ7mBjEhEJkRJHsqIliY1TBd7d5m07K3GISFpR4khW9FDiiQOgbSf1cYhIWlHiSEZlhTfIYaKd46Aah4ikHSWOZCQzpHoV1ThEJM0ocSQjmWljq6jGISJpRokjGdESb5lMjaNNJ+9yXI2QKyJpQokjGYcTR5I1Dlepu8dFJG0ocSQjqMQB6ucQkbShxJGMQDrHdfe4xC9WUUlJeQznXKpDkVYoM9UBpLVAOsc7ekt1kEsDymOV/OXdQhYs2cIHn+6notLROS+bKYO6cPHk/ozq0yHVIUorocSRjPJib5kdRFOVEofU7aPPDnD14+/x8fYihvUs4NIpAyhok8n6HcW8sOoz/rp8K9+edBQ/OXMYuVmRVIcrLZwSRzKCbKpS4pA6vLFuF5c+spS22ZnMv2g804d2w8wOby8qG86d//iYB9/4hA8K9/PAnOPpkp+TwoilpVMfRzKC6BzPzodIDhTvCiYmaVGWbtzDdx9aQt+ObXn26smcOqz7F5IGQH5OJj/9yjDmfXsca7Yf5NvzF7OvpDxFEUtroMSRjCD6OMwgvxsU7wwmJmkxNu0uZu7/LqN3hzY8PncS3QrqnzDsy8N7cP9F49mws5jvPrSE0mhFE0UqrY0SRzKCuAEQvMRRtD35eKTFKI9V8oNH36XSOR6cczyd8rIbVW7KoK7cff4Y3tu8j589vUpXXUkolDiSES3x5tOIZCV3nLxuULQjmJikRbj9H2tYtfUAt319NP27xDeI5oyRPbly2kAWLN3CE8sKQ4pQWjMljmQkM/tfdflKHPK5dzfv5b5/buBbE4/i1GHdEzrGtacNZtKATvzimdVs2VMScITS2ilxJKO8OLlLcavkd4eSXd4w7dKqVVQ6fvrXlXRrl8O/zxya8HEiGcZvvj4agOueeJ/KSjVZSXCUOJKR7CROVfK7eeNV6ZLcVu8Pb29i1dYD/PQrw8jPSe5q+b6d2vLTrwxl8Sd7eHzJ5oAiFFHiSE6QTVWgDvJWbl9JObe9sIbJA7tw5siegRzzG+P7MrF/J37z/Br2FOsSXQmGEkcykplvvLq8qsShfo7W7N5F6ykqi/GTrww94l6NRJkZ/3HuCIpKY/z6uY8COaaIEkcyoiUB1ziUOFqrz/aX8tCbGzlvTG+G9CgI9NiDu7fj4sn9WbB0C+9u3hvosaV1CjVxmNkZZrbGzNaZ2Y21bDczu8ffvsLMxjZU1sxuNrNPzWy5/5gZ5nuoV+CJQ01VrdXdL62l0jmuPW1wKMf/l+mD6Nouh//824e6t0OSFlriMLMI8FtgBjAMmG1mw2rsNgMY5D/mAvc2suydzrkx/uPZsN5Dg4LqHM/O9xKQ7h5vlTbvLuFPS7fwrQlH0bdTAF9EapGXk8kPTxvMsk17eX7VZ6GcQ1qPMGscE4B1zrkNzrly4I/AOTX2OQd4xHneBjqYWc9Glk29oDrHq4YdUY2jVZr/+gYyDC6fOjDU88wa14dB3fK55bmPKI9pqmJJXJiJozewpdrrQn9dY/ZpqOyVftPWg2bWsbaTm9lcM1tqZkt37gzpm3xQ93GAdy/HQX0TbG12F5Xxp6VbOHdMb3q0r38sqmRlRjL40cwhbNxdwmOLN4V6LmnZwkwctV0WUrNxta596it7L3AMMAbYBtxe28mdc/c558Y758Z37dq1UQHHLaimKoCCXnBgazDHkrTx8FubKI1W8v2TBzTJ+aYd240TBnTm7pfWcqA02iTnlJYnzMRRCPSt9roPUPM/Y1371FnWObfdOVfhnKsE7sdr1mp6sXKoKIPsdsEcr6A3HPgU1HHZapSUx3jkrY2cOrQ7A7sF9HvUADPj32cOZW9JlPte3dAk55SWJ8zEsQQYZGb9zSwbOB9YWGOfhcBF/tVVk4D9zrlt9ZX1+0CqnAesDPE91K28yFvm5AdzvILeECuFQ7pcsrVYsGQL+0qiXD61aWobVUb2ac9Zo3sx//UN7DhQ2qTnlpYhtMThnIsBVwLPAx8Cf3LOrTKzy8zsMn+3Z4ENwDq82sMP6ivrl7nVzD4wsxXANODasN5DvcoOesvsgBJHe78LZ79GM20NohWVzH/tE8Yf3ZFxR3dq8vP/62mDiVU47n5pbZOfW9JfqFPH+pfKPltj3bxqzx1wRWPL+usvDDjMxARe4+jjLQ98Cj1HBXNMabae/WAbn+47xM1nD0/J+ft1yeNbE4/i0cWb+d7k/gzoGtDvsbQKunM8UWV+4giqj0M1jlbDOce8VzcwsFs+04d0S1kcV50yiJzMDG57YU3KYpD0pMSRqHK/qSqoGkdeV8jI9Goc0qL9c+0uPtx2gLknDSAjI5gxqRLRtV0Ol04ZwLMffMbyLftSFoekHyWORB2ucQSUODIi0E6X5LYGv3t1Pd0LcjhnTK9Uh8KlJw2gc142tzynoUik8ZQ4EhV0Hwd4zVX7VeNoyVYU7uPN9bv53uT+5GRGUh0O+TmZXHXKQN7esIdXP9aQN9I4ShyJqqpx5AQ4kmlBbzigPo6W7HevbqBdbiazJxyV6lAO+9bEo+nbqQ2//vsazRQojaLEkaigL8cFaN/Hq3FoCtkWaeOuYp5buY1vTzqadrlZqQ7nsOzMDK47/Vg+3HaAhe+rqVQapsSRqPKDEMmGzOzgjtlpAFRGdWVVCzX/9Q1kZmTw3S/1S3UoRzhrVC+G9SzgthfWUBbTFxepnxJHosqKgq1tgJc4APZoKIiWZldRGU8sLeSrY3vTrSDcwQwTkZFh3DhjCIV7D/HYYs1PLvVT4khUeVGwHeOgxNGCPfTGRsorKrn0pKYdXiQeUwZ14cSBnfnvl9dxUAMgSj2UOBJVVhTczX9V2vWEzFwljhbmYGmUR97ayJeH9eCYZnyHtplxwxlD2FNczv2vfZLqcKQZU+JIVPnB4GscGRnQsT/s0R9tS/L4O5s5UBrj8qnHpDqUBo3q04EzR/Vk/msb2HFQAyBK7ZQ4EhVGHwd4zVWqcbQYZbEK5r/2CScO7Mzovh1SHU6jXHf6sZTHKvnvl9alOhRpppQ4ElV2EHJCmEOh8wDY+wlUamrPluDJdz9lx8EyLj853Glhg9S/Sx6zJxzF4+9sZuOu4lSHI82QEkeiykJoqgKvxhErhYO6nj7dVVQ6fvfqekb2bs+JAzunOpy4XDV9IFkRDYAotVPiSFTpPmhT63Tnyeky2Fvu+Cj4Y0uTem7lNjbuLuEHU4/BLHWDGSaiW7tcLp3Sn/9bsY0VhftSHY40M0ociYiWerWC3A7BH7u7Pz/D9tRMbCjBqKh03P3iWo7pmsfpw3ukOpyEXHrSALrkZ3PT06uo0FAkUo0SRyJK93nLNh2CP3abjt6kTttXNbyvNFv/t2Ira3cUcc2pg4mkcOj0ZLTLzeKnXxnG8i37eGzxplSHI82IEkciDu3zlmE0VYFX61CNI23FKiq5+8W1DOnRjjNH9kx1OEk5e3Qvpgzqwq1/X8N2zU8uPiWORFTVOMJoqgIvcez6GGJl4RxfQvX08q1s2FXMNacOTulETUEwM3557gjKKyr5xTOrUx2ONBNKHIk4XOPoEM7xuw+HypiXPCStlMUquPultQzvVcCXh3dPdTiBOLpzHldPH8TfPtjG31duS3U40gwocSTi0F5vGVZTVY+R3vKzD8I5voTmkTc3sXlPCf92xpC0u5KqPnNPGsDI3u350ZMfsENNVq2eEkciwm6q6jzQmyBqyzvhHF9CsbuojHteXsvUY7ty8uCuqQ4nUFmRDO785hgORSv4t7+s0DSzrZwSRyIO7QUMctuHc/yMCPQZr8SRZu56cS0l5RX85MyhqQ4lFAO75fPvM4eyaM1OHnlLV1m1ZkociSjaAW07e//gw9J3EuxY/Xl/ijRrKz/dz2PvbOaCiUcxsFsIQ9E0ExdOOppThnTjl39bzbJNe1MdjqSIEkciindCfrdwz9F3AuCgcEm455GkxSoqueEvK+iUl82/nnZsqsMJlZlx5zfG0KtDGy7/wzL1d7RSShyJKNoBeSG3YfedCJEcWP9yuOeRpD34xies2nqAX5w9nPZtm89c4mFp3zaL3104joOlMX7w6LuURjXVbGujxJGI4h3h1ziy20K/E2Hdi+GeR5KyfmcRd/zjY04f1p0zRqTn0CKJGNKjgN/MGsXSTXu5dsFyDUnSyihxJKJoJ+SFnDgABp7q3cuxVx2RzVFptIIrH3uPNlkR/uPcES3q8tvG+MqoXvzkzKE8t/Izbnp6pa60akWUOOJVVgTRYshvgsstj53hLVc/Hf65JG6/evZDPtx2gNu/MZruBbmpDiclLpkygO+fPIBHF2/m58+splI1j1ZBiSNeB/07Z/OboFmi0wDoNRZW/jn8c0lcFizZzCNvbeKSyf05ZUjLuEM8UTeeMYTvTe7PQ29u5EdPfqBmq1ZAiSNeVc1GHY9umvONnAXb3tdouc3Im+t28eOnVjJlUBdunDEk1eGknJnxkzOHcvUpA1mwdAuXPLyEA6XRVIclIVLiiNe+jd6yQxMljtHnQ2YbWDyvac4n9Vq2aQ+XPrKUAV3z+O0FY8mM6E8IvOTxw9OP5T/OHcFra3dx3m/fYN2Og6kOS0Ki3/p47d0EkWxo10TDZbftBGNmw/sLYP+nTXNOqdU7n+zhogfeoVtBLo9cPJGC3JZ/6W28Lpx0NH+4ZCJ7S6Kcec/r/P6NT9Tv0QIpccRr70Zo3xcymvCjO/EawMHLv2y6c8oX/HlZId+ev5juBbn8ce4kerRvnZ3hjTFpQGf+/i9T+NIxnfn5M6v55n1v8UHh/lSHJQFS4ojX9lXQtYnbtTseDRMvg/cf0w2BTay4LMZP/voB1z3xPuP7deQvl3+p1V5BFY9uBbk8OOd4bv3aKDbsLObs377OdU+8z4adRakOTQKgxBGP0v2wZz30Pq7pzz31R9DlWHhyLuxe3/Tnb2Wcczy/6jNm3P0ajy7ezKVT+vPwxRPomJed6tDShpnxjeP78sr1U7l0ygAWvr+V6Xe8yvf/dylvrNulq6/SmLWGm3bGjx/vli5dmvyBNiyCR86Bb//Fuzmvqe1aCw9+2essn/0Y9Bzd9DG0cNGKSl5YtZ35r2/gvc37OKZrHrd8bRTH9+uU6tDS3s6DZTz85kYeeWsjB0pj9CjI5azRPZl2bDfG9etITmaIg4ZKQsxsmXNu/BHrw0wcZnYGcDcQAeY7526psd387TOBEmCOc+7d+sqaWSdgAdAP2Ah8wzlX7zCdgSWOv/0rvPcoXL8WclI0Auq29+Gx86FkN3zpSph0BeR1Tk0sLcSh8gre3rCblz7azgurtrPjYBm9O7ThylMGMmtcH105FbDSaAX/WL2dp977lNfW7iRa4WibHWHc0R0Z1ac9o/p0YGiPAnp1yNVnn2JNnjjMLAJ8DJwGFAJLgNnOudXV9pkJXIWXOCYCdzvnJtZX1sxuBfY4524xsxuBjs65G+qLJZDEsb8Q7j0RjjkFZv0+uWMlq3gXPHeDd2NgRiYMmApHTYIeo6CgNxT08iaCimSmNs5mIFpRSXFZjIOlMYrLY+w6WM7W/YfYtq+UTbuLWbl1P+t2FFHpoG12hCmDuvDN4/ty8uBuRNJ8vvB0UFwW4631u3n1450s27SXNdsPHm7CyooYfTu25ajObeneLpcu7bLpnJdD5/xs2rfJIi8nk7bZEfKyM2mb4y2zIhlkZljaz/XeXKQicZwA3Oyc+7L/+kcAzrn/qrbP74BFzrnH/ddrgKl4tYlay1bt45zbZmY9/fL1jmWdaOJ46/c30GfLM2S5KJ3cXqJkclXbW9kUOQqAIz45V+/LI8byOXJ7zfKu3u39KjYxo2IRkyuXcLQ78lLdcjIpJZcysqnEqMQAo5IMnP+60u/mchz5hxZ+I2b9n0eSh8Nx5GdWXWaGkZOVQU5mhNysDNpkR8io5XNIqVY2/lWlc5TFKimPVRKtqKS8opJohaOiopJYnH0iZub9NA3/N58jPs/Gfrrp/GMoOf12hk36ckJl60ocYX4l7Q1sqfa6EK9W0dA+vRso2905tw3ATx61jjZoZnOBuQBHHXVUQm8gUtCD7W0HE7MsVkXa82bBTNrk9qP6NVU1f59qDnR35HYa2N7AL/YXVnTiXcbyLtCm4iDdyjbRIbaT9rGd5FQcItsdIruylKzKMgznP/z04SoxHBlUYtX/u8b5B9Lg7ke83/hOEO/xq8swIzPDf0QyiGQY2ZkZtMnyEkWk2f83aPn9jzVlAG38R03OeTXIspiXRCoqHbHKSm9Z4b2udA7noBJHpfO+rHlL74uYt6x+0KrFkZ+1O+JJev5E2rfND/yYYSaO2v4qa37ude3TmLL1cs7dB9wHXo0jnrJVJnztWuDaw6+nJ3KQJjU51QGIhMaAbP8hqRVmz1Mh0Lfa6z7A1kbuU1/Z7X4TFf5yR4Axi4hIA8JMHEuAQWbW38yygfOBhTX2WQhcZJ5JwH6/Gaq+sguB7/jPvwNozHERkSYUWlOVcy5mZlcCz+NdUvugc26VmV3mb58HPIt3RdU6vMtxv1tfWf/QtwB/MrPvAZuBWWG9BxEROZJuABQRkVrVdVWV7q4REZG4KHGIiEhclDhERCQuShwiIhKXVtE5bmY7gU0JFu8C7AownKAorvgorvgorvg017ggudiOds51rbmyVSSOZJjZ0tquKkg1xRUfxRUfxRWf5hoXhBObmqpERCQuShwiIhIXJY6G3ZfqAOqguOKjuOKjuOLTXOOCEGJTH4eIiMRFNQ4REYmLEoeIiMSlVSUOMzvDzNaY2Tp/vvKa283M7vG3rzCzsQ2VNbNOZvYPM1vrLzs2k7hmmdkqM6s0s4QuxQsprt+Y2Uf+/k+ZWYdmEtd/+PsuN7MXzKxXc4ir2vbrzMyZWZd44worNjO72cw+9T+z5WY2sznE5W+7yt+2ysxubQ5xmdmCap/VRjNb3kziGmNmb/txLTWzCQ0G4pxrFQ+84dnXAwPwJhF7HxhWY5+ZwHN4k41NAhY3VBa4FbjRf34j8OtmEtdQ4FhgETC+GX1epwOZ/vNfN6PPq6Ba+auBec0hLn97X7wpBjYBXZrRz/Jm4Lpm+Dc5DXgRyPFfd2sOcdUofztwU3OIC3gBmFGt/KKGYmlNNY4JwDrn3AbnXDnwR+CcGvucAzziPG8DHcybZbC+sucAD/vPHwbObQ5xOec+dM6tiTOWpojrBedczC//Nt7sjs0hrgPVyucR//TSYf1+AdwJ/FsCMTVFbMkIK67LgVucc2UAzrl4ZwkN9fMyMwO+ATzeTOJyQIH/vD1HztR6hNaUOHoDW6q9LvTXNWaf+sp2d96shfjLbs0krmQ1RVwX4307ahZxmdl/mtkW4ALgpuYQl5mdDXzqnHs/znhCj813pd8k8qDF30wbVlyDgSlmttjMXjWz45tJXFWmANudc2ubSVzXAL/xf/dvA37UUCCtKXFYLetqfoOra5/GlE1Uq4zLzH4MxIBHm0tczrkfO+f6+jFdmeq4zKwt8GPiT2Khx+Yv7wWOAcYA2/CaX5pDXJlAR7ymmuvxZgytbf+mjqvKbOKvbYQZ1+XAtf7v/rXAAw0F0poSRyFeW3GVPhxZJatrn/rKbvergvjLeKvFYcWVrNDiMrPvAF8BLnB+w2pziKuax4CvNYO4jgH6A++b2UZ//btm1qMZxIZzbrtzrsI5Vwncj9cckvK4/G1P+s017wCVeAP9pTouzCwT+CqwII54wo7rO8CT/vMnaMzPsTGdMi3hgfctZAPeH2JV59DwGvucyRc7lt5pqCzwG77YOX5rc4irWtlFJNY5HtbndQawGujazH6Og6qVvwr4c3OIq0b5jSTWOR7WZ9azWvlrgT82k7guA37hPx+M10RjqY6r2u//q83sd/9DYKr/fDqwrMFYEnkD6frAu2LgY7yrC35c7ZfsMv+5Ab/1t39AtX+4tZX113cGXgLW+stOzSSu8/C+ZZQB24Hnm0lc6/w/5OX+I66rl0KM6y/ASmAF8AzQuznEVeP4G0kgcYT4mf2vv+8KYCHVEkmK48oG/uD/PN8FTmkOcfnbHqo6RjP6OU4GluElk8XAuIbi0JAjIiISl9bUxyEiIgFQ4hARkbgocYiISFyUOEREJC5KHCIiEhclDmnVzKyoEftM8UdZXW5mbeI8/rlmNqza61+Y2amJxJoIM+tgZj9oqvNJ66DEIdKwC4DbnHNjnHOH4ix7LnA4cTjnbnLOvRhkcP7dyHXpAChxSKCUOEQAM5tqZovM7M/mzRfyqD+3wSV4I5neZGaP+vteb2ZL/MH9fl7tGBf56943s/81sy8BZ+MNILfczI4xs4fM7Ov+/tPN7D0z+8AfJDDHX7/RzH5uZu/624bUEu8cM3vCzJ4BXjCzfDN7qVqZqpFPbwGO8c//m/riF2ms+r6piLQ2xwHD8cbweQM40Tk338wmA//nnPuzmZ0ODMIbz8eAhWZ2ErAbb0DCE51zu8ysk3Nuj5ktrCoLUDXWnpnl4t1FPN0597GZPYI32Nxdfiy7nHNj/Wam64BLaon3BGCUf55M4Dzn3AHzJnt62z/3jcAI59wY/7y1xu+c+2cwH6G0BqpxiHzuHedcofMG7VsO9Ktln9P9x3t4w1kMwftHfAre+Fa7AJxzexo417HAJ865j/3XDwMnVdteNejcsjriAPhHtfMY8CszW4E3iVFvoHsc8Ys0mmocIp8rq/a8gtr/Pgz4L+fc776w0uxq4hvSvqFhvqtiqSsOgOJqzy8AuuKNMxT1R9PNreO8R8QvEg/VOETi8zxwsZnlA5hZbzPrhjfA5TfMrLO/vpO//0GgXS3H+QjoZ2YD/dcXAq8mEVd7YIefNKYBR9dx/rriF2k01ThE4uCce8HMhgJv+f0VRcC3nXOrzOw/gVfNrAKvKWgO3hSd9/s1kq9XO06pmX0XeMLvn1gCzEsitEeBZ8xsKV4z20f+eXab2RtmthJ4zjl3fW3xE/88MtKKaXRcERGJi5qqREQkLkocIiISFyUOERGJixKHiIjERYlDRETiosQhIiJxUeIQEZG4/H+zbnFQLC6BSwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pmf_control.plot(label='Control')\n",
    "pmf_treatment.plot(label='Treatment')\n",
    "\n",
    "plt.xlabel('Infection rate')\n",
    "plt.ylabel('PMF')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, it looks like the infection rate is about 10 times higher in the control group.\n",
    "\n",
    "We can use `div_dist` to compute the risk ratio."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "pmf_ratio = pmf_treatment.div_dist(pmf_control)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's the CDF of the risk ratio.  I cut it off at 1 because higher values have very low probabilities; that is, we are pretty sure the treatment is effective."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAchklEQVR4nO3de5Bc5Xnn8e8zN91HtxldGF1GEroggQSyJDAxNmCyIJKYOPauwU7ssLZZEhMTb20WlqpNtsq1tWaJa7PEGJXiEExVglIL2FG8Mph1uDgWGI0MkpDQiLkIaZA0N13mJmnU3c/+0S3UDKPu1kyfPj19fp+qqelz6Z6nX436N+e873mPuTsiIiIXUxZ2ASIiUtwUFCIikpGCQkREMlJQiIhIRgoKERHJqCLsAi5VTU2N19fXh12GiMiYsnPnzi53rx3Jc8dcUNTX19PQ0BB2GSIiY4qZvTfS5+rUk4iIZKSgEBGRjBQUIiKSkYJCREQyUlCIiEhGgQWFmT1hZh1m9vZFtpuZPWpmTWa228zWBlWLiIiMXJBHFE8Ct2XYvhFYmvq6B3g8wFpERGSEAruOwt1fNbP6DLvcATzlyXnOXzezaWY2192PBlWTJMUTTs/pc/SdjTEwGGdgMMbpwThn4wnOxRKcizuD8TjnYk4skXkaeif7NPXZZrLP+go5TIWfbY+sNRTgZ+T2GqOf9j97e2f5N83L+xhdDbnWkf01RvdecymhGNo7L42VQZgX3NUBh9OW21LrPhIUZnYPyaMOFixYUJDixrIT/YM0tvfS3NnHoe4B2k6epv3UGTr7znK8f5DeM7GwSxSRPDML7rXDDIrh3tawsejum4HNAOvWrdOdloZ4/+RpXmns5KEf7fnItqqKMuqmTWB29ThWz5vGzElVTJtYSfX4SiaPr2BSVQUTq8oZX1nOuMoyqsrLqKooo7K8jMpyo6KsLOsvYE6/n1lfI/MOufwnyLaLZXmRXN5H9rYYfaH5aO/Rvtfc2nv0/2ZZf0Ye2nu07ZmtLXN7jdH/jNGy74z8uWEGRRswP215HnAkpFrGnFMD59i6+wjP7Gxj1+GTH9r2XzauYMXcai6fNZm51eMpKwv+l1BESleYQbEVuM/MtgDXAqfUP5HdsVNnePzlJv6x4TBnziVYMWcKD25cwS1XzGJJ7eSC/GUiItESWFCY2dPAjUCNmbUBfwFUArj7JmAbcDvQBAwAdwdVSyk4G4uz6eUWHn+liVjc+ew1dXzl+npWXVatcBCRQAU56umuLNsd+EZQP7+U7DvSw/1b3uTdjj5+66q5PHDbChbMnBh2WSISEWNumvGo+fGb7/PAs7upnlDJ3929npuWzwq7JBGJGAVFEXv85WYefn4/1y6awWNfWkvN5HFhlyQiEaSgKFJ//fN3+e6LB/jMmsv4y3+7hqoKTcslIuFQUBShp147yHdfPMDvra3jkc+voVzDW0UkRPoztcj8sqmL/7Z1L7dcMZv/+bnVCgkRCZ2Cooj0nDnHn/2fXSyqmcSjd11NRbn+eUQkfDr1VET+x7Z3ONZzhmf/6HomVumfRkSKg/5kLRKvHujk6TcO8/VPLuaaBdPDLkdE5AMKiiLQc+YcDz67myW1k/jWLcvCLkdE5EN0fqMInD/l9MwfXc/4yvKwyxER+RAdUYRs53vHk6ecbljMWp1yEpEipKAIkbvz8E8bmTVlHPffsjTsckREhqWgCNEvm7p54+Bx7rv5co1yEpGipaAI0fdfbmJ29Ti+sH5+9p1FREKioAjJW4dPsr25m699YjHjKtSBLSLFS0ERkr/5RQvV4yu469oFYZciIpKRgiIEx06d4fm3j3HnhgVMHqe+CREpbgqKEDz9xiES7vz+tQvDLkVEJCsFRYHFE84/7jjMJ5fW6namIjImKCgK7NV3OznWc4a7Nmikk4iMDQqKAntmZxszJlVx84rZYZciIpITBUUB9Z45x//b187vrJ6rW5uKyJihT6sCemFvO2djCT5zdV3YpYiI5ExBUUDb9hylbtoE1i6YFnYpIiI5U1AUSN/ZGP/a1MWtq+Zgpvtgi8jYoaAokFcPdDIYS3DrKnVii8jYoqAokJ/tPcaMSVV8bKHuOSEiY4uCogAGYwl+vr+DT6+YRUW5mlxExhZ9ahXAr1q76T0T49ZVc8IuRUTkkikoCuDFfe1MqCznE0trwi5FROSSKSgC5u681NjB9UtmMr5S950QkbEn0KAws9vMrNHMmszswWG2TzWzfzazXWa218zuDrKeMLzXPcDh46f51PLasEsRERmRwILCzMqBx4CNwErgLjNbOWS3bwD73H0NcCPwXTOrCqqmMGxv7gbgNy7XaScRGZuCPKLYADS5e4u7DwJbgDuG7OPAFEtegTYZOA7EAqyp4LY3dzG7ehyLayaFXYqIyIgEGRR1wOG05bbUunTfA64AjgB7gPvdPTH0hczsHjNrMLOGzs7OoOrNu0TCea25m+uX1OhqbBEZs4IMiuE+GX3I8q3AW8BlwNXA98ys+iNPct/s7uvcfV1t7dg513+go5fu/kE+vmRm2KWIiIxYkEHRBqTfnWceySOHdHcDz3lSE9AKrAiwpoLa3pTsn7heQSEiY1iQQbEDWGpmi1Id1HcCW4fscwj4NICZzQaWAy0B1lRQ25u7WThzIvOm65anIjJ2VQT1wu4eM7P7gBeAcuAJd99rZvemtm8Cvg08aWZ7SJ6qesDdu4KqqZBi8QS/aunmt9fMDbsUEZFRCSwoANx9G7BtyLpNaY+PAP8myBrCsvdID71nY1y3WKedRGRs05XZAdlx8DgA1y5SUIjI2KagCEjDwRPMnzGBOVPHh12KiMioKCgC4O40vHeC9QtnhF2KiMioKSgCcPj4abr6znKNblIkIiVAQRGAnYeS/RMfW6CgEJGxT0ERgF2HTzGhspzlc6aEXYqIyKgpKALw5uGTXDVvKuVlmt9JRMY+BUWenY3FeedID9csmBZ2KSIieaGgyLO9R3oYjCe4Zv60sEsREckLBUWe7Tp8EoA1CgoRKREKijzb03aKWVPGMXfqhLBLERHJCwVFnu1qO8nqeVPDLkNEJG8UFHnUdzZGS1c/V9YpKESkdCgo8mjv+6dwR0cUIlJSFBR5tOf9UwBcVTct3EJERPJIQZFH+472UDtlHLVTxoVdiohI3igo8mjv+z1ceVl12GWIiOSVgiJPzpyL09TZx6rL1D8hIqVFQZEnjcd6iSeclTqiEJESo6DIk31HewBYpaAQkRKjoMiT/Ud7mFRVzvzpE8MuRUQkrxQUebLvaA8r5lZTpqnFRaTEKCjywN3Zf6yXK+bqRkUiUnoUFHlw9NQZes/EWD5H/RMiUnoUFHnQeKwXgBW69amIlCAFRR7sTwXFstkKChEpPQqKPDjQ3svcqeOZOqEy7FJERPJOQZEH+4/1slynnUSkRCkoRikWT9Dc0cdynXYSkRKloBilg939DMYTOqIQkZIVaFCY2W1m1mhmTWb24EX2udHM3jKzvWb2SpD1BEEd2SJS6iqCemEzKwceA34TaAN2mNlWd9+Xts804PvAbe5+yMxmBVVPUA4c66XM4PJZk8MuRUQkEEEeUWwAmty9xd0HgS3AHUP2+SLwnLsfAnD3jgDrCURjey/1MycxvrI87FJERAIRZFDUAYfTlttS69ItA6ab2ctmttPMvjzcC5nZPWbWYGYNnZ2dAZU7Mu+297F0to4mRKR0BRkUw82O50OWK4CPAb8F3Ar8VzNb9pEnuW9293Xuvq62tjb/lY7Q2Vicg939GvEkIiUtsD4KkkcQ89OW5wFHhtmny937gX4zexVYAxwIsK68ae3qJ+GwRP0TIlLCgjyi2AEsNbNFZlYF3AlsHbLPPwE3mFmFmU0ErgXeCbCmvGru6Adg6SwdUYhI6QrsiMLdY2Z2H/ACUA484e57zeze1PZN7v6OmT0P7AYSwA/c/e2gasq35s4+zGBRzaSwSxERCUyQp55w923AtiHrNg1ZfgR4JMg6gtLc2UfdtAlMqNKIJxEpXboyexSaO/tYUqv+CREpbQqKEUoknOaOfgWFiJQ8BcUIHes5w+lzcZbMUv+EiJQ2BcUINXf2AeiIQkRKnoJihJo7FBQiEg0KihFq7uynenwFNZOrwi5FRCRQCooRau7sY3HtZMyGm6lERKR0ZAwKM3sy7fFXAq9mDNHQWBGJimxHFGvSHt8fZCFjSe+Zc7T3nNWIJxGJhGxBMXS2VyE5GSCoI1tEoiHbFB7zzOxRklOGn3/8AXf/ZmCVFTENjRWRKMkWFH+W9rghyELGkuaOfirKjIUzJ4ZdiohI4DIGhbv/sFCFjCXNnX0smDmRynINGhOR0pf1k87MvmJmvzaz/tRXw8VuWRoVGvEkIlGS8YgiFQh/CvxH4Nck+yrWAo+YGe7+VOAVFplYPMHBrgFuXjE77FJERAoi2xHFHwOfdfeX3P2Uu590938BPpfaFjltJ04zGE+wpFZDY0UkGrIFRbW7Hxy6MrWuOoiCit0HI550n2wRiYhsQXF6hNtK1gdBUaOgEJFoyDY89goz2z3MegMWB1BP0Wvt6mfmpCqmTqwMuxQRkYLIFhRrgNnA4SHrFwJHAqmoyLV09lNfo/4JEYmObKee/hfQ4+7vpX8BA6ltkdPa1c8iBYWIREi2oKh394+cenL3BqA+kIqKWN/ZGB29ZxUUIhIp2YJifIZtE/JZyFhw8IPJABUUIhId2YJih5l9fehKM/sqsDOYkopXSyooFmnEk4hESLbO7D8FfmRmX+JCMKwDqoDPBlhXUWrt7McMTQYoIpGSbVLAduB6M7sJuDK1+v+mrs6OnNauPuZWj2d8ZXnYpYiIFEy2IwoA3P0l4KWAayl6rV39LFL/hIhEjObJzpG7a2isiESSgiJHx/sH6TkTo36mgkJEokVBkaOD3edHPCkoRCRaFBQ5aulUUIhINAUaFGZ2m5k1mlmTmT2YYb/1ZhY3s88HWc9oHOzup7zMmD9DQ2NFJFoCCwozKwceAzYCK4G7zGzlRfZ7GHghqFryobWrnwUzdJ9sEYmeID/1NgBN7t7i7oPAFuCOYfb7E+BZoCPAWkattWuAel1oJyIRFGRQ1PHh6cnbUus+YGZ1JK/w3pTphczsHjNrMLOGzs7OvBeajbtzsEvTi4tINAUZFDbMOh+y/FfAA+4ez/RC7r7Z3de5+7ra2tp81Zezzt6znD4XZ7GCQkQiKKcrs0eoDZiftjyPj97saB2wxcwAaoDbzSzm7j8OsK5Ldn4ywIW6hkJEIijIoNgBLDWzRcD7wJ3AF9N3cPdF5x+b2ZPAT4otJCDZkQ0aGisi0RRYULh7zMzuIzmaqRx4wt33mtm9qe0Z+yWKSWtXP1UVZVw2LXK34BARCfSIAnffBmwbsm7YgHD3PwyyltFo6exn4YyJlJcN1+0iIlLadFFADg52azJAEYkuBUUW8YRzqHtAQSEikaWgyOLIydMMxhMKChGJLAVFFi0a8SQiEaegyKK1sw9Ad7YTkchSUGTR2tXP5HEV1E4eF3YpIiKhUFBk0do9QH3NRFJXj4uIRI6CIouDXf0sqpkcdhkiIqFRUGRwNhan7YSGxopItCkoMjjUPUDCYYk6skUkwhQUGWhorIiIgiKjlk4FhYiIgiKD1q4+aqeMY8r4yrBLEREJjYIig5ZOTQYoIqKgyKClq1+3PxWRyFNQXMSJ/kGO9w+yWCOeRCTiFBQX0dKVnONpSa0uthORaFNQXERzR3LEk4JCRKJOQXERzV19VJWXMW+67pMtItGmoLiI5o5+Fs6cSEW5mkhEok2fghfR2tWnjmwRERQUwzoXT/Be9wCL1T8hIqKgGM6h4wPEEq6ObBERFBTDaupIDo29fJaCQkREQTGM80Gh6cVFRBQUw2ru7GN2tSYDFBEBBcWwmjv71T8hIpKioBgikXDebe9l2ewpYZciIlIUFBRDHDl1moHBuDqyRURSFBRDnO/IXqqgEBEBAg4KM7vNzBrNrMnMHhxm+5fMbHfqa7uZrQmynlwcaO8F0KknEZGUwILCzMqBx4CNwErgLjNbOWS3VuBT7r4a+DawOah6cnWgPXn70+mTqsIuRUSkKAR5RLEBaHL3FncfBLYAd6Tv4O7b3f1EavF1YF6A9eTkQHsvy2brtJOIyHlBBkUdcDhtuS217mK+Cvx0uA1mdo+ZNZhZQ2dnZx5L/LDkiKc+nXYSEUkTZFDYMOt82B3NbiIZFA8Mt93dN7v7OndfV1tbm8cSP6ztxGlOn4srKERE0lQE+NptwPy05XnAkaE7mdlq4AfARnfvDrCerN451gPAijkKChGR84I8otgBLDWzRWZWBdwJbE3fwcwWAM8Bf+DuBwKsJSf7j/ZiphFPIiLpAjuicPeYmd0HvACUA0+4+14zuze1fRPw58BM4PtmBhBz93VB1ZRNY3sPC2ZMZNK4IA+0RETGlkA/Ed19G7BtyLpNaY+/BnwtyBouxf5jvSzX0YSIyIfoyuyUgcEYrV39XDG3OuxSRESKioIiZf+xXtxh5WUKChGRdAqKlH1HkiOeVuqIQkTkQxQUKXuP9DB1QiXzpk8IuxQRkaKioEjZd+QUK+dWkxp9JSIiKQoK4GwszjtHe1k9f2rYpYiIFB0FBcn+icF4gmvmTwu7FBGRoqOgAN46fBKAq+dPD7cQEZEipKAgGRRzqsczZ+r4sEsRESk6CgqSQXG1TjuJiAwr8kFxvH+Q97oHuHrBtLBLEREpSpEPil0f9E9MC7UOEZFiFfmgePPwScoMrqrT0FgRkeEoKA6dYNnsKZpaXETkIiIdFOfiCX793gnW188IuxQRkaIV6aDY8/4p+gfjXLd4ZtiliIgUrUgHxWvNyVt0X7dYRxQiIhcT6aD4ZVMXK+ZMYebkcWGXIiJStCIbFAODMRoOnuCGpTVhlyIiUtQiGxTbm7oZjCf41LJZYZciIlLUIhsUP9/fzqSqcjYsUv+EiEgmkQyKeMJ5cV8HNy6fRVVFJJtARCRnkfyUfKP1OF19Z9l41ZywSxERKXqRDIqtu44wobKcm1eof0JEJJvIBcWZc3F+svsIt66azcQqTdshIpJN5ILiJ7uP0nsmxr9bPz/sUkRExoRIBYW788S/trJ01mQ+rmk7RERyEqmgeKmxg31He/j6DYsxs7DLEREZEyITFLF4god/2siCGRP53Wvqwi5HRGTMiExQ/M0vWmls7+Wh21fo2gkRkUsQiU/M15q7+e7PGtl45RxuXaVrJ0RELkWgQWFmt5lZo5k1mdmDw2w3M3s0tX23ma3Ndw0/f6edr/5wB/U1k/jO51arb0JE5BIFdiGBmZUDjwG/CbQBO8xsq7vvS9ttI7A09XUt8Hjq+4glEk5LVz9vHjrBi/va+dm+dlZdVs3f/eF6pk6oHM1Li4hEUpBXnG0Amty9BcDMtgB3AOlBcQfwlLs78LqZTTOzue5+9GIveqC9l5v+8mViiQSJBMQSCeIJiCcSxBPO2ViCs7EEANMnVvKtW5bxHz61mPGV5YG9URGRUhZkUNQBh9OW2/jo0cJw+9QBHwoKM7sHuAeg+rLFXFk3lYoyo8ws+b0s+b089X3p7MmsXTCdJbWTKSvTqSYRkdEIMiiG+4T2EeyDu28GNgOsW7fO//qua0ZfnYiI5CTIzuw2IH2ejHnAkRHsIyIiIQoyKHYAS81skZlVAXcCW4fssxX4cmr003XAqUz9EyIiUniBnXpy95iZ3Qe8AJQDT7j7XjO7N7V9E7ANuB1oAgaAu4OqR0RERibQebbdfRvJMEhftyntsQPfCLIGEREZnUhcmS0iIiOnoBARkYwUFCIikpGCQkREMrJkf/LYYWa9QGPYdRSJGqAr7CKKhNriArXFBWqLC5a7+5SRPDHQUU8BaXT3dWEXUQzMrEFtkaS2uEBtcYHa4gIzaxjpc3XqSUREMlJQiIhIRmMxKDaHXUARUVtcoLa4QG1xgdrighG3xZjrzBYRkcIai0cUIiJSQAoKERHJqGiDwsxuM7NGM2sysweH2W5m9mhq+24zWxtGnYWQQ1t8KdUGu81su5mtCaPOQsjWFmn7rTezuJl9vpD1FVIubWFmN5rZW2a218xeKXSNhZLD/5GpZvbPZrYr1RYlOVO1mT1hZh1m9vZFto/sc9Pdi+6L5LTkzcBioArYBawcss/twE9J3iXvOuBXYdcdYltcD0xPPd4Y5bZI2+9fSM5c/Pmw6w7x92IayXvUL0gtzwq77hDb4iHg4dTjWuA4UBV27QG0xSeBtcDbF9k+os/NYj2i2AA0uXuLuw8CW4A7huxzB/CUJ70OTDOzuYUutACytoW7b3f3E6nF10neKbAU5fJ7AfAnwLNARyGLK7Bc2uKLwHPufgjA3Uu1PXJpCwemmJkBk0kGRaywZQbP3V8l+d4uZkSfm8UaFHXA4bTlttS6S92nFFzq+/wqyb8YSlHWtjCzOuCzwCZKWy6/F8uA6Wb2spntNLMvF6y6wsqlLb4HXEHyVst7gPvdPVGY8orKiD43i3UKDxtm3dBxvLnsUwpyfp9mdhPJoPhEoBWFJ5e2+CvgAXePJ/94LFm5tEUF8DHg08AE4DUze93dDwRdXIHl0ha3Am8BNwNLgBfN7Bfu3hNwbcVmRJ+bxRoUbcD8tOV5JP8SuNR9SkFO79PMVgM/ADa6e3eBaiu0XNpiHbAlFRI1wO1mFnP3HxekwsLJ9f9Il7v3A/1m9iqwBii1oMilLe4GvuPJE/VNZtYKrADeKEyJRWNEn5vFeuppB7DUzBaZWRVwJ7B1yD5bgS+nevGvA065+9FCF1oAWdvCzBYAzwF/UIJ/LabL2hbuvsjd6929HngG+OMSDAnI7f/IPwE3mFmFmU0ErgXeKXCdhZBLWxwieWSFmc0GlgMtBa2yOIzoc7MojyjcPWZm9wEvkBzR8IS77zWze1PbN5Ec0XI70AQMkPyLoeTk2BZ/DswEvp/6SzrmJThjZo5tEQm5tIW7v2NmzwO7gQTwA3cfdtjkWJbj78W3gSfNbA/J0y8PuHvJTT9uZk8DNwI1ZtYG/AVQCaP73NQUHiIiklGxnnoSEZEioaAQEZGMFBQiIpKRgkJERDJSUIiISEYKComc1Kyyb5nZ26kZRael1l9mZs9keF79xWblHGEdDw1Z3p6v1xbJJw2Plcgxsz53n5x6/EPggLv/9xyeVw/8xN2vzPHnlLt7PJc6RIqZjigk6l4jNSla+hGDma0yszdSRx67zWxp+pPMbLGZvWlm64esv9HMXjKzfyA5+Rxm9uPUpHx7zeye1LrvABNSr//3qXV9qe9mZo+kjnj2mNkXAm4DkYyK8spskUIws3KS0zr87TCb7wX+t7v/fWpaiHJgdup5y0lOZX23u781zHM3AFe6e2tq+d+7+3EzmwDsMLNn3f1BM7vP3a8e5vm/B1xNcl6mmtRzXi3RKWpkDNARhUTRBDN7C+gGZgAvDrPPa8BDZvYAsNDdT6fW15KcQ+n3LxISAG+khQTAN81sF8l7hcwHlg7/tA98Anja3ePu3g68AqzP8hyRwCgoJIpOp/6SX0jyjmjfGLqDu/8D8BngNPCCmd2c2nSK5Hz+v5Hh9fvPPzCzG4FbgI+7+xrgTWB8lvpKen50GXsUFBJZ7n4K+Cbwn8ysMn2bmS0GWtz9UZIzbq5ObRoEfpfkDJxfzOHHTAVOuPuAma0gefvJ884N/bkprwJfMLNyM6sleXvLqE2HLUVEQSGR5u5vkrzH8p1DNn0BeDt1imoF8FTac/qB3wa+ZWbD3Yo13fNAhZntJjmD6etp2zYDu893Zqf5EckZX3eRvPf3f3b3Y5fyvkTyScNjRUQkIx1RiIhIRgoKERHJSEEhIiIZKShERCQjBYWIiGSkoBARkYwUFCIiktH/B+EI1F35qCJSAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pmf_ratio.make_cdf().plot()\n",
    "plt.xlim([0, 1])\n",
    "\n",
    "plt.xlabel('Risk ratio')\n",
    "plt.ylabel('CDF');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The median of the risk ratio is about 0.10.  Again, that's consistent with an effectiveness of 90%."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(0.10040984)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pmf_ratio.median()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To compute the distribution of effectiveness, we have to compute the distribution of `1-RR`, where `RR` is the risk ratio.  We can do that with `empiricaldist` by creating a deterministic `Pmf` with the quantity `1` and using `sub_dist` to subtract two `Pmf`s."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "effectiveness = Pmf.from_seq(1).sub_dist(pmf_ratio)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's the result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAdJ0lEQVR4nO3deXRcZ5nn8e8jyZL3XXa8y/uSYIdEJCFhcSANcaBJ03BIgBOWA5MJTVim58wkp08Df9BnJplMA8MQ2u2mQ5rTQLoPSYMZHBxoAgHSDpaJ49hO5MiSF3mRSrKt1ba2Z/6oq7iilKrKkm7dW6Xf55w6qrvUradey/XTvfe97zV3R0REZDglURcgIiLxpqAQEZGMFBQiIpKRgkJERDJSUIiISEZlURdwuebOnetVVVVRlyEiUlD27NnT4u6VI3ltwQVFVVUVNTU1UZchIlJQzOzoSF+rQ08iIpKRgkJERDJSUIiISEYKChERyUhBISIiGYUWFGb2iJk1m9n+YZabmX3TzOrMbJ+ZXRNWLSIiMnJh7lE8CtyaYfkWYHXwuBv4uxBrERGREQrtOgp3f8bMqjKscjvwPU+Oc77LzGaa2QJ3PxVWTSIixaStu5cXT7TReLabjgt9DLgz4DDgjqc8Hxjl3SSivOBuEXA8ZboxmPe6oDCzu0nudbB06dK8FCciEld7jp7h208f5jeHEvSNNgVyEGVQWJp5aT+xu28DtgFUV1frTksiMi65O3/71CG+9XQdldMq+NRblvO2NZUsnT2ZGZMnUFZilJhhBiVmwQPMDHtw5O8bZVA0AktSphcDJyOqRUQk9r7573V86+k67qhewlfet4HJ5fn5Co+ye+x24GNB76cbgDadnxARSe8XB5v4+i8P8efXLOKBD7whbyEBIe5RmNkPgc3AXDNrBL4CTABw963ADuA2oA7oBj4ZVi0iIoUs0XGR//qve9m4eAb/4/1vwCzdkfvwhNnr6cNZljvw2bDeX0SkWHztF4c439vP1++4mokTSvP+/royW0Qkxo6f6eZfa47zkeuWsrJyaiQ1KChERGLsH3/XQInBZzaviqwGBYWISEx1XezjR3saee/GhVwxY2JkdSgoRERi6mf7TtF5sY+PXh/thcYKChGRmHri+UZWzJ3CtctmRVqHgkJEJIaa2y/wXMMZ3nf1wrx3hx1KQSEiEkM/P3Aad3jvxgVRl6KgEBGJo18cbGLF3Cmsmjct6lIUFCIicdN1sY9d9a28c/28qEsBFBQiIrHz7OFWevudm9cpKEREJI3fvpJgcnlp5L2dBikoRERi5nd1LVy3fDYVZfkf1ykdBYWISIw0t1+gPtHFjSvnRF3KqxQUIiIxsqvhDAA3rFBQiIhIGrsbzjClvJQNC6ZHXcqrFBQiIjGy+8gZrlk2i7LS+Hw9x6cSEZFxruNCL7VNHVyzNB69nQYpKEREYmLv8XO4E5tusYMUFCIiMbH32DnM4OqlM6Mu5TUUFCIiMbH3+DlWVU5l+sQJUZfyGgoKEZEYcHdeaDzHxsUzoy7ldRQUIiIxcOLceVo6e7h6yYyoS3kdBYWISAwcONkOwJWLFBQiIpLG/hNtlJZYrC60G6SgEBGJgf0n2lhVOZWJE+IxEGAqBYWISAwcPNXOhoXx25sABYWISORaOi/S1H6RKxUUIiKSzkunkiey43h+AhQUIiKROxj0eFqvoBARkXRqT3dwxfSJzJpSHnUpaSkoREQidvBUO2uvmBZ1GcMKNSjM7FYzqzWzOjO7P83yGWb2UzN7wcwOmNknw6xHRCRuevsHOJzoZN2CcRgUZlYKPAxsATYAHzazDUNW+yxw0N03AZuBvzWzeO57iYiEoD7RRW+/s26c7lFcB9S5e7279wCPAbcPWceBaWZmwFTgDNAXYk0iIrFyqKkDgDXzx2dQLAKOp0w3BvNSfQtYD5wEXgS+4O4DQzdkZnebWY2Z1SQSibDqFRHJu0NNHZSWGCsrp0ZdyrDCDApLM8+HTL8b2AssBK4GvmVmr+sf5u7b3L3a3asrKyvHuk4Rkcgcaupg2ZzJsRy6Y1CYQdEILEmZXkxyzyHVJ4EnPKkOaADWhViTiEisvNLUyZp58T3sBOEGxW5gtZktD05Q3wlsH7LOMeCdAGY2H1gL1IdYk4hIbFzo7edIaxdr5sf3sBNAWVgbdvc+M7sX2AmUAo+4+wEzuydYvhX4KvComb1I8lDVfe7eElZNIiJx0tDSxYDDqhifyIYQgwLA3XcAO4bM25ry/CTwrjBrEBGJq8EeT6vnxXuPQldmi4hE5HBzJyUGKyqnRF1KRgoKEZGI1CU6WTp7MhVl8e3xBAoKEZHIvNLUyaqYH3YCBYWISCT6+gc40trFSgWFiIikc/zseXr7nVUxviJ7kIJCRCQCh5s7AXToSURE0tt95AwAK7RHISIi6Zzp6gFgxqQJEVeSnYJCRCQC9S1dXL98dtRl5ERBISISgcOJzoI47AQKChGRvDvT1cO57l5WxvyK7EEKChGRPKtPJHs8xflmRakUFCIieXZYQSEiIpnUJ7ooLy1h0axJUZeSEwWFiEieHWntYsnsSZSWpLtjdPwoKERE8qw+0VUwPZ5AQSEiklf9A87R1m5WzC2MHk+goBARyavGs9309A/E/mZFqRQUIiJ5VN/SBRROjydQUIiI5FV9IhkUVTr0JCIi6Rxt7WJaRRlzppRHXUrOFBQiInnU0NJF1dwpmBVG11hQUIiI5FV9oovlBXTYCRQUIiJ5c6G3n5Nt5xUUIiKS3rEz3bhTUF1jQUEhIpI3gz2etEchIiJpHWktvK6xoKAQEcmb+kQnc6eWM31i/O+TnUpBISKSJ0daugvusBMoKERE8qahtYuqOQqK1zCzW82s1szqzOz+YdbZbGZ7zeyAmf0mzHpERKLSdbGPRMfFgjs/AVAW1obNrBR4GPgToBHYbWbb3f1gyjozgW8Dt7r7MTObF1Y9IiJRqmtO3v60kIYXHxTmHsV1QJ2717t7D/AYcPuQdT4CPOHuxwDcvTnEekREInP8bDdQeD2eINygWAQcT5luDOalWgPMMrNfm9keM/tYug2Z2d1mVmNmNYlEIqRyRUTC8+qosTpH8RrpRrzyIdNlwLXAe4B3A18yszWve5H7NnevdvfqysrKsa9URCRkR1q6WDBjIpPKS6Mu5bKFdo6C5B7EkpTpxcDJNOu0uHsX0GVmzwCbgEMh1iUiknf1LYXZ4wnC3aPYDaw2s+VmVg7cCWwfss5PgLeaWZmZTQauB14KsSYRkUg0tHQV3BhPg0Lbo3D3PjO7F9gJlAKPuPsBM7snWL7V3V8ys58D+4AB4Dvuvj+smkREonC2q4e2870FebEdhHvoCXffAewYMm/rkOmHgIfCrENEJEoNrYV7Iht0ZbaISOiOtBTmYICDFBQiIiGrT3RRWmIsnT056lJGREEhIhKyI61dLJo5ifKywvzKLcyqRUQKSCHeJzuVgkJEJETuzpFWBYWIiAyjueMi3T39BXsNBSgoRERCNTjG04q5UyOuZOQyBoWZPZry/OOhVyMiUmQaXu0aW5g9niD7HsWmlOdfCLMQEZFi1NDSSUVZCQtnTIq6lBHLFhRDR3sVEZHL0BAMBlhSkm5A7cKQbQiPxWb2TZJDhg8+f5W7fz60ykREikB9oou1V0yLuoxRyRYU/y3leU2YhYiIFJve/gGOnelmyxuuiLqUUckYFO7+T/kqRESk2Bw7003fgLO8gHs8QQ7dY83s42b2RzPrCh41w92yVERELnm1a2wBX0MBWfYogkD4IvCXwB9Jnqu4BnjIzHD374VeoYhIgWpo6QRgZZHvUfwF8H53f9rd29z9nLv/CvhAsExERIZRn+hizpRyZkyeEHUpo5ItKKa7+5GhM4N508MoSESkWNQnCvf2p6myBcX5ES4TERn36ls6C3rojkHZuseuN7N9aeYbsCKEekREikLb+V5aOnuKYo8iW1BsAuYDx4fMXwacDKUiEZEiUJ9InsheUVn4exTZDj19HWh396OpD6A7WCYiImkUS9dYyB4UVe7+ukNP7l4DVIVSkYhIETic6KSsgO+TnSpbUEzMsKxwh0IUEQnZ4UQnVXOnMKG08G/7k+0T7Daz/zR0ppl9CtgTTkkiIoWvrrmTFQV8+9NU2U5mfxH4NzP7KJeCoRooB94fYl0iIgWrt3+Ao63dvPvKwh4McFC2QQGbgBvN7GbgqmD2z4Krs0VEJI2jrcnBAFcWQY8nyL5HAYC7Pw08HXItIiJFoa452TV21bziCIrCP8siIhIzh4NrKFYqKEREJJ265k6umD6RqRU5HbSJPQWFiMgYq2vuLJrDTqCgEBEZUwMDTl1zJ6vnKyhyYma3mlmtmdWZ2f0Z1nuTmfWb2QfDrEdEJGwnzp3nfG8/a+ZPi7qUMRNaUJhZKfAwsAXYAHzYzDYMs96DwM6wahERyZdXmjsAWK1DTzm5Dqhz93p37wEeA25Ps97ngMeB5hBrERHJi0NNyR5Pq7VHkZNFvHZ48sZg3qvMbBHJK7y3ZtqQmd1tZjVmVpNIJMa8UBGRsXKoqYP50yuYMamwb3+aKsygsDTzfMj0N4D73L0/04bcfZu7V7t7dWVl5VjVJyIy5l5p6iyq8xOQ45XZI9QILEmZXszrb3ZUDTxmZgBzgdvMrM/dfxxiXSIioRjs8XTndUuyr1xAwgyK3cBqM1sOnADuBD6SuoK7Lx98bmaPAv9PISEihaoYezxBiEHh7n1mdi/J3kylwCPufsDM7gmWZzwvISJSaGpPJ3s8KSgug7vvAHYMmZc2INz9E2HWIiISttqmwaAonq6xoCuzRUTGTO3pDhbNnMS0icXT4wkUFCIiY6b2dAdrryiuw06goBARGRM9fQMcTnQqKEREJL3DiU76Bpx1CgoREUnn5dPtAKxfMD3iSsaegkJEZAy8dKqD8rISVsydEnUpY05BISIyBl461c7qeVMpKy2+r9Xi+0QiInnm7hw82V6Uh51AQSEiMmqJjou0dvWwQUEhIiLpHDiZPJF95UIFhYiIpHHgZBsA6xUUIiKSzoGT7VTNmcz0Ihu6Y5CCQkRklPafbOPKhTOiLiM0CgoRkVE4193D8TPnuWqRgkJERNIYPJF91aLiPD8BCgoRkVHZ15g8kX2VDj2JiEg6L544x5LZk5g1pTzqUkKjoBARGYUXjrexcfHMqMsIlYJCRGSEWjovcuLceTYtLt7DTqCgEBEZsReOnwNgk/YoREQknRca2ygxeIP2KEREJJ3nj51lzfxpTC4vi7qUUCkoRERGYGDA2XvsHNcumxV1KaFTUIiIjMCh5g46LvZxzVIFhYiIpFFz5CwAb6qaHXEl4VNQiIiMQM2RM1ROq2DJ7ElRlxI6BYWIyAjUHD1L9bJZmFnUpYROQSEicplOt12g8ez5cXEiGxQUIiKX7bmGVgCuXz4n4kryQ0EhInKZnms4w9SKMtYvmBZ1KXkRalCY2a1mVmtmdWZ2f5rlHzWzfcHjWTPbFGY9IiJjYVd9K2+qmkVZ6fj4Wzu0T2lmpcDDwBZgA/BhM9swZLUG4O3uvhH4KrAtrHpERMZCU/sF6hNd3LBifBx2gnD3KK4D6ty93t17gMeA21NXcPdn3f1sMLkLWBxiPSIio/bbV1oAuGnV3IgryZ8wg2IRcDxlujGYN5xPAU+mW2Bmd5tZjZnVJBKJMSxRROTyPHu4hTlTyrlyYfHe+nSoMIMiXediT7ui2c0kg+K+dMvdfZu7V7t7dWVl5RiWKCKSO3fn2bpWblgxZ1xcPzEozCEPG4ElKdOLgZNDVzKzjcB3gC3u3hpiPSIio3LwVDun2y+wee34+oM1zD2K3cBqM1tuZuXAncD21BXMbCnwBHCXux8KsRYRkVF7+uVmADavnRdxJfkV2h6Fu/eZ2b3ATqAUeMTdD5jZPcHyrcCXgTnAt4PduD53rw6rJhGR0fj3l5vZtGQmldMqoi4lr0K924a77wB2DJm3NeX5p4FPh1mDiMhYaOm8yN7j5/gvt6yJupS8Gx9Xi4iIjNKvaxO4wzvWja/DTqCgEBHJya9ebmL+9Ipx1S12kIJCRCSLnr4BnjnUwjvWzRtX3WIHKShERLL4XV2Czot93LJ+ftSlREJBISKSxb89f5JZkyfw1tXj6/qJQQoKEZEMOi708tSB07x340LKy8bnV+b4/NQiIjnaeaCJi30D/NkbMw1VV9wUFCIiGfxk7wmWzp7MNUtnRl1KZBQUIiLDaGq/wO/rWvizqxeOy95OgxQUIiLD+MneEww43D6ODzuBgkJEJK2BAeeHfzjOtctmsbJyatTlREpBISKSxu/qWmho6eKuG5ZFXUrkFBQiImk88vsGKqdVsOUNV0RdSuQUFCIiQ+w/0cavaxPcdcMyKspKoy4ncgoKEZEhHtpZy8zJE/jETVVRlxILCgoRkRTP1bfym0MJPvP2lUyfOCHqcmJBQSEiEnB3/tfOWuZPr+DjN1ZFXU5sKChERAK/ermZPUfP8vl3rmbiBJ2bGKSgEBEhed3EQztrqZozmQ9VL4m6nFhRUIiIAD/4wzFePt3BX75rLRNK9dWYSq0hIuPesdZuHnjyZW5aNYc/3bgg6nJiR0EhIuNab/8AX/yX5zHggT/fOK4H/xtOWdQFiIhE6aGdtfzx2Dn+74ffyJLZk6MuJ5a0RyEi49bjexrZ9kw9d92wjD/dtDDqcmJLQSEi49IvDzZx3+P7uHHlHL703g1RlxNrCgoRGXd+vv80n/n+HjYsnM7f33XtuL0Xdq50jkJExg135x9+W8//fPJl3rhkJt/9xHVM0zAdWSkoRGRcONPVw/2P7+Opg01sueoKvvahq5lUrquvc6GgEJGidrGvn8f3nOB/P1VLx4Ve/vo96/nUW5arG+xlUFCISFFqar/AP+86yg+eO0ZrVw/Vy2bxN++/inVXTI+6tIKjoBCRotHXP0DN0bN8/7ljPPniKfrdeee6+XzixipuWjVHexEjFGpQmNmtwP8BSoHvuPsDQ5ZbsPw2oBv4hLv/McyaRKR4NLdfYP/JNg6caOeFxjb+0NBK+4U+plWU8fEbq/j4m6tYOkcX0Y1WaEFhZqXAw8CfAI3AbjPb7u4HU1bbAqwOHtcDfxf8FJFxwN3p7Xcu9vXT0zdAT/9A8mffAN09/XRe7KP9fC/nzvdypquHRMdFmtovcOLceY62dtN2vvfVba2YO4UtVy3gbWsq2by2kikVOmAyVsJsyeuAOnevBzCzx4DbgdSguB34nrs7sMvMZprZAnc/NdxGDzV1cMvXfjPsmyY3NbzMS3NbKds2stWQ2zayvT6H98i2jZwaI9s2Rt/eo/2suXyO0bZ3LlvJ/jlyeIdRtmdObRGT/yO9A05P30Au7/aqaRVlzJ8xkQUzJvLejQtYNW8qVy2awfoF05mqYAhNmC27CDieMt3I6/cW0q2zCHhNUJjZ3cDdANMXrmDt/GmZ3znLYchcjlJmO5aZbRu5HArNvo3R1ZDLSpbDVrJ9lrFpiyyfNVsNY9AYY/Nvlu31+Wjv0R+Hz8e/WVmpUVFaQsWEUspLSygvSz4qgp+Ty0uZUl7G9EkTmDFpArOnlOtmQhEJMyjS/ZoM/UMjl3Vw923ANoDq6mp/+KPXjL46ERHJSZjXrTcCqbeJWgycHME6IiISoTCDYjew2syWm1k5cCewfcg624GPWdINQFum8xMiIpJ/oR16cvc+M7sX2Emye+wj7n7AzO4Jlm8FdpDsGltHsnvsJ8OqR0RERibUbgLuvoNkGKTO25ry3IHPhlmDiIiMjsbWFRGRjBQUIiKSkYJCREQyUlCIiEhGlstwE3FiZh1AbdR1xMRcoCXqImJCbXGJ2uIStcUla909y7AW6RXi4Ci17l4ddRFxYGY1aosktcUlaotL1BaXmFnNSF+rQ08iIpKRgkJERDIqxKDYFnUBMaK2uERtcYna4hK1xSUjbouCO5ktIiL5VYh7FCIikkcKChERySi2QWFmt5pZrZnVmdn9aZabmX0zWL7PzIr2bkY5tMVHgzbYZ2bPmtmmKOrMh2xtkbLem8ys38w+mM/68imXtjCzzWa218wOmNnw9xAucDn8H5lhZj81sxeCtijKkarN7BEzazaz/cMsH9n3prvH7kFyWPLDwAqgHHgB2DBknduAJ0neJe8G4Lmo646wLW4EZgXPt4zntkhZ71ckRy7+YNR1R/h7MZPkPeqXBtPzoq47wrb4K+DB4HklcAYoj7r2ENribcA1wP5hlo/oezOuexTXAXXuXu/uPcBjwO1D1rkd+J4n7QJmmtmCfBeaB1nbwt2fdfezweQukncKLEa5/F4AfA54HGjOZ3F5lktbfAR4wt2PAbh7sbZHLm3hwDRL3lB8Ksmg6MtvmeFz92dIfrbhjOh7M65BsQg4njLdGMy73HWKweV+zk+R/IuhGGVtCzNbBLwf2Epxy+X3Yg0wy8x+bWZ7zOxjeasuv3Jpi28B60neavlF4AvuPpCf8mJlRN+bcR3Cw9LMG9qPN5d1ikHOn9PMbiYZFG8JtaLo5NIW3wDuc/f+5B+PRSuXtigDrgXeCUwC/sPMdrn7obCLy7Nc2uLdwF7gHcBK4Bdm9lt3bw+5trgZ0fdmXIOiEViSMr2Y5F8Cl7tOMcjpc5rZRuA7wBZ3b81TbfmWS1tUA48FITEXuM3M+tz9x3mpMH9y/T/S4u5dQJeZPQNsAootKHJpi08CD3jyQH2dmTUA64A/5KfE2BjR92ZcDz3tBlab2XIzKwfuBLYPWWc78LHgLP4NQJu7n8p3oXmQtS3MbCnwBHBXEf61mCprW7j7cnevcvcq4EfAXxRhSEBu/0d+ArzVzMrMbDJwPfBSnuvMh1za4hjJPSvMbD6wFqjPa5XxMKLvzVjuUbh7n5ndC+wk2aPhEXc/YGb3BMu3kuzRchtQB3ST/Iuh6OTYFl8G5gDfDv6S7vMiHDEzx7YYF3JpC3d/ycx+DuwDBoDvuHvabpOFLMffi68Cj5rZiyQPv9zn7kU3/LiZ/RDYDMw1s0bgK8AEGN33pobwEBGRjOJ66ElERGJCQSEiIhkpKEREJCMFhYiIZKSgEBGRjBQUUjSC0WL3pjzuD+a/NRgxdK+ZTTKzh4Lph0bwHn81ZPrZsapfJK7UPVaKhpl1uvvUNPO3khwl87vBdDtQ6e4Xx+o9RIqZ9iikqJnZp4EPAV82s++b2XZgCvCcmd1hZpVm9riZ7Q4eNwWvm2pm3zWzF4Nx+z9gZg8Ak4I9k+8H63UGP//FzG5Led9Hg9eUBnswu4Pt/Odg+eZgsL4fmdnLQW0WLLvWzH4TDOS3c3B0TzP7vJkdDLbzWDDv7Sl7UM+b2bS8Na6MH1GPn66HHmP1APpJDvw2+LgjmP8oKfelADpTnv8AeEvwfCnwUvD8QeAbKevNGvra1GmSI9b+U/C8nOQInZOAu4G/DuZXADXAcpJXz7aRHGunBPgPkoM5TgCeJbnHA3AHySuNITkmT0XwfGbw86fATcHzqUBZ1P8OehTfI5ZDeIiM0Hl3v/oyX3MLsCFlpNnpwV/lt5AcMwgAv3S/j+E8CXzTzCqAW4Fn3P28mb0L2GiX7rQ3A1gN9AB/cPdGADPbC1QB54CrSI5uCskhKQbH4tkHfN/Mfgz8OJj3e+BrwR7OE4PbExlLCgoZ70qAN7v7+dSZwWGgnE/gufsFM/s1yeGs7wB+OLgp4HPuvnPI9jcDqedI+kn+fzTggLu/Oc3bvIfkHczeB3zJzK509wfM7Gckx+/ZZWa3uPvLudYtkgudo5Dx7ing3sEJM7t6mPmzgqe9ZjZhmG09RnKQtbeSHKCO4OdnBl9jZmvMbEqGemqBSjN7c7D+BDO70sxKgCXu/jTw30ne5nSqma109xfd/UGSh7XW5faxRXKnoJBiMsle2z32gRxe83mgOjhBfBC4J5j/NyTvDrffzF4Abg7mbwP2DZ7MHuIpkn/x/9KTt+SE5D1CDgJ/tOQN7/+eDHvywes+CDwYvO9ekvdELwX+ORj99Hng6+5+DvhiSo3nKd67G0qE1D1WREQy0h6FiIhkpKAQEZGMFBQiIpKRgkJERDJSUIiISEYKChERyUhBISIiGf1/AQsGs3QMPXIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "effectiveness.make_cdf().plot()\n",
    "plt.xlim([0, 1])\n",
    "\n",
    "plt.xlabel('Effectiveness')\n",
    "plt.ylabel('CDF');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The posterior mean is about 89%."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.8949353341973734"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "effectiveness.mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And the 95% credible interval is between 81% and 95%.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.8099631 , 0.95352564])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "effectiveness.credible_interval(0.95)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If my guesses about the data are close enough, and the modeling decisions are good enough, it is unlikely that the effectiveness of the vaccine is less than 80%."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}