<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Think Bayes
This notebook presents example code and exercise solutions for Think Bayes.
Copyright 2018 Allen B. Downey
MIT License
Step1:  Comparing distributions
Let's get back to the Kim Rhode problem from Chapter 4
Step2:  Based on the data, the distribution for Rhode is slightly farther right than the distribution for Wei, but there is a lot of overlap.
Step3:  To compute the probability that Rhode actually has a higher value of p, there are two options
Step4:  Beta also provides Sample, which returns a NumPy array, so we an perform the comparisons using array operations
Step5:  The other option is to make Pmf objects that approximate the Beta distributions, and enumerate pairs of values
Step6:  Exercise
Step7:  flip returns True with probability p and False with probability 1-p
Now we can simulate 1000 rematches and count wins and losses.
Step8:  Or, realizing that the distribution of k is binomial, we can simplify the code using NumPy
Step9:  Alternatively, we can make a mixture that represents the distribution of k, taking into account our uncertainty about x
Step10:  Alternatively, we could use MakeMixture
Step11:  Here's how we use it.
Step12:  Exercise
Step13:  Or we could use Sample and NumPy
Step14:  Alternatively, we could compute the distribution of the sum by enumeration
Step15:  Here's how it's used
Step16:  The Pmf class provides a + operator that does the same thing.
Step17:  Exercise
Step18:  Distribution of maximum
Suppose Kim Rhode continues to compete in six more Olympics.  What should we expect her best result to be?
Once again, there are two ways we can compute the distribution of the maximum
Step19:  And here's a version using NumPy.  I'll generate an array with 6 rows and 10 columns
Step20:  Compute the maximum in each column
Step21:  And then plot the distribution of maximums
Step22:  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
Step23:  Exercise
<ASSISTANT_TASK:>
Python Code:
# Configure Jupyter so figures appear in the notebook
%matplotlib inline
# Configure Jupyter to display the assigned value after an assignment
%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'
# import classes from thinkbayes2
from thinkbayes2 import Hist, Pmf, Suite, Beta
import thinkplot
import numpy as np
Explanation: Think Bayes
This notebook presents example code and exercise solutions for Think Bayes.
Copyright 2018 Allen B. Downey
MIT License: https://opensource.org/licenses/MIT
End of explanation
rhode = Beta(1, 1, label='Rhode')
rhode.Update((22, 11))
wei = Beta(1, 1, label='Wei')
wei.Update((21, 12))
Explanation: Comparing distributions
Let's get back to the Kim Rhode problem from Chapter 4:
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.
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?
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.
First I create a Beta distribution for each of the competitors, and update it with the results.
End of explanation
thinkplot.Pdf(rhode.MakePmf())
thinkplot.Pdf(wei.MakePmf())
thinkplot.Config(xlabel='x', ylabel='Probability')
Explanation: Based on the data, the distribution for Rhode is slightly farther right than the distribution for Wei, but there is a lot of overlap.
End of explanation
iters = 1000
count = 0
for _ in range(iters):
    x1 = rhode.Random()
    x2 = wei.Random()
    if x1 > x2:
        count += 1
count / iters
Explanation: To compute the probability that Rhode actually has a higher value of p, there are two options:
Sampling: we could draw random samples from the posterior distributions and compare them.
Enumeration: we could enumerate all possible pairs of values and add up the "probability of superiority".
I'll start with sampling.  The Beta object provides a method that draws a random value from a Beta distribution:
End of explanation
rhode_sample = rhode.Sample(iters)
wei_sample = wei.Sample(iters)
np.mean(rhode_sample > wei_sample)
Explanation: Beta also provides Sample, which returns a NumPy array, so we an perform the comparisons using array operations:
End of explanation
def ProbGreater(pmf1, pmf2):
    total = 0
    for x1, prob1 in pmf1.Items():
        for x2, prob2 in pmf2.Items():
            if x1 > x2:
                total += prob1 * prob2
    return total
pmf1 = rhode.MakePmf(1001)
pmf2 = wei.MakePmf(1001)
ProbGreater(pmf1, pmf2)
pmf1.ProbGreater(pmf2)
pmf1.ProbLess(pmf2)
Explanation: The other option is to make Pmf objects that approximate the Beta distributions, and enumerate pairs of values:
End of explanation
import random
def flip(p):
    return random.random() < p
Explanation: Exercise: Run this analysis again with a different prior and see how much effect it has on the results.
Simulation
To make predictions about a rematch, we have two options again:
Sampling.  For each simulated match, we draw a random value of x for each contestant, then simulate 25 shots and count hits.
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.
I'll do it by sampling first.
End of explanation
iters = 1000
wins = 0
losses = 0
for _ in range(iters):
    x1 = rhode.Random()
    x2 = wei.Random()
    
    count1 = count2 = 0
    for _ in range(25):
        if flip(x1):
            count1 += 1
        if flip(x2):
            count2 += 1
            
    if count1 > count2:
        wins += 1
    if count1 < count2:
        losses += 1
        
