Think Bayes solutions: Chapter 4

This notebook presents solutions to exercises in Think Bayes.

Copyright 2016 Allen B. Downey

MIT License: https://opensource.org/licenses/MIT


In [3]:
from __future__ import print_function, division

import numpy as np

import thinkbayes2
from thinkbayes2 import Pmf, Cdf, Suite, Beta
import thinkplot

% matplotlib inline

The Euro problem

Here's a class that represents hypotheses about the probability a coin lands heads.


In [2]:
class Euro(Suite):

    def Likelihood(self, data, hypo):
        """Computes the likelihood of `data` given `hypo`.
        
        data: string 'H' or 'T'
        hypo: probability of heads, 0-100
        
        returns: float
        """
        if hypo == 
        
        x = hypo
        if data == 'H':
            return x/100
        else:
            return 1 - x/100


  File "<ipython-input-2-fbdaef23fd42>", line 11
    if hypo ==
               ^
SyntaxError: invalid syntax

We can make a uniform prior and update it with 140 heads and 110 tails:


In [ ]:
suite = Euro(range(0, 101))
dataset = 'H' * 140 + 'T' * 110

for data in dataset:
    suite.Update(data)

And here's what the posterior looks like.


In [ ]:
thinkplot.Pdf(suite)

We can summarize the posterior several ways, including the mean:


In [ ]:
suite.Mean()

Median:


In [ ]:
suite.Percentile(50)

The peak of the posterior, known as the Maximum Aposteori Probability (MAP)


In [ ]:
suite.MAP()

And a 90% credible interval


In [ ]:
suite.CredibleInterval(90)

We can look up a particular value in the posterior PMF, but the result doesn't mean much, because we could have divided the range (0-100) into as many pieces as we like, and the result would be different.


In [ ]:
suite.Prob(50)

Different priors

Let's see how that looks with different priors.

Here's a function that makes a uniform prior:


In [3]:
def UniformPrior(label='uniform'):
    """Makes a Suite with a uniform prior."""
    suite = Euro(range(0, 101), label=label)
    return suite

And another that makes a triangular prior.


In [4]:
def TrianglePrior(label='triangle'):
    """Makes a Suite with a triangle prior."""
    suite = Euro(label=label)
    for x in range(0, 51):
        suite[x] = x
    for x in range(51, 101):
        suite[x] = 100-x 
    suite.Normalize()
    return suite

Here's what they look like:


In [5]:
triangle = TrianglePrior()
uniform = UniformPrior()
suites = [triangle, uniform]

thinkplot.Pdfs(suites)
thinkplot.Config(xlabel='x', ylabel='Probability')


---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-5-0fbdc9ee8111> in <module>()
----> 1 triangle = TrianglePrior()
      2 uniform = UniformPrior()
      3 suites = [triangle, uniform]
      4 
      5 thinkplot.Pdfs(suites)

<ipython-input-4-7f8dd9492dd9> in TrianglePrior(label)
      1 def TrianglePrior(label='triangle'):
      2     """Makes a Suite with a triangle prior."""
----> 3     suite = Euro(label=label)
      4     for x in range(0, 51):
      5         suite[x] = x

NameError: name 'Euro' is not defined

If we update them both with the same data:


In [6]:
def RunUpdate(suite, heads=140, tails=110):
    """Updates the Suite with the given number of heads and tails.

    suite: Suite object
    heads: int
    tails: int
    """
    dataset = 'H' * heads + 'T' * tails
    for data in dataset:
        suite.Update(data)

In [7]:
for suite in suites:
    RunUpdate(suite)


---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-7-428a3b895282> in <module>()
----> 1 for suite in suites:
      2     RunUpdate(suite)

NameError: name 'suites' is not defined

The results are almost identical; the remaining difference is unlikely to matter in practice.


In [8]:
thinkplot.Pdfs(suites)
thinkplot.Config(xlabel='x', ylabel='Probability')


---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-8-bd09bae7fd70> in <module>()
----> 1 thinkplot.Pdfs(suites)
      2 thinkplot.Config(xlabel='x', ylabel='Probability')

NameError: name 'suites' is not defined

The binomial likelihood function

We can make the Euro class more efficient by computing the likelihood of the entire dataset at once, rather than one coin toss at a time.

If the probability of heads is p, we can compute the probability of k=140 heads in n=250 tosses using the binomial PMF.


