Bayesian Statistics
Made Simple
Allen B. Downey
Olin College
sites.google.com/site/simplebayes
Follow along at home
sites.google.com/site/simplebayes
The plan
From Bayes's Theorem to Bayesian inference.
A computational framework.
Work on example problems.
Goals
At 4:40, you should be ready to:
Think Bayes
This tutorial is based on my book,
Think Bayes
Bayesian Statistics in Python
Published by O'Reilly Media
and available under a
Creative Commons license from
thinkbayes.com
Probability
p(A): the probability that A occurs.
p(A|B): the probability that A occurs, given that B has occurred.
p(A and B) = p(A) p(B|A)
Bayes's Theorem
By definition of conjoint probability:
p(A and B) = p(A) p(B|A) = (1)
p(B and A) = p(B) p(A|B)
Equate the right hand sides
p(B) p(A|B) = p(A) p(B|A) (2)
Divide by p(B) and ...
Bayes's Theorem
Bayes's Theorem
One way to think about it:
Bayes's Theorem is an algorithm to get
from p(B|A) to p(A|B).
Useful if p(B|A), p(A) and p(B) are easier
than p(A|B).
OR ...
Diachronic interpretation
H: Hypothesis
D: Data
Given p(H), the probability of the hypothesis before you saw the data.
Find p(H|D), the probability of the hypothesis after you saw the data.
A cookie problem
Suppose there are two bowls of cookies.
Bowl #1 has 10 chocolate and 30 vanilla.
Bowl #2 has 20 of each.
Fred picks a bowl at random, and then picks a cookie at random. The cookie turns out to be vanilla.
What is the probability that Fred picked from Bowl #1?
from Wikipedia
Cookie problem
H: Hypothesis that cookie came from Bowl 1.
D: Cookie is vanilla.
Given p(H), the probability of the hypothesis before you saw the data.
Find p(H|D), the probability of the hypothesis after you saw the data.
Diachronic interpretation
p(H|D) = p(H) p(D|H) / p(D)
p(H): prior
p(D|H): conditional likelihood of the data
p(D): total likelihood of the data
Diachronic interpretation
p(H|D) = p(H) p(D|H) / p(D)
p(H): prior = 1/2
p(D|H): conditional likelihood of the data = 3/4
p(D): total likelihood of the data = 5/8
Diachronic interpretation
p(H|D) = (1/2)(3/4) / (5/8) = 3/5
p(H): prior = 1/2
p(D|H): conditional likelihood of the data = 3/4
p(D): total likelihood of the data = 5/8
A little intuition
p(H|D): posterior = 60%
p(H): prior = 50%
Vanilla cookie was more likely under H.
Slightly increases our degree of belief in H.
Computation
Pmf represents a Probability Mass Function
Maps from possible values to probabilities.
Diagram by yuml.me
Install test
How many of you got install_test.py running?
Don't try to fix it now!
Instead...
Partner up
Icebreaker
What was your first computer?
What was your first programming language?
What is the longest time you have spent finding a stupid bug?
Start your engines
$ python
>>> from thinkbayes import Pmf
Pmf
from thinkbayes import Pmf
# make an empty Pmf
d6 = Pmf()
# outcomes of a six-sided die
for x in [1,2,3,4,5,6]:
d6.Set(x, 1)
Pmf
d6.Print()
d6.Normalize()
d6.Print()
d6.Random()
Style
I know what you're thinking.
The Bayesian framework
1) Build a Pmf that maps from each hypothesis to a prior probability, p(H).
2) Multiply each prior probability by the likelihood of the data, p(D|H).
3) Normalize, which divides through by the total likelihood, p(D).
Prior
pmf = Pmf()
pmf.Set('Bowl 1', 0.5)
pmf.Set('Bowl 2', 0.5)
Update
p(Vanilla | Bowl 1) = 30/40
p(Vanilla | Bowl 2) = 20/40
pmf.Mult('Bowl 1', 0.75)
pmf.Mult('Bowl 2', 0.5)
Normalize
pmf.Normalize()
0.625 # return val is p(D)
print pmf.Prob('Bowl 1')
0.6
Exercise
What if we select another cookie, and it’s chocolate?
The posterior (after the first cookie) becomes the prior (before the second cookie).
Exercise
What if we select another cookie, and it’s chocolate?
pmf.Mult('Bowl 1', 0.25)
pmf.Mult('Bowl 2', 0.5)
pmf.Normalize()
pmf.Print()
Bowl 1 0.43�Bowl 2 0.573�
Summary
Bayes's Theorem,
Cookie problem,
Pmf class.
The dice problem
I have a box of dice that contains a 4-sided die, a 6-sided die, an 8-sided die, a 12-sided die and a 20-sided die.
Suppose I select a die from the box at random, roll it, and get a 6. What is the probability that I rolled each die?
Hypothesis suites
A suite is a mutually exclusive and collectively exhaustive set of hypotheses.
Represented by a Suite that maps
hypothesis → probability.
Suite
class Suite(Pmf):
"Represents a suite of hypotheses and
their probabilities."
def __init__(self, hypos):
"Initializes the distribution."
for hypo in hypos:
self.Set(hypo, 1)
self.Normalize()
Suite
def Update(self, data):
"Updates the suite based on data."
for hypo in self.Values():� like = self.Likelihood(data, hypo)� self.Mult(hypo, like)�
self.Normalize()�
self.Likelihood?
Suite
Likelihood is an abstract method.
Child classes inherit Update,
provide Likelihood.
Likelihood
Outcome: 6
What is the likelihood of getting n on an m-sided die?
Likelihood
# hypo is the number of sides on the die
# data is the outcome
class Dice(Suite):
def Likelihood(self, data, hypo):� # write this method!�
Write your solution in dice.py
Likelihood
# hypo is the number of sides on the die
# data is the outcome
class Dice(Suite):
def Likelihood(self, data, hypo):� if hypo < data:� return 0� else:� return 1.0/hypo�
Dice
# start with equal priors
suite = Dice([4, 6, 8, 12, 20])
# update with the data
suite.Update(6)
suite.Print()
Dice
Posterior distribution:
4 0.0
6 0.39
8 0.30
12 0.19
20 0.12
More data? No problem...
Dice
for roll in [8, 7, 7, 5, 4]:
suite.Update(roll)
suite.Print()
Dice
Posterior distribution:
4 0.0
6 0.0
8 0.92
12 0.080
20 0.0038
Summary
Dice problem,
Likelihood function,
Suite class.
Recess!
http://ksdcitizens.org/wp-content/uploads/2010/09/recess_time.jpg
Trains
The trainspotting problem:
Modify train.py to compute your answer.
Trains
Train
print suite.Mean()
print suite.MaximumLikelihood()�print suite.CredibleInterval(90)�
Trains
Trains
Tanks
Good time for questions
A Euro problem
"When spun on edge 250 times, a Belgian one-euro coin came up heads 140 times and tails 110. 'It looks very suspicious to me,' said Barry Blight, a statistics lecturer at the London School of Economics. 'If the coin were unbiased, the chance of getting a result as extreme as that would be less than 7%.' "
From "The Guardian" quoted by MacKay, Information Theory, Inference, and Learning Algorithms.
A Euro problem
MacKay asks, "But do these data give evidence that the coin is biased rather than fair?"
Assume that the coin has probability x of landing heads.
(Forget that x is a probability; just think of it as a physical characteristic.)
A Euro problem
Estimation: Based on the data (140 heads, 110 tails), what is x?
Hypothesis testing: What is the probability that the coin is fair?
Euro
We can use the Suite template again.
We just have to figure out the likelihood function.
Likelihood
# hypo is the prob of heads (1-100)
# data is a string, either 'H' or 'T'
class Euro(Suite):�� def Likelihood(self, data, hypo):� # one more, please!
Modify euro.py to compute your answer.
Likelihood
# hypo is the prob of heads (1-100)
# data is a string, either 'H' or 'T'
class Euro(Suite):�� def Likelihood(self, data, hypo):� x = hypo / 100.0� if data == 'H':� return x� else:� return 1-x
Prior
What do we believe about x before seeing the data?
Start with something simple; we'll come back and review.
Uniform prior: any value of x between 0% and 100% is equally likely.
Prior
suite = Euro(range(0, 101))
Update
Suppose we spin the coin once and get heads.
suite.Update('H')
What does the posterior distribution look like?
Hint: what is p(x=0% | D)?
Update
Suppose we spin the coin again, and get heads again.
suite.Update('H')
What does the posterior distribution look like?
Update
Suppose we spin the coin again, and get tails.
suite.Update('T')
What does the posterior distribution look like?
Hint: what's p(x=100% | D)?
Update
After 10 spins, 7 heads and 3 tails:
for outcome in 'HHHHHHHTTT':
suite.Update(outcome)
Update
And finally, after 140 heads and 110 tails:
evidence = 'H' * 140 + 'T' * 110
for outcome in evidence:
suite.Update(outcome)
Posterior
Posterior
Given the posterior distribution, what is the probability that x is 50%?
print suite.Prob(50)
And the answer is... 0.021
Hmm. Maybe that's not the right question.
Posterior
How about the most likely value of x?
print pmf.MaximumLikelihood()
And the answer is 56%.
Posterior
Or the expected value?
print suite.Mean()
And the answer is 55.95%.
Posterior
Credible interval?
print suite.CredibleInterval(90)
Posterior
The 5th percentile is 51.
The 95th percentile is 61.
These values form a 90% credible interval.
So can we say: "There's a 90% chance that x is between 51 and 61?"
Frequentist response
Thank you smbc-comics.com
Bayesian response
Yes, x is a random variable,
Yes, (51, 61) is a 90% credible interval,
Yes, x has a 90% chance of being in it.
Pro: Bayesian stats lends itself to decision analysis.
Con: The prior is subjective.
The prior is subjective
Remember the prior?
We chose it pretty arbitrarily, and reasonable people might disagree.
Is x as likely to be 1% as 50%?
Given what we know about coins, I doubt it.
Prior
How should we capture background knowledge about coins?
Try a triangle prior.
Posterior
What do you think the posterior distributions look like?
I was going to put an image here, but then I Googled "posterior". Never mind.
Swamp the prior
With enough data, reasonable people converge.
But if any p(Hi) = 0, no data will change that.
Swamp the prior
Priors can be arbitrarily low, but avoid 0.
"I beseech you, in the bowels of Christ, think it possible that you may be mistaken."�
Summary of estimation
Recess!
http://images2.wikia.nocookie.net/__cb20120116013507/recess/images/4/4c/Recess_Pic_for_the_Internet.png
Hypothesis testing
Remember the original question:
"But do these data give evidence that the coin is biased rather than fair?"
What does it mean to say that data give evidence for (or against) a hypothesis?
Hypothesis testing
D is evidence in favor of H if
p(H|D) > p(H)
which is true if
p(D|H) > p(D|~H)
or equivalently if
p(D|H) / p(D|~H) > 1
Hypothesis testing
This term
p(D|H) / p(D|~H)
is called the likelihood ratio, or Bayes factor.
It measures the strength of the evidence.
Hypothesis testing
F: hypothesis that the coin is fair
B: hypothesis that the coin is biased
p(D|F) is easy.
p(D|B) is hard because B is underspecified.
Bogosity
Tempting: we got 140 heads out of 250 spins, so B is the hypothesis that x = 140/250.
But,
We need some rules
Likelihood
def AverageLikelihood(suite, data):
total = 0
for hypo, prob in suite.Items():
like = suite.Likelihood(data, hypo)
total += prob * like
return total
Hypothesis testing
F: hypothesis that x = 50%.
B: hypothesis that x is not 50%, but might be any other value with equal probability.
Prior
fair = Euro()
fair.Set(50, 1)
Prior
bias = Euro()
for x in range(0, 101):
if x != 50:
bias.Set(x, 1)
bias.Normalize()
Bayes factor
data = 140, 110
like_fair = AverageLikelihood(fair, data)
like_bias = AverageLikelihood(bias, data)
ratio = like_bias / like_fair
Hypothesis testing
Read euro2.py.
Notice the new representation of the data, and corresponding Likelihood function.
Run it and interpret the results.
Hypothesis testing
And the answer is:
p(D|B) = 2.6 · 10-76
p(D|F) = 5.5 · 10-76
Likelihood ratio is about 0.47.
So this dataset is evidence against B.
Fair comparison?
bias = Euro()
for x in range(0, 49):
bias.Set(x, x)
for x in range(51, 101):
bias.Set(x, 100-x)
bias.Normalize()
Conclusion
Summary
Euro problem,
Bayesian estimation,
Bayesian hypothesis testing.
Recess!
http://ksdcitizens.org/wp-content/uploads/2010/09/recess_time.jpg
Word problem for geeks
ALICE: What did you get on the math SAT?
BOB: 760
ALICE: Oh, well I got a 780. I guess that means I'm smarter than you.
NARRATOR: Really? What is the probability that Alice is smarter than Bob?
Assume, define, quantify
Assume: each person has some probability, x, of answering a random SAT question correctly.
Define: "Alice is smarter than Bob" means xa > xb.
How can we quantify Prob { xa > xb } ?
Be Bayesian
Treat x as a random quantity.
Start with a prior distribution.
Update it.
Compare posterior distributions.
Prior?
Distribution of raw scores.
Likelihood
def Likelihood(self, data, hypo):� x = hypo� score = data� raw = self.exam.Reverse(score)�� yes, no = raw, self.exam.max_score - raw� like = x**yes * (1-x)**no� return like
Posterior
PmfProbGreater
def PmfProbGreater(pmf1, pmf2):
"""Returns the prob that a value from pmf1
is greater than a value from pmf2."""
PmfProbGreater
def PmfProbGreater(pmf1, pmf2):
"""Returns the prob that a value from pmf1
is greater than a value from pmf2."""
Iterate through all pairs of values.
Check whether the value from pmf1 is greater.
Add up total probability of successful pairs.
PmfProbGreater
def PmfProbGreater(pmf1, pmf2):
for x1, p1 in pmf1.Items():
for x2, p2 in pmf2.Items():
# FILL THIS IN!
PmfProbGreater
def PmfProbGreater(pmf1, pmf2):
total = 0.0
for x1, p1 in pmf1.Items():
for x2, p2 in pmf2.Items():
if x1 > x2:
total += p1 * p2
return total
And the answer is...
Alice: 780
Bob: 760
Probability that Alice is "smarter": 61%
Why this example?
Posterior distribution is often the input to the next step in an analysis.
Real world problems start (and end!) with modeling.
Modeling
Modeling
Modeling
Therefore:
Recess!
http://images2.wikia.nocookie.net/__cb20120116013507/recess/images/4/4c/Recess_Pic_for_the_Internet.png
Relax
The Red Line problem
Kendall Square to South Station, catch the 4:40 train.
Arrive Kendall/MIT at 4:15, usually arrive South Station by 4:30.
Counting passengers
Few passengers waiting: just missed a train. Wait 7 minutes.
Many passengers: it’s been a while. Wait 2 minutes.
Platform packed: something’s wrong. Taxi.
Three inferences
1. Ahead of time, watch passengers arrive, estimate arrival rate.
2. On the day, count passengers on platform, estimate time since last train.
3. Compute the distribution of remaining time, conditioned on time since last train.
Data
Three weeks of train arrivals at Kendall/MIT between 4pm and 5pm.
Observer bias
More people arrive during longer intervals, so the average interval reported by passengers is longer than the average reported by MBTA.
Example: Alternating 5 and 10 minute delays.
MBTA average: 7.5 minutes
Passenger average: 8.33 minutes
Data
Three weeks of train arrivals at Kendall/MIT between 4pm and 5pm.
z = time between trains
zb = time between trains as seen by random arrival
Wait time
z = time between trains
zb = time between trains as seen by random arrival
y = wait time for a random arrival
Example
z: Trains run every 7.5 minutes.
zb: On average, a random passenger arrives during an 8.33-minute interval.
y: How long should the passenger expect to wait?
The Bayesian part
x = time since the last train
prior x: distribution of x before I see the data.
posterior x: distribution of x after.
Example
When I arrive, I see 15 passengers in my measurement area.
Based on previous observation, I expect 2 passengers to arrive per minute.
How long has it been since the last train?
Prediction
z: time between trains.
x: time since last train.
y: time until the next train.
z = x + y
y = z - x
Decision analysis
Having observed num passengers,
1. Compute the predictive distribution of time until next train, y.
2. Compute the probability that y exceeds 15 minutes.
Summary
Bayesian analysis consistent with intuition.
Predictive distributions guide decision analysis.
Almost done!
Think Bayes
This tutorial is based on my book,
Think Bayes
Bayesian Statistics Made Simple
Published by O'Reilly Media
and available under a
Creative Commons license from
thinkbayes.com
Case studies
Think Stats
You might also like
Think Stats
Prob/Stat for Programmers
Published by O'Reilly Media
and available under a
Creative Commons license from
thinkstats.com
More reading
More reading
More reading (not free)
Howson and Urbach, Scientific Reasoning: The Bayesian Approach
Sivia, Data Analysis: A Bayesian Tutorial
Gelman et al, Bayesian Data Analysis
Need help?
I am always looking for interesting projects.
Remember
Where does this fit?
Usual approach:
Where does this fit?
Problem:
My theory