wins/iters, losses/iters
Explanation: flip returns True with probability p and False with probability 1-p
Now we can simulate 1000 rematches and count wins and losses.
End of explanation
rhode_rematch = np.random.binomial(25, rhode_sample)
thinkplot.Hist(Pmf(rhode_rematch))
wei_rematch = np.random.binomial(25, wei_sample)
np.mean(rhode_rematch > wei_rematch)
np.mean(rhode_rematch < wei_rematch)
Explanation: Or, realizing that the distribution of k is binomial, we can simplify the code using NumPy:
End of explanation
from thinkbayes2 import MakeBinomialPmf
def MakeBinomialMix(pmf, label=''):
    mix = Pmf(label=label)
    for x, prob in pmf.Items():
        binom = MakeBinomialPmf(n=25, p=x)
        for k, p in binom.Items():
            mix[k] += prob * p
    return mix
rhode_rematch = MakeBinomialMix(rhode.MakePmf(), label='Rhode')
wei_rematch = MakeBinomialMix(wei.MakePmf(), label='Wei')
thinkplot.Pdf(rhode_rematch)
thinkplot.Pdf(wei_rematch)
thinkplot.Config(xlabel='hits')
rhode_rematch.ProbGreater(wei_rematch), rhode_rematch.ProbLess(wei_rematch)
Explanation: Alternatively, we can make a mixture that represents the distribution of k, taking into account our uncertainty about x:
End of explanation
from thinkbayes2 import MakeMixture
def MakeBinomialMix2(pmf):
    binomials = Pmf()
    for x, prob in pmf.Items():
        binom = MakeBinomialPmf(n=25, p=x)
        binomials[binom] = prob
    return MakeMixture(binomials)
Explanation: Alternatively, we could use MakeMixture:
End of explanation
rhode_rematch = MakeBinomialMix2(rhode.MakePmf())
wei_rematch = MakeBinomialMix2(wei.MakePmf())
rhode_rematch.ProbGreater(wei_rematch), rhode_rematch.ProbLess(wei_rematch)
Explanation: Here's how we use it.
End of explanation
iters = 1000
pmf = Pmf()
for _ in range(iters):
    k = rhode_rematch.Random() + wei_rematch.Random()
    pmf[k] += 1
pmf.Normalize()
thinkplot.Hist(pmf)
Explanation: Exercise: Run this analysis again with a different prior and see how much effect it has on the results.
Distributions of sums and differences
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:
Sampling: We can draw samples from the distributions and add them up.
Enumeration: We can enumerate all possible pairs of values.
I'll start with sampling:
End of explanation
ks = rhode_rematch.Sample(iters) + wei_rematch.Sample(iters)
pmf = Pmf(ks)
thinkplot.Hist(pmf)
Explanation: Or we could use Sample and NumPy:
End of explanation
def AddPmfs(pmf1, pmf2):
    pmf = Pmf()
    for v1, p1 in pmf1.Items():
        for v2, p2 in pmf2.Items():
            pmf[v1 + v2] += p1 * p2
    return pmf
Explanation: Alternatively, we could compute the distribution of the sum by enumeration:
End of explanation
pmf = AddPmfs(rhode_rematch, wei_rematch)
thinkplot.Pdf(pmf)
Explanation: Here's how it's used:
End of explanation
pmf = rhode_rematch + wei_rematch
thinkplot.Pdf(pmf)
Explanation: The Pmf class provides a + operator that does the same thing.
End of explanation
# Solution
pmf = rhode_rematch - wei_rematch
thinkplot.Pdf(pmf)
# Solution
# On average, we expect Rhode to win by about 1 clay.
pmf.Mean(), pmf.Median(), pmf.Mode()
# Solution
# But there is, according to this model, a 2% chance that she could win by 10.
sum([p for (x, p) in pmf.Items() if x >= 10])
Explanation: 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?
End of explanation
iters = 1000
pmf = Pmf()
for _ in range(iters):
    ks = rhode_rematch.Sample(6)
    pmf[max(ks)] += 1
pmf.Normalize()
thinkplot.Hist(pmf)
Explanation: Distribution of maximum
Suppose Kim Rhode continues to compete in six more Olympics.  What should we expect her best result to be?
Once again, there are two ways we can compute the distribution of the maximum:
Sampling.
Analysis of the CDF.
Here's a simple version by sampling:
End of explanation
iters = 1000
ks = rhode_rematch.Sample((6, iters))
ks
Explanation: And here's a version using NumPy.  I'll generate an array with 6 rows and 10 columns:
End of explanation
maxes = np.max(ks, axis=0)
maxes[:10]
Explanation: Compute the maximum in each column:
End of explanation
pmf = Pmf(maxes)
thinkplot.Hist(pmf)
Explanation: And then plot the distribution of maximums:
End of explanation
pmf = rhode_rematch.Max(6).MakePmf()
thinkplot.Hist(pmf)
Explanation: 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: 
$ CDF_{max}(x) = CDF(x)^6 $
Pmf provides a method that computes and returns this Cdf, so we can compute the distribution of the maximum like this:
End of explanation
def Min(pmf, k):
    cdf = pmf.MakeCdf()
    cdf.ps = 1 - (1-cdf.ps)**k
    return cdf
pmf = Min(rhode_rematch, 6).MakePmf()
thinkplot.Hist(pmf)
Explanation: Exercise:  Here's how Pmf.Max works:
def Max(self, k):
    Computes the CDF of the maximum of k selections from this dist.
    k: int
    returns: new Cdf
    
    cdf = self.MakeCdf()
    cdf.ps **= k
    return cdf
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.
End of explanation
<END_TASK>