In [9]:
class Euro2(thinkbayes2.Suite):
    """Represents hypotheses about the probability of heads."""

    def Likelihood(self, data, hypo):
        """Computes the likelihood of the data under the hypothesis.

        hypo: integer value of x, the probability of heads (0-100)
        data: tuple of (number of heads, number of tails)
        """
        x = hypo / 100.0
        heads, tails = data
        like = x**heads * (1-x)**tails
        return like

I left out the binomial coefficient ${n}\choose{k}$ because it does not depend on p, so it's the same for all hypotheses.


In [10]:
suite = Euro2(range(0, 101))
dataset = 140, 110
suite.Update(dataset)


Out[10]:
2.6088321798736822e-76

Here's what the posterior looks like.


In [11]:
thinkplot.Pdf(suite)


The Beta distribution

The Beta distribution is a conjugate prior for the binomial likelihood function, which means that if you start with a Beta distribution and update with a binomial likelihood, the posterior is also Beta.

Also, given the parameters of the prior and the data, we can compute the parameters of the posterior directly. The following class represents a Beta distribution and provides a constant-time Update method.


In [12]:
from scipy import special

class Beta:
    """Represents a Beta distribution.

    See http://en.wikipedia.org/wiki/Beta_distribution
    """
    def __init__(self, alpha=1, beta=1, label=None):
        """Initializes a Beta distribution."""
        self.alpha = alpha
        self.beta = beta
        self.label = label if label is not None else '_nolegend_'

    def Update(self, data):
        """Updates a Beta distribution.

        data: pair of int (heads, tails)
        """
        heads, tails = data
        self.alpha += heads
        self.beta += tails

    def Mean(self):
        """Computes the mean of this distribution."""
        return self.alpha / (self.alpha + self.beta)

    def MAP(self):
        """Computes the value with maximum a posteori probability."""
        a = self.alpha - 1
        b = self.beta - 1
        return a / (a + b)

    def Random(self):
        """Generates a random variate from this distribution."""
        return random.betavariate(self.alpha, self.beta)

    def Sample(self, n):
        """Generates a random sample from this distribution.

        n: int sample size
        """
        size = n,
        return np.random.beta(self.alpha, self.beta, size)

    def EvalPdf(self, x):
        """Evaluates the PDF at x."""
        return x ** (self.alpha - 1) * (1 - x) ** (self.beta - 1)

    def MakePmf(self, steps=101, label=None):
        """Returns a Pmf of this distribution.

        Note: Normally, we just evaluate the PDF at a sequence
        of points and treat the probability density as a probability
        mass.

        But if alpha or beta is less than one, we have to be
        more careful because the PDF goes to infinity at x=0
        and x=1.  In that case we evaluate the CDF and compute
        differences.

        The result is a little funny, because the values at 0 and 1
        are not symmetric.  Nevertheless, it is a reasonable discrete
        model of the continuous distribution, and behaves well as
        the number of values increases.
        """
        if label is None and self.label is not None:
            label = self.label

        if self.alpha < 1 or self.beta < 1:
            cdf = self.MakeCdf()
            pmf = cdf.MakePmf()
            return pmf

        xs = [i / (steps - 1) for i in range(steps)]
        probs = [self.EvalPdf(x) for x in xs]
        pmf = Pmf(dict(zip(xs, probs)), label=label)
        return pmf

    def MakeCdf(self, steps=101):
        """Returns the CDF of this distribution."""
        xs = [i / (steps - 1) for i in range(steps)]
        ps = special.betainc(self.alpha, self.beta, xs)
        cdf = Cdf(xs, ps)
        return cdf

    def Percentile(self, ps):
        """Returns the given percentiles from this distribution.

        ps: scalar, array, or list of [0-100]
        """
        ps = np.asarray(ps) / 100
        xs = special.betaincinv(self.alpha, self.beta, ps)
        return xs

Here's how we use it.


In [13]:
beta = Beta()
beta.Update((140, 110))
beta.Mean()


Out[13]:
0.5595238095238095

And here's the posterior.


In [14]:
thinkplot.Pdf(beta.MakePmf())


Amazing, no?

Exercise: One way to construct priors is to make a Beta distribution and adjust the parameters until it has the shape you want. Then when you do an update, the data get added to the parameters of the prior. Since the parameters of the prior play the same mathematical role as the data, they are sometimes called "precounts".

Suppose you believe that most coins are fair or unlikely to deviate from 50% by more than a few percentage points. Construct a prior that captures this belief and update it with the Euro data. How much effect does it have on the posterior, compared to the uniform prior?

Hint: A Beta distribution with parameters (1, 1) is uniform from 0 to 1.


