18  Introduction to Bayesian Statistics

Course textbook for F2026 under active development

📄 Download the R code from this chapter

📝 Download the class handout for this chapter

We have spent this entire course learning how to make inferences about populations from samples of data. Every tool we have built – confidence intervals, hypothesis tests, standard errors – follows the same fundamental logic. We consider our test statistic to be drawn from a random variable (the sampling distribution) that has a known shape. We posit something about that sampling distribution (that it is centered on 50% for a poll, or on 0 for the difference in means, or whatever we want), and then ask: how surprising is our test statistic given our assumption?

This approach to statistics is called frequentist statistics. It is so called because the sampling distribution describes the hypothetical distribution of all possible test statistics if we estimated them in many samples. This is the dominant paradigm in the social sciences, and everything you have learned in this course falls under that umbrella. But it is not the only way to think about uncertainty. In this chapter we are going to explore an alternative: Bayesian statistics.

This chapter is not going to make you a Bayesian statistician. That would take a whole other course (or two). What I want to do is give you an intuition for how Bayesian reasoning works, show you how it differs from what we’ve been doing, and let you see that the two approaches are answering fundamentally different questions. Whether one is “better” than the other is a debate that has been raging for decades and I am not going to settle it here. But I think understanding both makes you a better consumer of statistics regardless of which camp you end up in.

18.1 The Fundamental Difference

Let’s think through how we have been thinking about hypothesis testing so far. We have some sample of data, and we ask:

Frequentist question: If the null hypothesis is true, how likely are these data?

We calculate a p-value, which is the probability of observing data as extreme or more extreme than what we got, assuming the null is true. If that probability is small enough (below our \(\alpha\) threshold), we reject the null.

Bayesian statistics flips this around entirely. Instead of asking how likely the data are given some hypothesis, a Bayesian asks:

Bayesian question: Given the data I have observed, how likely is the hypothesis?

Read those two questions again. They are not the same question. This is a point that trips up a lot of people – including a lot of professional researchers. When we calculate a p-value of 0.03, we are not saying there is a 3% chance the null hypothesis is true. We are saying there is a 3% chance of seeing data this extreme if the null were true. Those are very different statements.

The Bayesian approach gives you something that, frankly, most people think frequentist statistics gives them: a direct probability statement about the hypothesis. “There is an 85% probability that the true effect is positive.” “There is a 94% probability that Candidate A is ahead.” These are Bayesian statements. You cannot make them in the frequentist framework – even though people do it all the time without realizing they are breaking the rules.

18.2 Making use of our probability math

The center of Bayesian statistics is Bayes theorem. The good news is that we can derive this equation completely with the tools we built in the probability section at the start of the semester.

The conditional probability of A given B is equal to:

\[P(A|B) = \frac{P(A \& B)}{P(B)}\]

The probability of A and B both occurring divided by the probability of B occurring alone.

Likewise we could write out:

\[ \begin{aligned} P(B|A) &= \frac{P(B \& A)}{P(A)}\\ P(A\&B) &= P(B|A)*P(A) \end{aligned} \]

Subbing in what we just showed for \(P(A\&B)\) we get Bayes’ Theorem:

\[P(A|B) = \frac{P(B|A) \times P(A)}{P(B)}\]

That’s great that the math gives us that answer, but why is this an important formula?

Here is the most basic sort of problem we can think about using Bayes theory:

Let’s say you take a medical test for a disease and get a positive result. Given that every test has a potential for false positives (saying you have the disease when you don’t), how can we think about the probability you have the disease (A), given that you have a positive test (B).

Subbing those things in:

\[P(Disease|Pos) = \frac{P(Pos|Disease) \times P(Disease)}{P(Pos)}\]

Let’s take each of these things in turn.

What is the probability you test positive given that you have the disease, \(P(Pos|Disease)\)? This is the true positive rate of the test. People can look at data of real outcomes and ask: how often do people who actually have the disease test positive? Let’s say this is 95% or .95:

\[P(Disease|Pos) = \frac{.95 \times P(Disease)}{P(Pos)}\]

What is the overall probability of having the disease \(P(Disease)\)? This is just the overall prevalence of the disease in the population. Let’s say that is 1% or .01.

\[P(Disease|Pos) = \frac{.95 \times .01}{P(Pos)}\]

What is the probability that you test positive, \(P(Pos)\). This one is a little more tricky because you can test positive through two routes.

First, you can have the disease and test positive. 1% of people have the disease and there is 95% chance they test positive, so this probability is:

\[.01*.95= 0.0095\]

Second, you can not have the disease and test positive (a false positive). 99% of people do not have the disease and there is (let’s imagine) a 4% false positive rate for this test. This means that:

\[.99*.04 = 0.0396\].

If we add these two together we get \(.0095+.0396 = .0491\) which is the overall probability of testing positive (either because they have the disease and test positive, or they don’t have the disease and falsely test positive).

Subbing this in:

$$ \[\begin{align} P(Disease|Pos) &= \frac{.95 \times .01}{.0491}\\ P(Disease|Pos) &= .1935\\ \end{align}\] $$

After one positive there is a 19% chance that you actually have the disease.

This seems low (I mean, it’s right, but let’s think about why). Thinking through the numerator and denominator we can see why it’s low. The probability would go up if the numerator was higher. The numerator could get higher if the true positive rate of the test was higher (although it’s already pretty high), or the overall probability of having the disease was higher. The probability would go up if the denominator was lower, which we can primarily think of occurring through a lower false-positive rate.

One of the key things about Bayes theory is that it is iterative. What is the probability of having the disease if you test positive a second time? When we are going in naively our prior belief about whether you have the disease was just the population proportion, 1%. But given this test we now think that there is a 19% chance you have the disease.

We can think about this like having entered a new population where the base rate of having the disease is 19.35% instead of the former 1%. Everywhere we used that 1% before we can sub in this new prior belief.

So \(P(Disease)=.1935\), \(P(Pos|Disease)=.95\) (Unchanged), The probability of having the disease and testing positive is now \(.1935*.95=.1838\) and the probability of not having the disease and testing positive is now \(.8065*.04=.0323\), and therefore the total probability of testing positive for people with one positive test is \(.1838+.0323=.2161\).

Plugging this all into Bayes:

\[\begin{align} P(Disease|Pos_2) &= \frac{.95 \times .1935}{.2161}\\ P(Disease|Pos_2) &= .851\\ \end{align}\]

A dramatically higher probability of actually having the disease!

18.2.1 Actually Two Hypothesis

One important step in seeing how we move forward is thinking about the probability that we calculated above:

\[ \begin{align} P(Disease|Pos) &= \frac{P(Pos|Disease) \times P(Disease)}{P(Pos)}\\ P(Disease|Pos) &= \frac{.95 * .01}{(.95 * .01) + (.99*.04)}\\ P(Disease|Pos) &=.1935 \end{align} \]

The numerator is the true positive rate (how likely you are to test positive conditional on having the disease) multiplied by the probability of having the disease. The denominator represents the different ways you can test positive: (1) having the disease and testing positive (true positive); (2) not having the disease and testing positive (false positive).

The question I want hit you with is: what is \(1-P(Disease|Pos) = .8065\)? What does this complement represent?

It represents, \(P(NoDisease|Pos)\). We can see this by putting this same thing through Bayes:

\[P(NoDisease|Pos) = \frac{P(Pos|NoDisease) \times P(NoDisease)}{P(Pos)}\]

