09_predict.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Prediction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "This notebook is part of [Bite Size Bayes](https://allendowney.github.io/BiteSizeBayes/), an introduction to probability and Bayesian statistics using Python.\n",
    "\n",
    "Copyright 2020 Allen B. Downey\n",
    "\n",
    "License: [Attribution-NonCommercial-ShareAlike 4.0 International (CC BY-NC-SA 4.0)](https://creativecommons.org/licenses/by-nc-sa/4.0/)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "The following cell downloads `utils.py`, which contains some utility function we'll need."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from os.path import basename, exists\n",
    "\n",
    "def download(url):\n",
    "    filename = basename(url)\n",
    "    if not exists(filename):\n",
    "        from urllib.request import urlretrieve\n",
    "        local, _ = urlretrieve(url, filename)\n",
    "        print('Downloaded ' + local)\n",
    "\n",
    "download('https://github.com/AllenDowney/BiteSizeBayes/raw/master/utils.py')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "If everything we need is installed, the following cell should run with no error messages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Review\n",
    "\n",
    "[In the previous notebook](https://colab.research.google.com/github/AllenDowney/BiteSizeBayes/blob/master/08_soccer.ipynb) I introduced the World Cup problem:\n",
    "\n",
    ">In the 2018 FIFA World Cup final, France defeated Croatia 4 goals to 2.  Based on this outcome:\n",
    ">\n",
    ">1. How confident should we be that France is the better team?\n",
    ">\n",
    ">2. If the same teams played again, what is the chance Croatia would win?\n",
    "\n",
    "I started with the assumption that for any team against any other team there is some unknown goal-scoring rate, λ.\n",
    "\n",
    "And I showed that if we know λ, we can compute the probability of scoring $k$ goals in a game:\n",
    "\n",
    "$f(k; λ) = λ^k \\exp(-λ) ~/~ k!$\n",
    "\n",
    "This function is the [Poisson distribution](https://en.wikipedia.org/wiki/Poisson_distribution), and we can use SciPy to compute it.\n",
    "\n",
    "For example, if we know that λ is 1.4, we can compute the distribution of $k$ like this:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import poisson\n",
    "\n",
    "λ = 1.4\n",
    "xs = np.arange(11)\n",
    "ys = poisson.pmf(xs, λ)\n",
    "pmf_goals = pd.Series(ys, index=xs)\n",
    "\n",
    "pmf_goals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In theory, the values of $k$ go to infinity, but I cut it off at 10 because higher values have very low probability.\n",
    "\n",
    "Here's what the distribution of $k$ looks like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "pmf_goals.plot.bar(color='C0', alpha=0.5, label='Poisson distribution')\n",
    "\n",
    "plt.xlabel('Number of goals')\n",
    "plt.ylabel('Probability')\n",
    "plt.title('Distribution of goals scored')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because the values of $k$ are discrete integers, I use a bar chart to plot the distribution.\n",
    "\n",
    "Next I used a gamma distributon to represent the prior distribution of λ.  I set the parameter of the gamma distribution, α, to 1.4, which is the average number of goals each team scores in World Cup play.\n",
    "\n",
    "I broke the range of values for λ into 101 discrete possibilities and used SciPy to compute the prior probabilities:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import gamma\n",
    "\n",
    "α = 1.4\n",
    "qs = np.linspace(0, 10, 101)\n",
    "ps = gamma.pdf(qs, α)\n",
    "\n",
    "prior = pd.Series(ps, index=qs)\n",
    "prior /= prior.sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what the prior distribution looks like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "prior.plot(label='prior', color='gray')\n",
    "\n",
    "plt.xlabel('Goal scoring rate')\n",
    "plt.ylabel('Probability')\n",
    "plt.title('Prior distribution')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because the values of $λ$ are continuous, I use a line to plot the distribution.\n",
    "\n",
    "Next we did a Bayesian update, using the Poisson distribution to compute the likelihood of the data, $k$, for each possible goal scoring rate, $λ$\n",
    "\n",
    "$ f(k; λ) = λ^k \\exp(-λ) ~/~ k! $\n",
    "\n",
    "Since France scored 4 goals, the likelihood of the data is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "λs = prior.index\n",
    "k = 4\n",
    "likelihood = poisson.pmf(k, λs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can use the following function to do the update:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bayes_update(pmf, likelihood):\n",
    "    \"\"\"Do a Bayesian update.\n",
    "    \n",
    "    pmf: Series that represents the prior\n",
    "    likelihood: sequence of likelihoods\n",
    "    \n",
    "    returns: float probability of the data\n",
    "    \"\"\"\n",
    "    pmf *= likelihood\n",
    "    prob_data = pmf.sum()\n",
    "    pmf /= prob_data\n",
    "    return prob_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "france = prior.copy()\n",
    "bayes_update(france, likelihood)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can do the same thing for Croatia, given that they scored 2 goals:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "k = 2\n",
    "λs = prior.index\n",
    "likelihood = poisson.pmf(k, λs)\n",
    "\n",
    "croatia = prior.copy()\n",
    "bayes_update(croatia, likelihood)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what the results look like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "prior.plot(label='prior', color='gray')\n",
    "france.plot(label='France', color='C0')\n",
    "croatia.plot(label='Croatia', color='C3')\n",
    "\n",
    "plt.xlabel('Goal scoring rate')\n",
    "plt.ylabel('Probability')\n",
    "plt.title('Posterior distributions for France and Croatia')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the previous notebook we used the following function to compute the mean of a PMF."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pmf_mean(pmf):\n",
    "    \"\"\"Compute the mean of a PMF.\n",
    "    \n",
    "    pmf: Series representing a PMF\n",
    "    \n",
    "    return: float\n",
    "    \"\"\"\n",
    "    return np.sum(pmf.index * pmf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use it to compute the posterior means for France and Croatia."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "pmf_mean(france), pmf_mean(croatia)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Based on the outcome of the game, we have some reason to think France is the better team.\n",
    "\n",
    "But if we look at the posterior distribution of λ for France and Croatia, there is still a lot of overlap, which means we are still not certain which is the better team.\n",
    "\n",
    "In the next section, we'll use the posterior distributions to compute the \"probability of superiority\"."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Probability of superiority\n",
    "\n",
    "Now that we have a posterior distribution for each team, we can answer the first question: How confident should we be that France is the better team?\n",
    "\n",
    "In the model, \"better\" means having a higher goal-scoring rate against the opponent.  We can use the posterior distributions to compute the probability that a random value drawn from France's distribution exceeds a value drawn from Croatia's.\n",
    "\n",
    "The following function takes a Series that represents a PMF and draws a sample from it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pmf_choice(pmf, n):\n",
    "    \"\"\"Draw a random sample from a PMF.\n",
    "    \n",
    "    pmf: Series representing a PMF\n",
    "    n: number of values to draw\n",
    "    \n",
    "    returns: NumPy array\n",
    "    \"\"\"\n",
    "    return np.random.choice(pmf.index,\n",
    "                            size=n, \n",
    "                            replace=True, \n",
    "                            p=pmf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`pmf_choice` uses `np.random.choice`, which chooses `n` values from the PMF with replacement, so the same value can appear more than once.  It uses the probabilities from the PMF as weights, so the number of times each value appears is proportional to its probability.\n",
    "\n",
    "Here's a sample from the posterior distribution for France."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_france = pmf_choice(france, 1000)\n",
    "sample_france.mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here's a sample for Croatia."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_croatia = pmf_choice(croatia, 1000)\n",
    "sample_croatia.mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To estimate the probability of superiority, we can count the number of times the value from France's distribution exceeds the value from Croatia's distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.mean(sample_france > sample_croatia)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On the basis of one game, we have moderate confidence that France is actually the better team."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Predicting the rematch\n",
    "\n",
    "Now we can take on the second question: If the same teams played again, what is the chance Croatia would win?\n",
    "\n",
    "To answer this question, we'll generate a sample from the \"posterior predictive distribution\", which is the number of goals we expect a team to score.\n",
    "\n",
    "If we knew the goal scoring rate, λ, the distribution of goals would be a Poisson distributon with parameter λ.\n",
    "\n",
    "Since we don't know λ, we can use the sample we generated in the previous section to generate a sample of goals, like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "goals_france = np.random.poisson(sample_france)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`np.random.poisson` uses each element of `sample_france` to generate a random value; that is, each value in the result is based on a different value of λ.\n",
    "\n",
    "To see what the resulting distribution looks like, we'll use this function from a previous notebook:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pmf_from_seq(seq):\n",
    "    \"\"\"Make a PMF from a sequence of values.\n",
    "    \n",
    "    seq: sequence\n",
    "    \n",
    "    returns: Series representing a PMF\n",
    "    \"\"\"\n",
    "    pmf = pd.Series(seq).value_counts(sort=False).sort_index()\n",
    "    pmf /= pmf.sum()\n",
    "    return pmf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "pmf_france = pmf_from_seq(goals_france)\n",
    "pmf_france.plot.bar(color='C0', label='France')\n",
    "\n",
    "plt.xlabel('Goals scored')\n",
    "plt.ylabel('Probability')\n",
    "plt.title('Predictive distribution')\n",
    "plt.legend()\n",
    "\n",
    "goals_france.mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This distribution represents two sources of uncertainty: we don't know the actual value of λ, and even if we did, we would not know the number of goals in the next game."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Generate and plot the predictive distribution for Croatia."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In a sense, these distributions represent the outcomes of 1000 simulated games."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Compute the fraction of simulated rematches Croatia would win, how many France would win, and how many would end in a tie."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assuming that Croatia wins half of the ties, their chance of winning the rematch is about 33%."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "In this notebook, we finished off the World Cup problem:\n",
    "\n",
    "* We used posterior distributions to generate samples of goal-scoring rates.\n",
    "\n",
    "* We compared samples to compute a \"probability of superiority\".\n",
    "\n",
    "* We used samples and `np.random.poisson` to generate samples of goals score and to estimate their distributions.\n",
    "\n",
    "* We used those distributions to compute the probabilities of winning, losing, and tying in a rematch.\n",
    "\n",
    "The goal distributions we computed are called \"[posterior predictive distributions](https://en.wikipedia.org/wiki/Posterior_predictive_distribution)\" because they use posterior distribution to make predictions.\n",
    "\n",
    "[In the next notebook](https://colab.research.google.com/github/AllenDowney/BiteSizeBayes/blob/master/10_joint.ipynb) we'll take a break from Bayesian methods and learn about an important idea in probability: joint distributions.\n",
    "\n",
    "But first, here's an exercise where you can practice what you learned in this notebook."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise\n",
    "\n",
    "**Exercise:** Here's a variation on the World Cup Problem:\n",
    "\n",
    ">In the 2014 FIFA World Cup, Germany played Brazil in a semifinal match. Germany scored after 11 minutes and again at the 23 minute mark. At that point in the match, how many goals would you expect Germany to score after 90 minutes? What was the probability that they would score 5 more goals (as, in fact, they did)?\n",
    "\n",
    "In this version, notice that the data is not the number of goals in a fixed period of time, but the time between goals.\n",
    "\n",
    "To compute the likelihood of data like this, we can take advantage of the theory of [Poisson processes](https://en.wikipedia.org/wiki/Poisson_point_process).  In our model of a soccer game, we assume that each team has a goal-scoring rate, λ, in goals per game.  And we assume that λ is constant, so the chance of scoring a goal in the same at any moment of the game."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Under these assumptions, the distribution of goals follows a Poisson distribution, as we've already seen.  Also, the time between goals follows an [exponential distribution](https://en.wikipedia.org/wiki/Exponential_distribution).\n",
    "\n",
    "If the goal-scoring rate is λ, the probability of seeing an interval between goals of $t$ is proportional to the PDF of the exponential distribution:\n",
    "\n",
    "$f(t; λ) = λ~\\exp(-λ t)$\n",
    "\n",
    "Because $t$ is a continuous quantity, the value of this expression is not really a probability; technically it is a [probability density](https://en.wikipedia.org/wiki/Probability_density_function).  However, it is proportional to the probability of the data, so we can use it as a likelihood in a Bayesian update.\n",
    "\n",
    "To see what the exponential distribution looks like, let's assume again that λ is 1.4; we can compute the distribution of $t$ like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "def expo_pdf(t, λ):\n",
    "    \"\"\"Compute the PDF of the exponential distribution.\n",
    "    \n",
    "    t: time\n",
    "    λ: rate\n",
    "    \n",
    "    returns: probability density\n",
    "    \"\"\"\n",
    "    return λ * np.exp(-λ * t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "λ = 1.4\n",
    "ts = np.linspace(0, 4, 101)\n",
    "ys = expo_pdf(ts, λ)\n",
    "pmf_time = pd.Series(ys, index=ts)\n",
    "pmf_time /= pmf_time.sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "pmf_time.plot(label='λ = 1.4')\n",
    "\n",
    "plt.xlabel('Time between goals (games)')\n",
    "plt.ylabel('Probability')\n",
    "plt.title('Distribution of time between goals')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is counterintuitive, but true, that the most likely time to score a goal is immediately.  After that, the probability of each possible interval is a little lower.\n",
    "\n",
    "With a goal-scoring rate of 1.4, it is possible that a team will take more than one game to score a goal, but it is unlikely that they will take more than two games."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we're ready to solve the problem.  Here are the steps I recommend:\n",
    "\n",
    "1. Starting with the same gamma prior we used in the previous problem, compute the likelihood of scoring a goal after 11 minutes for each possible value of λ.  Don't forget to convert all times into units of games.\n",
    "\n",
    "2. Use `bayes_update` to compute the posterior distribution of λ for Germany after the first goal.  If the total probability of the data is greater than 1, don't panic; because the likelihoods are not probabilities, the return value from `bayes_update` is not a probability either.  But the posterior distribution is still valid.\n",
    "\n",
    "3. Compute the likelihood of scoring another goal after 12 more minutes and do another update.  Plot the prior, posterior after one goal, and posterior after two goals.\n",
    "\n",
    "4. Use `pmf_choice` to generate a sample of 1000 values from the posterior distribution of goal scoring rate.\n",
    "\n",
    "5. Use `np.random.poisson` to compute a sample of goals Germany might score during the remaining time in the game, `90-23` minutes.  Note: you will have to think about how to generate predicted goals for a fraction of a game.\n",
    "\n",
    "6. Compute and plot the PMF of possible goals scored and select from it the probability of scoring 5 more goals.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "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.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}