In [15]:
beta2 = Beta(200,200)
thinkplot.Pdf(beta2.MakePmf())



In [16]:
beta2.Update((140, 110))
thinkplot.Pdf(beta2.MakePmf())
beta2.Mean()


Out[16]:
0.5230769230769231

In [17]:
# Solution goes here

In [18]:
# Solution goes here

In [19]:
# Solution goes here

Exercise: 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 skeets, sending the match into sudden death. In the first round, both hit 1 of 2 skeets. In the next two rounds, they each hit 2 skeets. 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?

As always, you will have to make some modeling decisions, but one approach is to estimate, for each shooter, the probability of hitting a skeet. Then, to estimate the probability that Rhode is a better shooter, you can draw samples from the two posterior distributions and compare them. To estimate the probability of winning a rematch, you could draw samples from the posterior distributions and simulate a round of 25 shots.


In [4]:
betaRhode = Beta(4,2)
betaWei = Beta(4,2)
thinkplot.Pdf(betaRhode.MakePmf())



In [21]:
betaRhode.Update((15+1+2+2, 10+1+0))
betaWei.Update((15+1+2+1, 10+1+1))

In [22]:
thinkplot.Pdf(betaRhode.MakePmf())
thinkplot.Pdf(betaWei.MakePmf())


Now, we will estimate the probability that Rhode is actually a better shooter than Wei. To do this, we will make a cumulative distribution for Rhode and Wei's beta distributions.


In [23]:
rhodeSample = betaRhode.MakeCdf(500).Sample(1000)
weiSample = betaWei.MakeCdf(500).Sample(1000)

In [24]:
np.mean(rhodeSample > weiSample)


Out[24]:
0.61599999999999999

With a relatively small number of steps in the cdf (500) and small number of samples (1000), we see that Rhode is not that likely to be a better shooter than Wei. We assign a ~54% chance that Rhode is a better shooter then Wei.

Next, we calculate the chance that Rhode will win if the two play again.


In [25]:
rhodeRematch = np.random.binomial(25, rhodeSample)

In [26]:
print(rhodeRematch)