What is: \(P(Pos|NoDisease)\), the probability of testing positive if you do not have the disease? That’s just the false positive rate again, \(.04\). And what is \(P(NoDisease)\)? That’s the complement of the proportion having the disease, so .99.

The denominator is exactly the same as what we have above (a really critical part of what’s to come!).

So calculating this out:

\[ \begin{align} P(NoDisease|Pos) &= \frac{P(Pos|NoDisease) \times P(NoDisease)}{P(Pos)}\\ P(NoDisease|Pos) &= \frac{.04*.99}{(.95 * .01) + (.99*.04)}\\ P(NoDisease|Pos) &= \frac{.0396}{.0491}\\ P(NoDisease|Pos) &= .8065 \end{align} \]

Great, so the complement of our initial probability \(P(Disease|Pos)\) is just \(P(NoDisease|Pos)\). This is really important for two reasons.

First, when we are doing Bayesian style analysis, we cannot just test a single hypothesis. We have to have (at least) two hypotheses because we have to have things that sum to 1. Down below, we will see that we cannot just calculate \(p(FairCoin|7in10Heads)\), for example.

Second (and this is the most important bit), the denominator for the two hypotheses are identical, and in fact, are just the two numerators summed together. This is the mechanism that’s going to allow us to scale up to many hypotheses in what we do next.

18.2.2 Generalizing

The components of Bayes theory all have names that you will commonly see:

• \(P(A|B), P(Hypothesis|Data)\) is the posterior. Above this was the probability of having (or not) the disease given the data point of a positive test.
• \(P(B|A), P(Data|Hypothesis)\) is the likelihood. This is actually what we calculate with frequentist statistics. Above this was the probability of testing positive given that you have (or don’t have) the disease.
• \(P(A), P(Hypothesis)\) is the prior. This is our prior belief of the probability of the hypothesis. Above this was the probability of having the disease (or not) in the population.
• \(P(B), P(Data)\) is the marginal likelihood. This is the overall probability of the data across all hypotheses. But remember: this is just the sum of all of the likelihoods multipled by all the priors across all hypotheses.

Because for each hypothesis the Bayes equation divides by the same number \(P(B)\) this means that all it is doing is rescaling all of the numerators so that they sum to 1. For this reason, you will sometimes see Bayes equation written as:

\[ posterior \propto prior*likelihood \]

The posterior of our beliefs across the hypothesis is proportional to the prior times the likelihood.

18.3 Moving towards data: another coin flip example

Here’s our old problem: we flip a coin 10 times and get 7 heads. What can we conclude?

The way we deal with this in a frequentist framework should be old hat to you now:

We posit something true about the world (this is a fair coin), and then determine the likelihood of getting the data we got if that were true.

x <- seq(0,10)
plot(x, dbinom(x,10, prob=.5))
abline(v=7, lty=2)

And make some sort of pronouncement of what the probability is of observing this data given this hypothesis.

Again: Bayesian statistics is going to flip this on it’s head. Now we are going to answer \(P(Hypothesis|Data)\): what is the probability of our initial hypothesis given the data. Above we determined: what is the probability that I have the disease (hypothesis) given a positive test. Now we are going to answer: what is the probability this is a certain type of coin given the data that I have.

Writing out the Bayes formula for this:

\[P(\text{Hypothesis} | \text{Data}) = \frac{P(\text{Data} | \text{Hypothesis}) \times P(\text{Hypothesis})}{P(\text{Data})}\]

As we discussed above, this is not as easy as doing \(P(FairCoin|7in10Flips)\). We actually need to choose the alternative hypotheses to test across.

Let’s do three hypotheses: the coin comes up heads 25% of the time, the coin comes up heads 50% of the time, and the coin comes up heads 75% of the time. Let’s further say that we think each of these outcomes is equally likely. If we were to graphically display that prior belief about those coins it would look like:

pi <- c(.25,.50,.75)
prior <- rep(1/3,3)

plot(pi, prior, main="Prior Belief", pch=16)

We are going to calculate with Bayes rule for each of these equations, but remember that we showed above that the denominator for all three equations are identical, and in fact are just the sum of the three numerators. Because of this, we are just going to focus on the three numerators.

So numerator one is going to be:

\[ P(\pi=.25|7in10Heads) = P(7in10Heads|\pi=.25) *P(\pi=.25) \]

What is \(P(7in10Heads|\pi=.25)\) the probability of getting 7 in 10 heads for a coin that comes up heads 25% of the time? This should be a very familiar question: it’s just the binomial.

dbinom(7,10,.25)

[1] 0.003089905

A very small number. It’s extremely unlikely to get this.

What’s \(P(\pi=.25)\) the probability that the coin comes up heads 25% of the time? That’s just our prior belief about the prevalence of these coins, which we set above as \(1/3\).

So to calculate this numerator:

dbinom(7,10,.25)*(1/3)

[1] 0.001029968

Or to calculate all three at once:

post <- dbinom(7,10,pi)*prior
plot(pi,post, pch=16, main="Posterior, Un-normalized")

Now these numbers are not giving us real probabilities yet because we haven’t divided by the denominator \(P(Data)\), the overall probability of getting 7 in 10 heads across all three hypotheses. Remember that \(posterior \propto prior*likelihood\). Our actual posterior will be proportional to these numbers because they will all be divided by the same number. So before we even get there, we can say that the hypothesis with the highest probability is that this coin comes up heads 75% of the time.

To actually get those probabilities we divide each by the sum of the numerators:

post.norm <- post/sum(post)
plot(pi, post.norm, pch=16, main="Posterior, Normalized")

post.norm

[1] 0.008338481 0.316244595 0.675416924

Our posterior beliefs are that there is a less than 1% chance this coin comes up heads 25% of the time, around a 30% belief it comes up heads 50% of the time, and a 67% chance that it comes up heads 75% of the time.

18.3.1 Changing our priors

It should be very clear that Bayesian analysis is sensitive to your priors. We (somewhat ridiculously) posited a prior distribution that all three coins were equally likely. But that’s not true. Unfair coins are pretty rare. So what happens if we look at the same data, and same hypotheses, but instead have a prior that there is a 95% chance this is a fair coin, with just a 2.5% chance for each of the unfair coins?

#Set priors
prior <- c(.025,.975,.025)
plot(pi, prior, pch=16, main="Prior")

#Calculate likelihoods*priors
post <- dbinom(7,10,pi)*prior

#Normalise by dividing by the sum
post.norm <- post/sum(post)

#Display posterior
plot(pi, prior, pch=16, main="Posterior", col="gray80")
points(pi, post.norm, pch=16)

post.norm

[1] 0.0006405694 0.9474733096 0.0518861210

Now this data has barely moved us off our (much stronger) priors. Because we went in with the belief that unfair coins are extremely uncommon, having an observation of 7/10 heads isn’t very persuasive evidence that we are dealing with an unfair coin.

However, if we got more data that would begin to change.

What if we kept flipping and we got 70 in 100 heads?

#Calculate likelihoods*priors (This is where the change in the data is)
post <- dbinom(70,100,pi)*prior

#Normalise by dividing by the sum
post.norm <- post/sum(post)

#Display posterior
plot(pi, prior, pch=16, main="Posterior", col="gray80")
points(pi, post.norm, pch=16)

post.norm

[1] 8.065957e-20 1.936790e-02 9.806321e-01

That data overwhelmed our strong priors and we end up a high degree of certainty that this is an unfair coin.

18.3.2 Scaling this up

