{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Think Bayes\n",
    "\n",
    "This notebook presents example code and exercise solutions for Think Bayes.\n",
    "\n",
    "Copyright 2018 Allen B. Downey\n",
    "\n",
    "MIT License: https://opensource.org/licenses/MIT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Configure Jupyter so figures appear in the notebook\n",
    "%matplotlib inline\n",
    "\n",
    "# Configure Jupyter to display the assigned value after an assignment\n",
    "%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'\n",
    "\n",
    "# import classes from thinkbayes2\n",
    "from thinkbayes2 import Hist, Pmf, Suite, Beta\n",
    "import thinkplot\n",
    "\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparing distributions\n",
    "\n",
    "Let's get back to the Kim Rhode problem from Chapter 4:\n",
    "\n",
    "> At the 2016 Summer Olympics in the Women's Skeet event, Kim Rhode faced Wei Meng in the bronze medal match. They each hit 15 of 25 targets, sending the match into sudden death. In the first round, both hit 1 of 2 targets. In the next two rounds, they each hit 2 targets. Finally, in the fourth round, Rhode hit 2 and Wei hit 1, so Rhode won the bronze medal, making her the first Summer Olympian to win an individual medal at six consecutive summer games.\n",
    "\n",
    ">But after all that shooting, what is the probability that Rhode is actually a better shooter than Wei? If the same match were held again, what is the probability that Rhode would win?\n",
    "\n",
    "I'll start with a uniform distribution for `x`, the probability of hitting a target, but we should check whether the results are sensitive to that choice.\n",
    "\n",
    "First I create a Beta distribution for each of the competitors, and update it with the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "rhode = Beta(1, 1, label='Rhode')\n",
    "rhode.Update((22, 11))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "wei = Beta(1, 1, label='Wei')\n",
    "wei.Update((21, 12))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Based on the data, the distribution for Rhode is slightly farther right than the distribution for Wei, but there is a lot of overlap."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "thinkplot.Pdf(rhode.MakePmf())\n",
    "thinkplot.Pdf(wei.MakePmf())\n",
    "thinkplot.Config(xlabel='x', ylabel='Probability')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To compute the probability that Rhode actually has a higher value of `p`, there are two options:\n",
    "\n",
    "1. Sampling: we could draw random samples from the posterior distributions and compare them.\n",
    "\n",
    "2. Enumeration: we could enumerate all possible pairs of values and add up the \"probability of superiority\".\n",
    "\n",
    "I'll start with sampling.  The Beta object provides a method that draws a random value from a Beta distribution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "iters = 1000\n",
    "count = 0\n",
    "for _ in range(iters):\n",
    "    x1 = rhode.Random()\n",
    "    x2 = wei.Random()\n",
    "    if x1 > x2:\n",
    "        count += 1\n",
    "\n",
    "count / iters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`Beta` also provides `Sample`, which returns a NumPy array, so we an perform the comparisons using array operations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "rhode_sample = rhode.Sample(iters)\n",
    "wei_sample = wei.Sample(iters)\n",
    "np.mean(rhode_sample > wei_sample)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The other option is to make `Pmf` objects that approximate the Beta distributions, and enumerate pairs of values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def ProbGreater(pmf1, pmf2):\n",
    "    total = 0\n",
    "    for x1, prob1 in pmf1.Items():\n",
    "        for x2, prob2 in pmf2.Items():\n",
    "            if x1 > x2:\n",
    "                total += prob1 * prob2\n",
    "    return total"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pmf1 = rhode.MakePmf(1001)\n",
    "pmf2 = wei.MakePmf(1001)\n",
    "ProbGreater(pmf1, pmf2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pmf1.ProbGreater(pmf2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pmf1.ProbLess(pmf2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Run this analysis again with a different prior and see how much effect it has on the results."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulation\n",
    "\n",
    "To make predictions about a rematch, we have two options again:\n",
    "\n",
    "1. Sampling.  For each simulated match, we draw a random value of `x` for each contestant, then simulate 25 shots and count hits.\n",
    "\n",
    "2. Computing a mixture.  If we knew `x` exactly, the distribution of hits, `k`, would be binomial.  Since we don't know `x`, the distribution of `k` is a mixture of binomials with different values of `x`.\n",
    "\n",
    "I'll do it by sampling first."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import random\n",
    "\n",
    "def flip(p):\n",
    "    return random.random() < p"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`flip` returns True with probability `p` and False with probability `1-p`\n",
    "\n",
    "Now we can simulate 1000 rematches and count wins and losses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "iters = 1000\n",
    "wins = 0\n",
    "losses = 0\n",
    "\n",
    "for _ in range(iters):\n",
    "    x1 = rhode.Random()\n",
    "    x2 = wei.Random()\n",
    "    \n",
    "    count1 = count2 = 0\n",
    "    for _ in range(25):\n",
    "        if flip(x1):\n",
    "            count1 += 1\n",
    "        if flip(x2):\n",
    "            count2 += 1\n",
    "            \n",
    "    if count1 > count2:\n",
    "        wins += 1\n",
    "    if count1 < count2:\n",
    "        losses += 1\n",
    "        \n",
    "wins/iters, losses/iters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or, realizing that the distribution of `k` is binomial, we can simplify the code using NumPy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "rhode_rematch = np.random.binomial(25, rhode_sample)\n",
    "thinkplot.Hist(Pmf(rhode_rematch))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "wei_rematch = np.random.binomial(25, wei_sample)\n",
    "np.mean(rhode_rematch > wei_rematch)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "np.mean(rhode_rematch < wei_rematch)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, we can make a mixture that represents the distribution of `k`, taking into account our uncertainty about `x`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from thinkbayes2 import MakeBinomialPmf\n",
    "\n",
    "def MakeBinomialMix(pmf, label=''):\n",
    "    mix = Pmf(label=label)\n",
    "    for x, prob in pmf.Items():\n",
    "        binom = MakeBinomialPmf(n=25, p=x)\n",
    "        for k, p in binom.Items():\n",
    "            mix[k] += prob * p\n",
    "    return mix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "rhode_rematch = MakeBinomialMix(rhode.MakePmf(), label='Rhode')\n",
    "wei_rematch = MakeBinomialMix(wei.MakePmf(), label='Wei')\n",
    "thinkplot.Pdf(rhode_rematch)\n",
    "thinkplot.Pdf(wei_rematch)\n",
    "thinkplot.Config(xlabel='hits')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "rhode_rematch.ProbGreater(wei_rematch), rhode_rematch.ProbLess(wei_rematch)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, we could use MakeMixture:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from thinkbayes2 import MakeMixture\n",
    "\n",
    "def MakeBinomialMix2(pmf):\n",
    "    binomials = Pmf()\n",
    "    for x, prob in pmf.Items():\n",
    "        binom = MakeBinomialPmf(n=25, p=x)\n",
    "        binomials[binom] = prob\n",
    "    return MakeMixture(binomials)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's how we use it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "rhode_rematch = MakeBinomialMix2(rhode.MakePmf())\n",
    "wei_rematch = MakeBinomialMix2(wei.MakePmf())\n",
    "rhode_rematch.ProbGreater(wei_rematch), rhode_rematch.ProbLess(wei_rematch)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Run this analysis again with a different prior and see how much effect it has on the results."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Distributions of sums and differences\n",
    "\n",
    "Suppose we want to know the total number of targets the two contestants will hit in a rematch.  There are two ways we might compute the distribution of this sum:\n",
    "\n",
    "1. Sampling: We can draw samples from the distributions and add them up.\n",
    "\n",
    "2. Enumeration: We can enumerate all possible pairs of values.\n",
    "\n",
    "I'll start with sampling:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "iters = 1000\n",
    "pmf = Pmf()\n",
    "for _ in range(iters):\n",
    "    k = rhode_rematch.Random() + wei_rematch.Random()\n",
    "    pmf[k] += 1\n",
    "pmf.Normalize()\n",
    "thinkplot.Hist(pmf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or we could use `Sample` and NumPy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ks = rhode_rematch.Sample(iters) + wei_rematch.Sample(iters)\n",
    "pmf = Pmf(ks)\n",
    "thinkplot.Hist(pmf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, we could compute the distribution of the sum by enumeration:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def AddPmfs(pmf1, pmf2):\n",
    "    pmf = Pmf()\n",
    "    for v1, p1 in pmf1.Items():\n",
    "        for v2, p2 in pmf2.Items():\n",
    "            pmf[v1 + v2] += p1 * p2\n",
    "    return pmf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's how it's used:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pmf = AddPmfs(rhode_rematch, wei_rematch)\n",
    "thinkplot.Pdf(pmf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `Pmf` class provides a `+` operator that does the same thing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pmf = rhode_rematch + wei_rematch\n",
    "thinkplot.Pdf(pmf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:**  The Pmf class also provides the `-` operator, which computes the distribution of the difference in values from two distributions.  Use the distributions from the previous section to compute the distribution of the differential between Rhode and Wei in a rematch.  On average, how many clays should we expect Rhode to win by?  What is the probability that Rhode wins by 10 or more?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Distribution of maximum\n",
    "\n",
    "Suppose Kim Rhode continues to compete in six more Olympics.  What should we expect her best result to be?\n",
    "\n",
    "Once again, there are two ways we can compute the distribution of the maximum:\n",
    "\n",
    "1. Sampling.\n",
    "\n",
    "2. Analysis of the CDF.\n",
    "\n",
    "Here's a simple version by sampling:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "iters = 1000\n",
    "pmf = Pmf()\n",
    "for _ in range(iters):\n",
    "    ks = rhode_rematch.Sample(6)\n",
    "    pmf[max(ks)] += 1\n",
    "pmf.Normalize()\n",
    "thinkplot.Hist(pmf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here's a version using NumPy.  I'll generate an array with 6 rows and 10 columns:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "iters = 1000\n",
    "ks = rhode_rematch.Sample((6, iters))\n",
    "ks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the maximum in each column:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "maxes = np.max(ks, axis=0)\n",
    "maxes[:10]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And then plot the distribution of maximums:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pmf = Pmf(maxes)\n",
    "thinkplot.Hist(pmf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or we can figure it out analytically.  If the maximum is less-than-or-equal-to some value `k`, all 6 random selections must be less-than-or-equal-to `k`, so: \n",
    "\n",
    "$ CDF_{max}(x) = CDF(x)^6 $\n",
    "\n",
    "`Pmf` provides a method that computes and returns this `Cdf`, so we can compute the distribution of the maximum like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pmf = rhode_rematch.Max(6).MakePmf()\n",
    "thinkplot.Hist(pmf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:**  Here's how Pmf.Max works:\n",
    "\n",
    "    def Max(self, k):\n",
    "        \"\"\"Computes the CDF of the maximum of k selections from this dist.\n",
    "\n",
    "        k: int\n",
    "\n",
    "        returns: new Cdf\n",
    "        \"\"\"\n",
    "        cdf = self.MakeCdf()\n",
    "        cdf.ps **= k\n",
    "        return cdf\n",
    "\n",
    "Write a function that takes a Pmf and an integer `n` and returns a Pmf that represents the distribution of the minimum of `k` values drawn from the given Pmf.  Use your function to compute the distribution of the minimum score Kim Rhode would be expected to shoot in six competitions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def Min(pmf, k):\n",
    "    cdf = pmf.MakeCdf()\n",
    "    cdf.ps = 1 - (1-cdf.ps)**k\n",
    "    return cdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "pmf = Min(rhode_rematch, 6).MakePmf()\n",
    "thinkplot.Hist(pmf)"
   ]
  }
 ],
 "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}