[ 8 17 17 24 17 19 15 19 22 12 15 14 19 13 14 16 22 16 19 11 13 16 12 16 15
 15 21 13 14 17 17 15 20 12 18 14 19 15 19 16 11 18 19 16 21 20 13 11 18 16
 17 16 18 18 19 11  8 15 19 14 19 15 20 19 21 17 18 12 21 10 20 16 13 20 17
 11 21 16 19 18 21 12 15 12  8 12 14 15 17 19 21 15 22 13 13 18 14 21 13 15
 18 22 19 16 14 17 22 18 16 15 15 17 13 18 18 16 16 13 13 19 19 20 14 15 12
 16 15 11 18 16 16 14 18 16 15 15 18 17 15 15 14 20 15 13 19 12 16 14 17 16
 20 15 21 22 15 18 17 17 11 16 17 16 18 16 16 18 19 15 15 15 11 17 18 16 18
 20 15 18 20 22  9 14 12 17 21 20 19 16 18 21 15 18 10 17 11 19 12 19 18 19
 21 16 14 21 19 14 17 20 15 21 17 15 21 18 17 18 18 17 12 17 19 16 21 18 14
 16 18 21 13 16 17 11 19 20 19 14 16 19 14 14 21 22 14 17 18 14 15 18 10 21
 13 22 18 17 11 19 10 14 14 17 18 16 19 14 19 19 10 21 20 14 20 12 17 16 13
 15 16 18 14 18 12 18 17 15 18 16 22 17 17 10 19 19 18 21 16 17 15 17 17 22
 17 19 18 20 16 18 18 15 20 13 16  8 16 18 22 15 17 16 13 18 16  8 18 12 17
 17 19 15 20 14 15 17 22 19 21 19 17 18 16 13 11 16 17 19 17 11 15 15 17 16
 16 11 16 10 14 17 16 19 19 14 13 14  8 21 19 16 15 13 20 14 18 23 14 12 13
 17 10 17 18 15 16 14 19 11 19 14 20  9 18 14 12 16 16 12 14 13 15 14 21 21
 18 14 13 16 18 13 19 16 19 15 16 15 17 19 20 17 15 15 18 14 11 22 17 18 12
 13 17 13 16 10 16 14 16 14 15 12 22 20 19 18 13 19 13 15 17 16 17 15 15 20
 17 11 12 15 17  9 12 19 19  7 10 17 15 15 17 15 13 18 17 17 19 12 14 16 16
 19 14 25 14 14 18 20 23 20 16 19 19 12 19 14 16 17 12 15 18 19 18 15 14 15
 14 13 16 17 12 13 15 18 16 16 17 12 17 18 16 19 15 18 14 12 16 21 19 17 20
 17 13 17 15 18 20 21 12 16 15 18  9 16 18 11 10  8 13 17 21 15 11 11 15 20
 16 17 14 11 16 15 10 18 13 15 18 14 19 11 18 17 18 24 20 15 14 16 15 11 12
 18 17 15 11 11 19 15 14 17 16 19 21 13 17 14 19 19 15 18 19 21 22 13 16 12
 19 18 17 12 20 24 16 14 16 11 15 13 19 19 10 15 18 16 14 17 15 17 14 17 20
 19 16 21 19 17 13 15 15 12 14 19 14 18 17 15 15 18 13 15 14 10 17 19 15 20
 16 15 24 18 15 13 14 17 17 15 18 16 15 18 18 12 21 14 19 17 14 19 18 12 13
 21 15 18 16 19 16 16 15 16 20 15  8 15 12 15 18 21 17 12 18 18 20 16 19 14
 17 22 16 16 12 16 18 19  9 19 12 20 17  9 18 20 15 18 20 10 11 21 23 15 19
 15  9 17 20 12  9 15 16 19 13 12 17 12 15 18 20 20 17 14 16 15 12 15 15 22
 14 14 13 17 13 18 22 23 16 15 20 19 18 17 13 13 14 19 16 12 16 17 21 19 14
 16 16 25 21 19 21 14 21 17 14 17 10 19 19 13 23 13 14 18 14 13 18 18 22 17
 15 16 16 19 20 14 12 10 17 11 18 18 18 18 18 15 15 19 23 12 21 17 14 17 16
 17 15 12 17 17 17 15 14 19 21 17 16 19 12 20 17 15 19 21 20 12 13 17 16 21
 20 20 13 18 16 19 19 17 14 15 14 13 17 17 19 18 14 14 16 14 16 17 12 14 13
 13 19 12 14 14 15 17 17 11 16 20 13 19 14 13 15 14 20 14 12 16 20 18 12 17
 17 19 12 17 22 16 20 13 11 12 17 14 19 18 18 14 16 18 16 18 15 15 18 16  9
 14  9 21 20 20 14 12  9 20 17 16 17 19 20 12 23 18 18 19 15 15 17 18 19 17
 15 18 17 19 14 16 18 18 19 19 17 12 16 18 14 15 15 17 17 11 13 15 21 20 21
 10 10 13 20 17 16 20 20 13 13 17 16 12 15 12 23 15 17 16 14 13 16 17 14 16]

Exercise Suppose that instead of observing coin tosses directly, you measure the outcome using an instrument that is not always correct. Specifically, suppose there is a probability y that an actual heads is reported as tails, or actual tails reported as heads.

Write a class that estimates the bias of a coin given a series of outcomes and the value of y.

How does the spread of the posterior distribution depend on y?


In [ ]:
# Solution goes here

In [ ]:
# Solution goes here

In [ ]:
# Solution goes here

In [ ]:
# Solution goes here

In [ ]:
# Solution goes here

In [ ]:
# Solution goes here

In [ ]:
# Solution goes here

Exercise This exercise is inspired by a question posted by a “redditor” named dominosci on Reddit’s statistics “subreddit” at http://reddit.com/r/statistics.

Reddit is an online forum with many interest groups called subreddits. Users, called redditors, post links to online content and other web pages. Other redditors vote on the links, giving an “upvote” to high-quality links and a “downvote” to links that are bad or irrelevant.

A problem, identified by dominosci, is that some redditors are more reliable than others, and Reddit does not take this into account.

The challenge is to devise a system so that when a redditor casts a vote, the estimated quality of the link is updated in accordance with the reliability of the redditor, and the estimated reliability of the redditor is updated in accordance with the quality of the link.

One approach is to model the quality of the link as the probability of garnering an upvote, and to model the reliability of the redditor as the probability of correctly giving an upvote to a high-quality item.

Write class definitions for redditors and links and an update function that updates both objects whenever a redditor casts a vote.


In [ ]:
# Solution goes here

In [ ]:
# Solution goes here

In [ ]:
# Solution goes here

In [ ]:
# Solution goes here

In [ ]:
# Solution goes here

In [ ]:
# Solution goes here

In [ ]:
# Solution goes here

In [ ]:
# Solution goes here

In [ ]: