workshop01.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Bayesian Statistics Made Simple\n",
    "===\n",
    "\n",
    "Code and exercises from my workshop on Bayesian statistics in Python.\n",
    "\n",
    "Copyright 2016 Allen Downey\n",
    "\n",
    "MIT License: https://opensource.org/licenses/MIT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from __future__ import print_function, division\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "import math\n",
    "import numpy as np\n",
    "from scipy.special import gamma\n",
    "\n",
    "from thinkbayes2 import Pmf, Suite\n",
    "import thinkplot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Working with Pmfs\n",
    "---\n",
    "Create a Pmf object to represent a six-sided die."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "d6 = Pmf()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A Pmf is a map from possible outcomes to their probabilities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for x in [1,2,3,4,5,6]:\n",
    "    d6[x] = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initially the probabilities don't add up to 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "d6.Print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`Normalize` adds up the probabilities and divides through.  The return value is the total probability before normalizing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "d6.Normalize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the Pmf is normalized."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "d6.Print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can compute its mean (which only works if it's normalized)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "d6.Mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`Random` chooses a random value from the Pmf."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "d6.Random()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`thinkplot` provides methods for plotting Pmfs in a few different styles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "thinkplot.Hist(d6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 1:**  The Pmf object provides `__add__`, so you can use the `+` operator to compute the Pmf of the sum of two dice.\n",
    "\n",
    "Compute and plot the Pmf of the sum of two 6-sided dice."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 2:** Suppose I roll two dice and tell you the result is greater than 3.\n",
    "\n",
    "Plot the Pmf of the remaining possible outcomes and compute its mean."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The cookie problem\n",
    "---\n",
    "Create a Pmf with two equally likely hypotheses.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cookie = Pmf(['Bowl 1', 'Bowl 2'])\n",
    "cookie.Print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Update each hypothesis with the likelihood of the data (a vanilla cookie)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cookie['Bowl 1'] *= 0.75\n",
    "cookie['Bowl 2'] *= 0.5\n",
    "cookie.Normalize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Print the posterior probabilities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cookie.Print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 3:** Suppose we put the first cookie back, stir, choose again from the same bowl, and get a chocolate cookie.\n",
    "\n",
    "Hint: The posterior (after the first cookie) becomes the prior (before the second cookie)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 4:** Instead of doing two updates, what if we collapse the two pieces of data into one update?\n",
    "\n",
    "Re-initialize `Pmf` with two equally likely hypotheses and perform one update based on two pieces of data, a vanilla cookie and a chocolate cookie.\n",
    "\n",
    "The result should be the same regardless of how many updates you do (or the order of updates)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The dice problem\n",
    "---\n",
    "Create a Suite to represent dice with different numbers of sides."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "suite = Suite([4, 6, 8, 12])\n",
    "suite.Print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 5:** We'll solve this problem two ways.  First we'll do it \"by hand\", as we did with the cookie problem; that is, we'll multiply each hypothesis by the likelihood of the data, and the renormalize.\n",
    "\n",
    "In the space below, update `suite` based on the likelihood of the data (rolling a 6), then normalize and print the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 6:**  Now let's do the same calculation using `Suite.Update`.\n",
    "\n",
    "Write a definition for a new class called `Dice` that extends `Suite`.  Then define a method called `Likelihood` that takes `data` and `hypo` and returns the probability of the data (the outcome of rolling the die) for a given hypothesis (number of sides on the die).\n",
    "\n",
    "Hint: What should you do if the outcome exceeds the hypothetical number of sides on the die?\n",
    "\n",
    "Here's an outline to get you started:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class Dice(Suite):\n",
    "        # hypo is the number of sides on the die\n",
    "        # data is the outcome\n",
    "        def Likelihood(self, data, hypo):\n",
    "            return 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can create a `Dice` object and update it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dice = Dice([4, 6, 8, 12])\n",
    "dice.Update(6)\n",
    "dice.Print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we get more data, we can perform more updates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for roll in [8, 7, 7, 5, 4]:\n",
    "    dice.Update(roll)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here are the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dice.Print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The German tank problem\n",
    "---\n",
    "The German tank problem is actually identical to the dice problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "class Tank(Suite):\n",
    "    # hypo is the number of tanks\n",
    "    # data is an observed serial number\n",
    "    def Likelihood(self, data, hypo):\n",
    "        if data > hypo:\n",
    "            return 0\n",
    "        else:\n",
    "            return 1 / hypo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here are the posterior probabilities after seeing Tank #37."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "tank = Tank(range(100))\n",
    "tank.Update(37)\n",
    "thinkplot.Pdf(tank)\n",
    "tank.Mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 7:**  Suppose we see another tank with serial number 17.  What effect does this have on the posterior probabilities?\n",
    "\n",
    "Update the suite again with the new data and plot the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Euro problem\n",
    "---\n",
    "\n",
    "**Exercise 8:**  Write a class definition for `Euro`, which extends `Suite` and defines a likelihood function that computes the probability of the data (heads or tails) for a given value of `x` (the probability of heads).\n",
    "\n",
    "Note that `hypo` is in the range 0 to 100.  Here's an outline to get you started."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class Euro(Suite):\n",
    "    \n",
    "    def Likelihood(self, data, hypo):\n",
    "        \"\"\" \n",
    "        hypo is the prob of heads (0-100)\n",
    "        data is a string, either 'H' or 'T'\n",
    "        \"\"\"\n",
    "        return 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll start with a uniform distribution from 0 to 100."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "euro = Euro(range(101))\n",
    "thinkplot.Pdf(euro)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can update with a single heads:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "euro.Update('H')\n",
    "thinkplot.Pdf(euro)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another heads:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "euro.Update('H')\n",
    "thinkplot.Pdf(euro)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And a tails:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "euro.Update('T')\n",
    "thinkplot.Pdf(euro)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Starting over, here's what it looks like after 7 heads and 3 tails."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "euro = Euro(range(101))\n",
    "\n",
    "for outcome in 'HHHHHHHTTT':\n",
    "    euro.Update(outcome)\n",
    "\n",
    "thinkplot.Pdf(euro)\n",
    "euro.MaximumLikelihood()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The maximum posterior probability is 70%, which is the observed proportion.\n",
    "\n",
    "Here are the posterior probabilities after 140 heads and 110 tails."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "euro = Euro(range(101))\n",
    "\n",
    "evidence = 'H' * 140 + 'T' * 110\n",
    "for outcome in evidence:\n",
    "    euro.Update(outcome)\n",
    "    \n",
    "thinkplot.Pdf(euro)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The posterior mean s about 56%"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "euro.Mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So is the value with maximum posterior probability."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "euro.MaximumLikelihood()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The posterior credible interval has a 90% chance of containing the true value (provided that the prior distribution truly represents our background knowledge)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "euro.CredibleInterval(90)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The World Cup Problem\n",
    "---\n",
    "We'll use λ to represent the hypothetical goal-scoring rate in goals per game.\n",
    "\n",
    "To compute prior probabilities for values of λ, we'll use a Gamma distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def EvalGammaPdf(lam, a):\n",
    "    \"\"\"Computes the Gamma PDF.\n",
    "\n",
    "    lam: where to evaluate the PDF\n",
    "    a: parameter of the gamma distribution\n",
    "\n",
    "    returns: float probability\n",
    "    \"\"\"\n",
    "    return lam**(a-1) * math.exp(-lam) / gamma(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`MakeGammaPmf` takes a list of hypothetical rates and returns a Pmf object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def MakeGammaPmf(lams, a):\n",
    "    \"\"\"Makes a PMF discrete approx to a Gamma distribution.\n",
    "\n",
    "    lam: parameter lambda in events per unit time\n",
    "    xs: upper bound of the Pmf\n",
    "\n",
    "    returns: normalized Pmf\n",
    "    \"\"\"\n",
    "    pmf = Pmf()\n",
    "    for lam in lams:\n",
    "        pmf[lam] = EvalGammaPdf(lam, a)\n",
    "        \n",
    "    pmf.Normalize()\n",
    "    return pmf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what the prior looks like.  I chose the parameter to yield a mean close to 1.3, which is the average number of goals per team per game in World Cup play."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "xs = np.linspace(0, 8, 101)\n",
    "pmf = MakeGammaPmf(xs, 1.3)\n",
    "thinkplot.Pdf(pmf)\n",
    "thinkplot.Config(xlabel='Goals per game')\n",
    "pmf.Mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 9:**  Write a class called `Soccer` that extends `Suite` and defines `Likelihood`, which should compute the probability of the data (the time between goals in minutes) for a hypothetical goal-scoring rate, `lam`, in goals per game.\n",
    "\n",
    "Hint: For a given value of `lam`, the time between goals is distributed exponentially.\n",
    "\n",
    "Here's an outline to get you started:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class Soccer(Suite):\n",
    "    \"\"\"Represents hypotheses about goal-scoring rates.\"\"\"\n",
    "\n",
    "    def Likelihood(self, data, hypo):\n",
    "        \"\"\"Computes the likelihood of the data under the hypothesis.\n",
    "\n",
    "        hypo: scoring rate in goals per game\n",
    "        data: interarrival time in minutes\n",
    "        \"\"\"\n",
    "        return 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can create a `Soccer` object and initialize it with the prior Pmf:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "soccer = Soccer(pmf)\n",
    "thinkplot.Pdf(soccer)\n",
    "thinkplot.Config(xlabel='Goals per game')\n",
    "soccer.Mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's the update after first goal at 11 minutes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "thinkplot.Pdf(soccer, color='0.7')\n",
    "soccer.Update(11)\n",
    "thinkplot.Pdf(soccer)\n",
    "thinkplot.Config(xlabel='Goals per game')\n",
    "soccer.Mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's the update after the second goal at 23 minutes (the time between first and second goals is 12 minutes).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "thinkplot.Pdf(soccer, color='0.7')\n",
    "soccer.Update(12)\n",
    "thinkplot.Pdf(soccer)\n",
    "thinkplot.Config(xlabel='Goals per game')\n",
    "soccer.Mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This distribution represents our belief about `lam` after two goals.\n",
    "\n",
    "### Predictive distributions\n",
    "\n",
    "Now to predict the number of goals in the remaining 67 minutes.\n",
    "\n",
    "The distribution of additional goals (after the first two) is Poisson, with a parameter that is the product of `lam` and the remaining time (in games)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def EvalPoissonPmf(k, lam):\n",
    "    \"\"\"Computes the Poisson PMF.\n",
    "\n",
    "    k: number of events\n",
    "    lam: parameter lambda in events per unit time\n",
    "\n",
    "    returns: float probability\n",
    "    \"\"\"\n",
    "    return lam ** k * math.exp(-lam) / gamma(k+1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`MakePoissonPmf` makes a Pmf that represents a Poisson distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def MakePoissonPmf(lam, high, step=1):\n",
    "    \"\"\"Makes a PMF discrete approx to a Poisson distribution.\n",
    "\n",
    "    lam: parameter lambda in events per unit time\n",
    "    high: upper bound of the Pmf\n",
    "\n",
    "    returns: normalized Pmf\n",
    "    \"\"\"\n",
    "    pmf = Pmf()\n",
    "    for k in range(0, high + 1, step):\n",
    "        p = EvalPoissonPmf(k, lam)\n",
    "        pmf[k] = p\n",
    "    pmf.Normalize()\n",
    "    return pmf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we assume that `lam` is the mean of the posterior, we can generate a predictive distribution for the number of goals in the remainder of the game."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "lam = soccer.Mean()\n",
    "rem_time = 90 - 23\n",
    "lt = lam * rem_time / 90\n",
    "pred = MakePoissonPmf(lt, 10)\n",
    "thinkplot.Hist(pred)\n",
    "thinkplot.Config(title='Option 1', \n",
    "                 xlabel='Expected goals',\n",
    "                 xlim=[-0.5, 10.5])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The predictive mean is about 2 goals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pred.Mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And the chance of scoring 5 more goals is still small."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pred.ProbGreater(4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But that answer is only approximate because it does not take into account our uncertainty about `lam`.\n",
    "\n",
    "The correct method is to compute a weighted mixture of Poisson distributions, one for each possible value of `lam`.\n",
    "\n",
    "The following figure shows the different predictive distributions for the different values of `lam`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for lam, prob in soccer.Items():\n",
    "    lt = lam * rem_time / 90\n",
    "    pred = MakePoissonPmf(lt, 14)\n",
    "    thinkplot.Pdf(pred, color='gray', alpha=0.3, linewidth=0.5)\n",
    "\n",
    "thinkplot.Config(xlabel='Expected goals')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can compute the mixture of these distributions by making a Meta-Pmf that maps from each Poisson Pmf to its probability."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "metapmf = Pmf()\n",
    "\n",
    "for lam, prob in soccer.Items():\n",
    "    lt = lam * rem_time / 90\n",
    "    pred = MakePoissonPmf(lt, 15)\n",
    "    metapmf[pred] = prob"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`MakeMixture` takes a Meta-Pmf (a Pmf that contains Pmfs) and returns a single Pmf that represents the weighted mixture of distributions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def MakeMixture(metapmf, label='mix'):\n",
    "    \"\"\"Make a mixture distribution.\n",
    "\n",
    "    Args:\n",
    "      metapmf: Pmf that maps from Pmfs to probs.\n",
    "      label: string label for the new Pmf.\n",
    "\n",
    "    Returns: Pmf object.\n",
    "    \"\"\"\n",
    "    mix = Pmf(label=label)\n",
    "    for pmf, p1 in metapmf.Items():\n",
    "        for x, p2 in pmf.Items():\n",
    "            mix[x] += p1 * p2\n",
    "    return mix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's the result for the World Cup problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "mix = MakeMixture(metapmf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here's what the mixture looks like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "thinkplot.Hist(mix)\n",
    "thinkplot.Config(title='Option 2', \n",
    "                 xlabel='Expected goals',\n",
    "                 xlim=[-0.5, 10.5])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 10:** Compute the predictive mean and the probability of scoring 5 or more additional goals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}