Because of the way that we set this up this is extremely easy to scale up to nearly-continuous hypotheses. What if we want to test the probability of many possible \(pi\) from 0 to 1?

pi <- seq(0,1,.001)

Now we have 1001 priors from 0 to 1.

Again, let’s assign these equal probabilities, which is just 1 divided by the number of hypotheses:

prior <- rep(1/1001,1001)
plot(pi, prior, type="l", main="Prior")

And then let’s run this through Bayes formula with the data that we have 7 heads in 10 flips:

post <- dbinom(7,10,pi)*prior
post.norm <- post/sum(post)

plot(pi, post.norm, type="l", main="Posterior")

Now we have a posterior distribution across our possible hypotheses.

When we were doing frequentist statistics we could only accept or reject the null hypothesis. The nice thing here is that we now have a distribution from which we can make actual probability claims on.

So for example: given these data, what is the probability that this coin comes up heads more than 50% of the time? That’s just adding up the probability mass in the curve to the right of .5:

sum(post.norm[pi>.5])

[1] 0.8860734

There is an 88.6% chance that this coin comes up heads more than 50% of the time given our data.

What about a better prior that more accurately shows that unfair coins are rare. For these sort of probability questions we can use the beta distribution, which is a random variable that assigns density to values between 0 and 1. We are going to generate it and divide by the sum so that the numbers add up to 1. (Don’t worry too much about how the beta distribution works, just think about this more continuously describing the prior distribution of possibilities.)

prior <- dbeta(pi, 100,100)
prior <- prior/sum(prior)
plot(pi, prior, type="l", main="Strong Prior")

Applying Bayes rule and our observation of 7 heads in 10 flips:

post <- dbinom(7,10,pi)*prior
post.norm <- post/sum(post)
plot(pi, prior, type="l", col="gray80", main="Posterior")
points(pi, post.norm, type="l")

With this small amount of data our posterior distribution is barely different from our prior distribution.

But if we scale up to a world where we have 70 heads in 100 flips:

post <- dbinom(70,100,pi)*prior
post.norm <- post/sum(post)
plot(pi, prior, type="l", col="gray80", main="Posterior")
points(pi, post.norm, type="l")

Out posterior changes more significantly.

Again, we can answer questions with the distribution that we could never with frequentist statistics. What’s the probability that the probability of the fair coin is between 55% and 60%?

sum(post.norm[pi<.6]) - sum(post.norm[pi<.55])

[1] 0.6016363

About 60%.

What are the two points of this distribution where 95% of the probability mass lies?

#Calculate cumulative probability
cumu.post <- cumsum(post.norm)
#Find the boundaries
pi[min(which(cumu.post>=.975))]

[1] 0.622

pi[max(which(cumu.post<=.025))]

[1] 0.509

The 95% credible interval (the Bayesian equivalent of a confidence interval) runs from 50.9% to 62.2%.

18.4 Real world application

How could we apply this to the real world?

The easiest application is to a scenario where we are interrogating a probability, because that maps cleanly on to the coin flip example that we have been working with.

Let’s say that we have the following poll of 1000 people on whether they support a Democratic (1) or Republican (0) candidate

x <- rbinom(1000,1,.52)
mean(x)

[1] 0.487

where 51.4% support the Democratic candidate.

In hypothesis testing world we do a t.test with a null hypothesis that

t.test(x, mu=.50)

    One Sample t-test

data:  x
t = -0.82206, df = 999, p-value = 0.4112
alternative hypothesis: true mean is not equal to 0.5
95 percent confidence interval:
 0.4559676 0.5180324
sample estimates:
mean of x 
    0.487 

And we see that the probability of observing these data if the truth was that the race is tied is .376 or around 38%. As such we almost certainly would fail to reject the null hypothesis.

How would we deal with this in a Bayesian fashion?

Let’s first define a prior. Maybe this race is a toss-up, so we can center our prior on 50%, and have it cover most probabilities between 40 and 60%:

dem.support <- seq(0,1,.01)
prior <- dbeta(dem.support, 50,50)
prior <- prior/sum(prior)

plot(dem.support, prior, type="l", main="Prior")

And then we consider Bayes formula for all our hypotheses: that the true dem support is 0%, 0.01%, 0.02%, filtered through our data where 514 out of 1000 people support the dem candidate:

post <- dbinom(514,1000,dem.support)*prior
post.norm <- post/sum(post)

plot(dem.support, prior, type="l", main="Prior", col="gray80")
points(dem.support, post.norm, type="l")
abline(v=.5, lty=2)

Unlike frequentist statistics where we can only say “We cannot reject the idea that the race is tied”, here we can give way more specific guidance.

What is the probability that the democratic candidate will win comfortably, with over 55% of the vote:

sum(post.norm[dem.support>=.55])

[1] 0.01442396

Not large. Just 1.4% probability of that occuring.

Or calculate the credible interval:

#Calculate cumulative probability
cumu.post <- cumsum(post.norm)
#Find the boundaries
dem.support[min(which(cumu.post>=.975))]

[1] 0.54

dem.support[max(which(cumu.post<=.025))]

[1] 0.47

What gets more helpful with Bayesian statistics is that we can actually have more informative priors.

What if there has already been some polling in this race: there have been three polls all with n=1000 with the results 51%, 55%, and 49.5%.

You don’t need to know this part, but we could combine these into a prior by inputting the number of successes as the first paramater and the number of failures as the second paramater:

prior <- dbeta(dem.support, 510+550+495, 490+450+505)
prior <- prior/sum(prior)
plot(dem.support, prior, type="l", main="Prior")
abline(v=.5,lty=2)

Now with this more informative prior we actually get a more certain answer from adding in our poll:

post <- dbinom(514,1000,dem.support)*prior
post.norm <- post/sum(post)
plot(dem.support, prior, type="l", main="Posterior", col="gray80")
points(dem.support, post.norm, type="l")
abline(v=.5,lty=2)

Now our credible interval is:

#Calculate cumulative probability
cumu.post <- cumsum(post.norm)
#Find the boundaries
dem.support[min(which(cumu.post>=.975))]

[1] 0.53

dem.support[max(which(cumu.post<=.025))]

[1] 0.49

Which is narrower despite having the same data!

18.5 Getting the posterior a different way: sampling

Every posterior we have built has come the same way. We lay out a grid of candidate values, we compute the likelihood times the prior at each one, and we divide by the sum. Enumerate every hypothesis, weigh each one, normalize. For a single parameter between 0 and 1 that is easy, and it is exactly right.

I am now going to show you a second way to get the very same posterior, and I want to be honest up front: for a problem this simple it is completely unnecessary. The grid is fine. But this second method is the only one that survives when we move to realistic models, so I want to introduce it here, in the simplest possible setting, where you can check its answer against the grid you already trust.

The idea is to sample from the posterior instead of enumerating it. Rather than visiting every candidate value of \(p\) and recording its height, we take a random walk across the values of \(p\), a walk that lingers where the posterior is tall and hurries through where it is short. If we design the walk correctly, then the values it visits, collected up, are the posterior. Where the posterior is twice as tall, the walk spends twice as long, so twice as many collected values land there.

The recipe for the walk is called the Metropolis algorithm, and it is startlingly simple. At each step:

1. Stand at your current value of \(p\).
2. Propose a nearby value by taking a small random step.
3. Compare the posterior height at the proposed value to your current one. If the proposal is taller, move to it. If it is shorter, move to it only sometimes, with probability equal to how much shorter it is (a proposal half as tall gets accepted half the time).
4. Write down wherever you now stand. Repeat.

Step 3 is the whole engine. Always taking uphill steps pulls the walk toward the tall part of the posterior. Sometimes taking downhill steps stops it from getting stuck at the very peak and lets it explore the whole shape, spending time in each region in proportion to that region’s height.

Let’s build it for our poll. Same data (514 of 1000 supporting the Democrat), same toss-up prior. First we need a function that returns the posterior height at any candidate value of \(p\). That is just likelihood times prior, the same thing we have been computing all along:

post.poll <- function(p) {
  if (p <= 0 || p >= 1) return(0)  #p has to be a probability
  dbinom(514, 1000, p) * dbeta(p, 50, 50)
}

Now the walk. This is the entire sampler:

set.seed(19104)
n.iter <- 20000
p.draws <- rep(NA, n.iter)

p.cur    <- 0.5                    #start somewhere
post.cur <- post.poll(p.cur)
step     <- 0.02                   #how big a step we propose

for (i in 1:n.iter) {
  p.prop    <- p.cur + rnorm(1, 0, step)   #propose a nearby value
  post.prop <- post.poll(p.prop)
  #move there with probability equal to the ratio of the two heights
  if (runif(1) < post.prop / post.cur) {
    p.cur    <- p.prop
    post.cur <- post.prop
  }
  p.draws[i] <- p.cur                       #write down where we are
}

#throw away the early "burn-in" steps that depend on where we started
p.draws <- p.draws[2001:n.iter]

We now have 18,000 values of \(p\). Did the walk actually reproduce the posterior? Let’s put a histogram of the sampled values against the grid posterior (I’ll use a fine grid here as the benchmark):

#Fine-grid posterior as the benchmark
p.grid <- seq(0, 1, 0.001)
prior  <- dbeta(p.grid, 50, 50); prior <- prior / sum(prior)
grid.post <- dbinom(514, 1000, p.grid) * prior
grid.post <- grid.post / sum(grid.post)

hist(p.draws, breaks = 40, freq = FALSE, xlim = c(.45, .58),
     main = "MCMC draws vs. the grid posterior",
     xlab = "p (Democratic support)")
#overlay the grid posterior (rescaled to a density for comparison)
points(p.grid, grid.post / mean(diff(p.grid)), type = "l", col = "firebrick", lwd = 2)

The histogram of sampled values traces out the same curve the grid gave us. We never computed the posterior at every point. We wandered around following the accept-reject rule, and the values we happened to visit piled up into the exact shape we were after.

So what did we sample? A big vector of plausible values of \(p\). Each draw is one level of Democratic support that the data and prior find believable, and the believable values show up in the vector more often than the iffy ones. To answer any question we summarize that vector, the same way we summarized the grid, and the same way we summarized bootstrap resamples back in the Power chapter:

mean(p.draws)                     #posterior mean

[1] 0.5126977

quantile(p.draws, c(.025, .975))  #95% credible interval

     2.5%     97.5% 
0.4836261 0.5422909 

mean(p.draws > 0.5)               #P(Democrat ahead)

[1] 0.8002778

A posterior mean of about 0.513, a credible interval of roughly 0.48 to 0.54, and about an 80% probability the Democrat is ahead. These are the same numbers the grid gave us, because the two methods compute the same thing, one by enumeration and one by sampling.

For a single parameter this really was pointless: the grid was easier and gave an identical answer. The reason it matters comes now, when the hypothesis stops being a single number.

18.6 More than one parameter: regression

Now the real payoff, and the thing that makes any of this worth the trouble. Everything above had a single unknown: the coin’s \(\pi\), the poll’s \(p\). But almost every model worth fitting has several. Take the simplest regression, fitting a line \(y = \alpha + \beta x\). Here a hypothesis is not one number, it is a pair: a candidate intercept \(\alpha\) and a candidate slope \(\beta\). Every pair is one particular line drawn through the data.

Nothing about the recipe changes. We still posit a possibility, score how well it fits the data, and let the walker shift around toward the possibilities that fit best. The only difference is that a possibility is now a pair of numbers instead of one, so the walker roams a two-dimensional plane of candidate \((\alpha, \beta)\) pairs instead of a one-dimensional line of candidate \(p\) values.

Let’s make dead-simple data where we know the answer, an intercept of 2 and a slope of 0.5:

set.seed(19104)
n <- 50
x <- runif(n, 0, 10)
y <- 2 + 0.5 * x + rnorm(n, mean = 0, sd = 1)
plot(x, y, pch = 16, main = "Our data")

We still need a way to score a candidate line: given this intercept and slope, how probable is the data? Same idea as always: a line through the middle of the points fits well and scores high, a line that misses them scores low. One small difference from the poll is worth a word. There the score came from a single dbinom, a tidy number. A line’s score is the probability of all 50 points at once, which is the product of 50 little densities, and that product is a very small number. So we work with its log, which keeps it a readable size. The score below is that log: it comes out negative, and a bigger (less negative) score means a better-fitting line.

#Log-score of a candidate line: bigger (less negative) = better fit
log.post <- function(alpha, beta) {
  sum(dnorm(y, mean = alpha + beta * x, sd = 1, log = TRUE)) +
  dnorm(alpha, 0, 10, log = TRUE) + dnorm(beta, 0, 10, log = TRUE)
}

You can watch the score sort good lines from bad. Here are three guesses; the green line runs through the data and scores highest (least negative), the red and blue miss and score far lower:

plot(x, y, pch = 16, main = "Three candidate lines")
abline(a = 0, b = 1,   col = "firebrick")   #(0, 1)
abline(a = 4, b = 0,   col = "dodgerblue")  #(4, 0)
abline(a = 2, b = 0.5, col = "forestgreen") #(2, 0.5)

log.post(0, 1)    #red

[1] -119.7835

log.post(4, 0)    #blue

[1] -140.0518

log.post(2, 0.5)  #green

[1] -74.44312

Now we turn the walker loose. It is the identical sampler from the poll, with one change: each step proposes a small move in the intercept and the slope at the same time, so it wanders the two-dimensional plane instead of the line. The accept-reject rule is the same “move with probability equal to the ratio of the two heights” as before; because our scores are now logs, that ratio is exp of their difference.

set.seed(19104)
n.iter <- 20000
draws <- matrix(NA, nrow = n.iter, ncol = 2)
colnames(draws) <- c("alpha", "beta")

alpha.cur <- 0; beta.cur <- 0
lp.cur <- log.post(alpha.cur, beta.cur)
step <- 0.15

for (i in 1:n.iter) {
  alpha.prop <- alpha.cur + rnorm(1, 0, step)
  beta.prop  <- beta.cur  + rnorm(1, 0, step)
  lp.prop    <- log.post(alpha.prop, beta.prop)
  #move with probability equal to the ratio of the two heights
  if (runif(1) < exp(lp.prop - lp.cur)) {
    alpha.cur <- alpha.prop; beta.cur <- beta.prop; lp.cur <- lp.prop
  }
  draws[i, ] <- c(alpha.cur, beta.cur)
}
draws <- draws[2001:n.iter, ]

Here is where the walker spent its time, plotted in the plane of intercept-slope pairs:

plot(draws[, "alpha"], draws[, "beta"], pch = 16, cex = .3,
     col = rgb(0, 0, 0, .1),
     xlab = "Intercept (alpha)", ylab = "Slope (beta)",
     main = "Where the walker went")

That cloud is the posterior over \((\alpha, \beta)\). It leans on a diagonal because a steeper slope has to be paired with a lower intercept to keep the line running through the data, so the two are not independent.

18.6.1 The posterior is a cloud of lines

What did we actually sample? On the poll each draw was a single number and the whole thing was a vector. Here each draw is a pair, an intercept and a slope, which is to say each draw is a whole line. So the posterior is a distribution over lines. We can see that literally by drawing a few hundred of the sampled lines faintly over the data:

plot(x, y, pch = 16, main = "The posterior is a cloud of lines")
for (i in seq(1, nrow(draws), length.out = 200)) {
  abline(a = draws[i, "alpha"], b = draws[i, "beta"], col = rgb(.2, .2, .8, .05))
}
points(x, y, pch = 16)

Where the lines bunch tightly the posterior is confident; where they fan out it is uncertain. That fan is the Bayesian answer to “what is the relationship between x and y”: not one best line, but the full set of lines the data find plausible.

From here it is all summarizing columns of draws, the same as always. The intercept draws are the posterior for \(\alpha\), the slope draws the posterior for \(\beta\):

colMeans(draws)                            #posterior means

    alpha      beta 
1.8811654 0.5418772 

quantile(draws[, "beta"], c(.025, .975))   #95% credible interval for the slope

     2.5%     97.5% 
0.4414345 0.6426614 

mean(draws[, "beta"] > 0)                  #P(slope is positive)

[1] 1

The slope posterior is centered near 0.54, with a credible interval of about 0.44 to 0.64 and essentially a 100% probability the slope is positive. Least squares gives the same interval:

confint(lm(y ~ x))

                2.5 %   97.5 %
(Intercept) 1.3682705 2.432993
x           0.4453066 0.634123

Same answer, different meaning, exactly as in the coin example.

18.6.2 Adding more parameters changes nothing

Stop and notice how little had to change to go from one parameter to two. On the poll we proposed a single new number each step; in the regression we proposed two. That was the whole difference. If we also wanted to estimate \(\sigma\) we would propose three numbers each step. A regression with ten predictors has twelve unknowns (ten slopes, an intercept, and \(\sigma\)), so the walk would propose twelve numbers each step. Nothing about the accept-reject rule changes. The walker does not know or care how many dimensions it is walking in, because at every step it only ever compares two heights, and a height is a single number no matter how many dimensions produced it.

This is the whole reason we sample instead of enumerate. A grid explodes as you add parameters: a hundred candidate values per unknown is ten thousand points in two dimensions, a hundred trillion in seven, and hopeless past that. The random walk shrugs that off. It never tries to hold the whole posterior in view the way a grid does, it just takes one step at a time, comparing where it is to where it might go. That is what lets real Bayesian software fit models with dozens or hundreds of parameters that no grid could ever touch.

18.6.3 This is what the software does

What you just wrote, in about fifteen lines, is a baby version of what that software does. When people fit Bayesian models in practice they use tools like Stan, brms, or PyMC. Those tools use much cleverer walkers than our little random step, and they sample \(\sigma\) and every other unknown along with the slopes. But the thing they hand back to you is the object you just built: a big table of draws, one row per plausible set of parameter values, that you summarize with means and quantiles.

If that sounds familiar, it should. Back in the bootstrap chapter we made a sampling distribution by resampling our data thousands of times and collecting a statistic from each resample, ending up with a big vector of numbers we summarized with quantile() and mean(). MCMC does the mirror-image thing: it collects thousands of plausible parameter values and hands you a big table you summarize the same way. In both cases the hard part is producing the pile of draws. Once you have it, reading it is something you already know how to do.

18.7 Where does this go?

END OF WHAT IVE WRITTEN

In the context of statistical inference, we replace “A” with a hypothesis (or parameter value) and “B” with our observed data:

\[P(\text{Hypothesis} | \text{Data}) = \frac{P(\text{Data} | \text{Hypothesis}) \times P(\text{Hypothesis})}{P(\text{Data})}\]

Each piece of this equation has a name:

• \(P(\text{Hypothesis} | \text{Data})\): the posterior. This is what we want. It is our updated belief about the hypothesis after seeing the data.
• \(P(\text{Data} | \text{Hypothesis})\): the likelihood. This is how probable the data are under a specific hypothesis. This is actually the thing that frequentist statistics focuses on!
• \(P(\text{Hypothesis})\): the prior. This is what we believed about the hypothesis before seeing the data.
• \(P(\text{Data})\): the marginal likelihood (or evidence). This is the overall probability of the data across all possible hypotheses. It serves as a normalizing constant.

The Bayesian recipe is simple in words:

\[\text{Posterior} \propto \text{Likelihood} \times \text{Prior}\]

Your updated belief is proportional to what the data tell you multiplied by what you believed going in. That is the whole game.

18.8 A Simple Example: Is This Coin Fair?

Let’s make this concrete. Suppose I hand you a coin and you want to figure out whether it is fair. You flip it 10 times and get 7 heads.

18.8.1 The Frequentist Approach

We have done exactly this kind of problem before. We set up a hypothesis test:

• \(H_0\): \(\pi = 0.5\) (the coin is fair)
• \(H_A\): \(\pi \neq 0.5\) (the coin is not fair)

We calculate the probability of getting 7 or more heads (or 3 or fewer, for the two-tailed test) under the null:

set.seed(19104)

#Two-tailed p-value for 7 heads in 10 flips, fair coin
p.value <- pbinom(2, 10, 0.5) + (1 - pbinom(6, 10, 0.5))
p.value

[1] 0.2265625

With a p-value of about 0.34, we fail to reject the null. We cannot rule out that this coin is fair. Note what we are not saying: we are not saying the coin is fair. We are not saying there is a 66% chance the coin is fair. We are saying the data are not surprising enough to reject the hypothesis of fairness. That is a pretty unsatisfying answer, but it is the only answer frequentism gives us.

18.8.2 The Bayesian Approach

A Bayesian thinks about this problem differently. Instead of testing a single null hypothesis, a Bayesian considers all possible values of \(\pi\) (the true probability of heads) and asks: given my data, which values of \(\pi\) are most plausible?

To do this, we need three things:

1. A prior distribution: What did we believe about \(\pi\) before flipping the coin?
2. A likelihood: How probable are 7 heads in 10 flips for each possible value of \(\pi\)?
3. Apply Bayes’ theorem to get the posterior distribution.

Let’s start simple. Suppose we have no idea what \(\pi\) is – we think every value between 0 and 1 is equally likely. That is a uniform prior:

#Define a grid of possible pi values
pi.vals <- seq(0, 1, 0.001)

#Uniform prior: every value equally likely
prior <- rep(1, length(pi.vals))
#Normalize so it sums to 1
prior <- prior / sum(prior)

plot(pi.vals, prior, type="l",
     xlab="Pi (Probability of Heads)", ylab="Prior Belief",
     main="Uniform Prior", ylim=c(0, max(prior)*3))

Now we need the likelihood. For each possible value of \(\pi\), what is the probability of observing 7 heads in 10 flips? We already know this – it is the binomial probability:

#Likelihood: P(7 heads in 10 flips | pi) for each possible pi
likelihood <- dbinom(7, 10, pi.vals)

plot(pi.vals, likelihood, type="l",
     xlab="Pi (Probability of Heads)", ylab="Likelihood",
     main="Likelihood of 7 Heads in 10 Flips")

This should look familiar – we have been working with the binomial distribution all semester. The likelihood is maximized at \(\pi = 0.7\), which makes sense: the value of \(\pi\) that makes our data most likely is the one that exactly matches our sample proportion.

Now we apply Bayes’ theorem. The posterior is proportional to the prior times the likelihood:

#Posterior = Prior x Likelihood (then normalize)
posterior <- prior * likelihood
posterior <- posterior / sum(posterior)

plot(pi.vals, posterior, type="l",
     xlab="Pi (Probability of Heads)", ylab="Posterior Belief",
     main="Posterior Distribution After 7/10 Heads")
abline(v=0.5, lty=2, col="gray50")
abline(v=0.7, lty=2, col="firebrick")
legend("topleft", c("Pi = 0.5", "Pi = 0.7"),
       lty=2, col=c("gray50","firebrick"))

This is the posterior distribution. It tells us, given our data and our prior beliefs, how plausible each value of \(\pi\) is. The most probable value is 0.7 (our sample proportion), but there is still quite a bit of probability mass around 0.5. With only 10 flips, we are pretty uncertain.

Now here is the cool part. Unlike the frequentist framework, we can make direct probability statements:

#What is the probability that the coin is biased towards heads?
prob.biased <- sum(posterior[pi.vals > 0.5])
prob.biased

[1] 0.8860734

We can say: “Given our data, there is about a 89% probability that the coin is biased towards heads.” That is a Bayesian statement. It directly answers the question most people actually care about. Compare that to the frequentist answer: “We fail to reject the null hypothesis that \(\pi = 0.5\) at the 0.05 significance level.” Both are correct, but they are answering different questions.

18.9 The Prior Matters (and That’s the Controversial Part)

You might have noticed something. The Bayesian answer depends on the prior. We used a uniform prior above – we assumed every value of \(\pi\) was equally likely before seeing the data. But what if we had different prior beliefs?

Let’s say I told you this coin came straight from the U.S. Mint and I just pulled it out of a fresh roll. You would probably believe, quite strongly, that it is close to fair. We can encode that belief as a prior that is concentrated around 0.5.

For the mathematically curious: we are going to use a Beta distribution as our prior. The Beta distribution is defined on the interval [0, 1], which makes it perfect for modeling beliefs about probabilities. It has two parameters, \(a\) and \(b\), and its mean is \(\frac{a}{a+b}\). When \(a = b\), the distribution is symmetric around 0.5. Larger values of \(a\) and \(b\) make the distribution more concentrated.

#Strong prior: we think the coin is probably fair
strong.prior <- dbeta(pi.vals, 20, 20)
strong.prior <- strong.prior / sum(strong.prior)

#Weak prior (uniform)
weak.prior <- rep(1, length(pi.vals))
weak.prior <- weak.prior / sum(weak.prior)

plot(pi.vals, strong.prior, type="l", col="dodgerblue", lwd=2,
     xlab="Pi", ylab="Prior Belief", main="Two Different Priors")
points(pi.vals, weak.prior, type="l", col="firebrick", lwd=2)
legend("topright", c("Strong prior (probably fair)", "Weak prior (uniform)"),
       col=c("dodgerblue","firebrick"), lwd=2)

Now let’s see what happens when we update both priors with the same data (7 heads in 10 flips):

#Posterior with strong prior
posterior.strong <- strong.prior * likelihood
posterior.strong <- posterior.strong / sum(posterior.strong)

#Posterior with weak prior
posterior.weak <- weak.prior * likelihood
posterior.weak <- posterior.weak / sum(posterior.weak)

plot(pi.vals, posterior.weak, type="l", col="firebrick", lwd=2,
     xlab="Pi", ylab="Posterior Belief",
     main="Posteriors Under Different Priors",
     ylim=c(0, max(c(posterior.weak, posterior.strong))))
points(pi.vals, posterior.strong, type="l", col="dodgerblue", lwd=2)
abline(v=0.7, lty=3)
abline(v=0.5, lty=3)
legend("topleft", c("Weak prior (uniform)", "Strong prior (probably fair)"),
       col=c("firebrick","dodgerblue"), lwd=2)

Look at what happened. With the weak (uniform) prior, the posterior is centered near 0.7 – our data dominate. With the strong prior, the posterior is pulled towards 0.5 – our prior belief that the coin is fair is pulling the estimate towards 0.5, and 10 flips is not enough data to overcome that belief.

This is the feature – or bug, depending on who you ask – of Bayesian statistics. The prior matters. Different priors lead to different conclusions. Frequentists see this as a problem: statistics should be objective, and the prior introduces subjectivity. Bayesians see this as a feature: of course your prior knowledge should matter! If you know the coin came from the Mint, that information is relevant and you should use it.

Here is the important thing though: as data accumulate, the prior matters less and less. Let’s see this in action. Instead of 10 flips, let’s say we flip the coin 1000 times and get 700 heads:

#Likelihood with more data
likelihood.big <- dbinom(700, 1000, pi.vals)

#Posterior with strong prior and lots of data
posterior.strong.big <- strong.prior * likelihood.big
posterior.strong.big <- posterior.strong.big / sum(posterior.strong.big)

#Posterior with weak prior and lots of data
posterior.weak.big <- weak.prior * likelihood.big
posterior.weak.big <- posterior.weak.big / sum(posterior.weak.big)

plot(pi.vals, posterior.weak.big, type="l", col="firebrick", lwd=2,
     xlab="Pi", ylab="Posterior Belief",
     main="Posteriors Converge With More Data (n=1000)",
     ylim=c(0, max(c(posterior.weak.big, posterior.strong.big))))
points(pi.vals, posterior.strong.big, type="l", col="dodgerblue", lwd=2)
legend("topleft", c("Weak prior", "Strong prior"),
       col=c("firebrick","dodgerblue"), lwd=2)

With enough data, the two posteriors are essentially identical. The data overwhelm the prior. This is an incredibly important property: reasonable people with different priors will converge to the same conclusion given enough evidence. The prior only matters when data are sparse.

18.10 Bayesian Updating: Learning Sequentially

One of the most elegant features of Bayesian statistics is that updating is sequential. Today’s posterior becomes tomorrow’s prior. Let’s watch a Bayesian learner update their beliefs about a coin one flip at a time.

set.seed(19104)

#True probability of heads
true.pi <- 0.6

#Generate 100 flips
flips <- rbinom(100, 1, true.pi)

#Start with a uniform prior (Beta(1,1))
a <- 1
b <- 1

#Track the posterior mean after each flip
posterior.means <- rep(NA, 100)

par(mfrow=c(2,3))
for(i in 1:100){
  if(flips[i] == 1){
    a <- a + 1  #Heads: increase a
  } else {
    b <- b + 1  #Tails: increase b
  }

  posterior.means[i] <- a / (a + b)

  #Plot the posterior at selected points
  if(i %in% c(1, 5, 10, 25, 50, 100)){
    curve(dbeta(x, a, b), from=0, to=1,
          xlab="Pi", ylab="Density",
          main=paste("After", i, "flips"))
    abline(v=0.6, lty=2, col="firebrick")
  }
}

par(mfrow=c(1,1))

After 1 flip, we know almost nothing. After 5 flips, we have a rough idea. After 25 flips, the posterior is starting to concentrate. After 100 flips, we are tightly centered around the true value of 0.6.

We can also track the posterior mean over time – watch how our best guess evolves:

plot(1:100, posterior.means, type="l",
     xlab="Number of Flips", ylab="Posterior Mean (Best Guess for Pi)",
     main="Bayesian Updating of Coin Bias Estimate")
abline(h=0.6, lty=2, col="firebrick")
legend("topright", "True Pi = 0.6", lty=2, col="firebrick")

This should remind you of something. Remember the Law of Large Numbers plot from the Sampling chapter? We showed that as the sample size grows, the sample mean converges to the population mean. The Bayesian posterior mean does the same thing – it converges to the truth. The two frameworks are approaching the same target, just from different philosophical starting points.

18.11 A Political Example: Updating Beliefs About an Election

Let’s bring this back to politics, where we have spent most of the semester. Suppose it’s September before a presidential election, and we want to estimate the Democratic candidate’s share of the two-party vote.

Before any polling, we have some prior belief. Presidential elections in the United States tend to be pretty close – they almost always fall between 45% and 55% of the two-party vote. A reasonable prior might be centered on 50% (no prior advantage to either party) with most of the mass between 45% and 55%.

set.seed(19104)

pi.vals <- seq(0, 1, 0.001)

#Prior: centered on 0.50, most mass between .45 and .55
#Beta(50, 50) has mean 0.5 and is concentrated
prior <- dbeta(pi.vals, 50, 50)
prior <- prior / sum(prior)

plot(pi.vals, prior, type="l", lwd=2,
     xlab="Democratic Two-Party Vote Share",
     ylab="Prior Belief",
     main="Prior Belief About Election Outcome",
     xlim=c(.3,.7))

Now suppose a poll of 1,000 likely voters comes in, and 53% support the Democrat. In the frequentist framework, we would construct a confidence interval around 0.53 and do a hypothesis test against the null of 0.50. Let’s first do that, since we know how:

#Frequentist approach
n.poll <- 1000
p.hat <- 0.53
se <- sqrt(p.hat * (1 - p.hat) / n.poll)

#95% confidence interval
ci.lower <- p.hat - 1.96 * se
ci.upper <- p.hat + 1.96 * se
cat("Frequentist 95% CI:", round(ci.lower, 3), "to", round(ci.upper, 3), "\n")

Frequentist 95% CI: 0.499 to 0.561 

#Hypothesis test: H0: pi = 0.5
z <- (p.hat - 0.5) / se
p.val <- 2 * pnorm(-abs(z))
cat("P-value for H0: pi = 0.5:", round(p.val, 4))

P-value for H0: pi = 0.5: 0.0573

Now the Bayesian approach. We update our prior with the poll data:

#Likelihood: 530 out of 1000 support Democrat
likelihood <- dbinom(530, 1000, pi.vals)

#Posterior
posterior <- prior * likelihood
posterior <- posterior / sum(posterior)

plot(pi.vals, prior, type="l", lwd=2, col="gray60",
     xlab="Democratic Two-Party Vote Share",
     ylab="Density",
     main="Prior, Likelihood, and Posterior",
     xlim=c(.4, .6),
     ylim=c(0, max(c(prior, posterior))*1.1))
points(pi.vals, likelihood / max(likelihood) * max(posterior),
       type="l", lwd=2, col="darkorange", lty=2)
points(pi.vals, posterior, type="l", lwd=2, col="dodgerblue")
legend("topleft", c("Prior", "Likelihood (scaled)", "Posterior"),
       col=c("gray60","darkorange","dodgerblue"), lwd=2, lty=c(1,2,1))

The posterior (blue) is a compromise between our prior (gray) and the data (orange). It is pulled slightly towards 0.50 compared to the raw poll result because our prior belief was that the election would be close. With 1,000 respondents the data are doing most of the work, but the prior still has a small effect.

Now here is the payoff. We can make direct probability statements about the election:

#Probability that the Democrat is ahead
prob.dem.ahead <- sum(posterior[pi.vals > 0.5])
cat("Probability Democrat is ahead:", round(prob.dem.ahead, 3), "\n")

Probability Democrat is ahead: 0.962 

#Bayesian "credible interval" -- central 95% of the posterior
#Find the 2.5th and 97.5th percentiles of the posterior
cumulative <- cumsum(posterior)
ci.bayes.lower <- pi.vals[min(which(cumulative >= 0.025))]
ci.bayes.upper <- pi.vals[min(which(cumulative >= 0.975))]
cat("Bayesian 95% Credible Interval:", ci.bayes.lower, "to", ci.bayes.upper)

Bayesian 95% Credible Interval: 0.498 to 0.557

That first number – the probability the Democrat is ahead – is exactly the kind of statement election forecasters make. When you see a model say “Candidate X has a 78% chance of winning,” that is a Bayesian probability. You literally cannot produce that number in the frequentist framework.

The second output is a credible interval. It looks a lot like a confidence interval, and numerically the two are often very similar. But the interpretation is completely different:

• Frequentist 95% CI: If we repeated this poll many times, 95% of the resulting intervals would contain the true value. We cannot say there is a 95% chance the truth is in this particular interval.
• Bayesian 95% credible interval: There is a 95% probability that the true value is in this interval.

The Bayesian interpretation is what most people think a confidence interval means. It is, frankly, a more intuitive statement.

18.12 Seeing the Two Frameworks Side by Side

Let’s do a final simulation to really see how the two frameworks compare. We are going to “play God” one more time. We will create a population where the true support for the Democrat is 52%. We will then take a sample, and produce both frequentist and Bayesian answers.

set.seed(19104)

#God says: true support is 52%
true.pi <- 0.52
n <- 500

#Take a sample
sample.data <- rbinom(n, 1, true.pi)
p.hat <- mean(sample.data)

cat("Sample proportion:", round(p.hat, 3), "\n")

Sample proportion: 0.518 

cat("Sample size:", n, "\n")

Sample size: 500 

#----- Frequentist -----
se.freq <- sqrt(p.hat * (1 - p.hat) / n)

#95% CI
freq.lower <- p.hat - 1.96 * se.freq
freq.upper <- p.hat + 1.96 * se.freq

#Test H0: pi = 0.5
z.stat <- (p.hat - 0.5) / se.freq
p.val <- 2 * pnorm(-abs(z.stat))

cat("--- Frequentist Results ---\n")

--- Frequentist Results ---

cat("95% CI:", round(freq.lower, 3), "to", round(freq.upper, 3), "\n")

95% CI: 0.474 to 0.562 

cat("P-value (H0: pi = 0.5):", round(p.val, 4), "\n")

P-value (H0: pi = 0.5): 0.4205 

cat("Reject null at alpha = 0.05?", ifelse(p.val < 0.05, "Yes", "No"), "\n")

Reject null at alpha = 0.05? No 

#----- Bayesian -----
#Prior: Beta(50, 50) -- elections are usually close
a.prior <- 50
b.prior <- 50

#Posterior: Beta(a + successes, b + failures)
successes <- sum(sample.data)
failures <- n - successes
a.post <- a.prior + successes
b.post <- b.prior + failures

#Posterior mean
post.mean <- a.post / (a.post + b.post)

#95% credible interval
cred.lower <- qbeta(0.025, a.post, b.post)
cred.upper <- qbeta(0.975, a.post, b.post)

#P(Democrat ahead)
prob.ahead <- 1 - pbeta(0.5, a.post, b.post)

cat("--- Bayesian Results ---\n")

--- Bayesian Results ---

cat("Posterior mean:", round(post.mean, 3), "\n")

Posterior mean: 0.515 

cat("95% Credible Interval:", round(cred.lower, 3), "to", round(cred.upper, 3), "\n")

95% Credible Interval: 0.475 to 0.555 

cat("P(Democrat ahead):", round(prob.ahead, 3), "\n")

P(Democrat ahead): 0.769 

#Visualize
curve(dbeta(x, a.post, b.post), from=0.4, to=0.65,
      xlab="Democratic Vote Share", ylab="Density",
      main="Bayesian Posterior with Frequentist CI",
      lwd=2, col="dodgerblue")
abline(v=true.pi, lty=2, col="firebrick", lwd=2)
abline(v=cred.lower, lty=3, col="dodgerblue")
abline(v=cred.upper, lty=3, col="dodgerblue")
abline(v=freq.lower, lty=3, col="darkorange")
abline(v=freq.upper, lty=3, col="darkorange")
legend("topright",
       c("Posterior", "True value", "Bayesian 95% CI", "Frequentist 95% CI"),
       col=c("dodgerblue","firebrick","dodgerblue","darkorange"),
       lwd=c(2,2,1,1), lty=c(1,2,3,3))

In many cases the frequentist confidence interval and the Bayesian credible interval will be close to one another, especially with moderate to large sample sizes. The real difference is not the numbers – it is what those numbers mean.

18.13 When Does the Prior Really Matter?

I have been saying that the prior becomes irrelevant with enough data. Let’s see where it does matter: small samples. This is where Bayesian methods can really shine.

Imagine you are a campaign manager and you commission a poll of a small state. You can only afford 50 respondents. Twenty-eight (56%) say they will vote for your candidate.

n.small <- 50
y.small <- 28
p.hat.small <- y.small / n.small

#Frequentist CI
se.small <- sqrt(p.hat.small * (1 - p.hat.small) / n.small)
cat("Frequentist 95% CI:",
    round(p.hat.small - 1.96*se.small, 3), "to",
    round(p.hat.small + 1.96*se.small, 3), "\n")

Frequentist 95% CI: 0.422 to 0.698 

That CI runs from about 0.42 to 0.70. Not very useful. It includes values where you are losing and values where you are winning in a landslide.

Now suppose you know, from historical data, that this state is competitive but leans slightly toward the other party. Your prior is centered around 0.48 with moderate certainty:

#Prior: Beta(24, 26) has mean ~0.48
a.prior.small <- 24
b.prior.small <- 26

a.post.small <- a.prior.small + y.small
b.post.small <- b.prior.small + (n.small - y.small)

post.mean.small <- a.post.small / (a.post.small + b.post.small)
cred.lower.small <- qbeta(0.025, a.post.small, b.post.small)
cred.upper.small <- qbeta(0.975, a.post.small, b.post.small)
prob.ahead.small <- 1 - pbeta(0.5, a.post.small, b.post.small)

cat("Bayesian posterior mean:", round(post.mean.small, 3), "\n")

Bayesian posterior mean: 0.52 

cat("Bayesian 95% Credible Interval:",
    round(cred.lower.small, 3), "to", round(cred.upper.small, 3), "\n")

Bayesian 95% Credible Interval: 0.422 to 0.617 

cat("P(Your candidate ahead):", round(prob.ahead.small, 3), "\n")

P(Your candidate ahead): 0.656 

The Bayesian credible interval is narrower because the prior is contributing information. The posterior mean is pulled from 0.56 (the raw data) towards 0.48 (the prior), ending up somewhere in between. This is sometimes called shrinkage – the estimate is “shrunk” towards the prior. Whether this is a good thing depends entirely on whether the prior is reasonable. If the state really does lean towards the other party, incorporating that information produces better estimates. If your prior is wrong, it makes things worse.

This is the core trade-off. Bayesian methods let you incorporate prior knowledge, which can improve estimates when data are sparse. But they require you to have prior knowledge (or at least be willing to specify one), and if that knowledge is wrong, it can hurt.

18.14 The Connection to What You Already Know

Let me tie this all together by connecting Bayesian ideas back to things we covered this semester.

18.14.1 Sampling Distributions and the Posterior

In the frequentist framework, we imagine repeating the study many times and looking at the distribution of estimates across those hypothetical repetitions. That is the sampling distribution. In the Bayesian framework, the posterior distribution plays a similar role – it tells us where we think the true parameter is. But the posterior is a statement about the parameter, not about hypothetical repetitions.

18.14.2 Confidence Intervals and Credible Intervals

A 95% frequentist confidence interval says: across many repetitions, 95% of intervals constructed this way would contain the truth. A 95% Bayesian credible interval says: there is a 95% probability the truth is in this interval. Same idea, different interpretation. And honestly, the Bayesian interpretation is what most people mean when they report a confidence interval.

18.14.3 Hypothesis Testing and Posterior Probabilities

Instead of p-values, Bayesians compute the posterior probability that a hypothesis is true (or that a parameter falls in some range). This is more directly useful for decision-making. “There is a 92% chance our candidate is ahead” is actionable information. “We reject the null at the 0.05 level” is… technically informative but harder to translate into action.

18.14.4 Power and Sample Size

Remember our long discussion about power – the probability of detecting a true effect? Bayesians don’t do power analysis in the same way, but the intuition carries over. With small samples, the posterior is wide (uncertain). With large samples, it is narrow (precise). The posterior width is analogous to what the standard error does in frequentist statistics. The more data you have, the less the prior matters, and the more precisely you can estimate the truth.

18.15 Why I Teach Frequentist Statistics First

At this point you might be asking: if Bayesian statistics gives more intuitive answers, why did we spend an entire semester doing frequentist statistics?

A few reasons:

1. The social sciences overwhelmingly use frequentist methods. If you read a paper in political science, economics, or sociology, you will see p-values and confidence intervals. You need to understand these tools to read the literature.

2. Bayesian computation used to be extremely difficult. It is only in the last 20-30 years that computers have become fast enough to do serious Bayesian analysis routinely. The frequentist framework was developed in part because the math was tractable with pencil and paper.

3. The prior is genuinely controversial. Two researchers with different priors can reach different conclusions from the same data. In some fields, this is seen as unacceptable. In others, it is embraced.

4. With enough data, the two approaches agree. For most of the problems we have tackled this semester – with hundreds or thousands of observations – the frequentist and Bayesian answers are nearly identical. The differences matter most with small samples or strong prior information.

My view is that every well-trained data scientist should understand both frameworks. The frequentist framework gives you a rigorous, prior-free methodology for evaluating evidence. The Bayesian framework gives you a principled way to incorporate prior knowledge and make direct probability statements. They are complementary tools, not competitors.

18.16 A Summary of Key Differences

	Frequentist	Bayesian
Parameters are…	Fixed but unknown	Random variables with distributions
Probability means…	Long-run frequency	Degree of belief
Inference uses…	Only the data	Data + prior beliefs
Key output	P-values, confidence intervals	Posterior distributions, credible intervals
“95% interval” means	95% of such intervals contain truth	95% probability truth is inside
Prior information	Not formally incorporated	Explicitly modeled
With large samples	Standard errors shrink	Posterior concentrates (prior irrelevant)
With small samples	Wide CIs, low power	Prior helps (or hurts)

The fundamental question of this course – how do we learn about populations from samples – does not change. What changes is the machinery we use to express and quantify uncertainty. Frequentist statistics asks how surprising the data are under a hypothesis. Bayesian statistics asks how probable the hypothesis is given the data. Both are valid. Both are useful. And now you know both exist.