11. Bayesian Statistics

Bayesian statistics treats unknown parameters as quantities with probability distributions, and uses Bayes' theorem (Section 3.4) to turn a prior belief and observed data into a posterior distribution that answers questions such as "how likely is it that the new version is better?" directly. This chapter builds the framework (priors, likelihoods, posteriors and the role of conjugacy), shows how to summarize a posterior, make predictions and compare models, and then covers the simulation algorithms (Monte Carlo, importance sampling, MCMC) that make Bayesian analysis practical for realistic models, ending with hierarchical models. These ideas underlie modern A/B testing tools, probabilistic programming languages such as Stan and PyMC, and the regularization methods of Chapter 10.

<a id="111-bayesian-framework"></a>

11.1 Bayesian Framework

This subchapter sets out what makes an analysis Bayesian, how prior and data combine into a posterior, and how to choose a prior that is defensible.

Bayesian vs. Frequentist Inference

Why it matters. Chapters 8 and 9 answered questions such as "what is the conversion rate?" with estimators, confidence intervals and p-values, whose guarantees describe what would happen over many repetitions of the study. Stakeholders usually want something else: "given these data, how probable is it that the rate exceeds 5%?" That is a probability statement about the parameter, and only the Bayesian approach makes it. Knowing which school a method belongs to tells you what its numbers mean, which questions it can answer, and where the two schools disagree in practice.

Two meanings of probability. Section 3.1 introduced two interpretations of probability: long-run relative frequency and degree of belief. Each school of statistics is built on one of them.

• Frequentist. A parameter $\theta$ (a conversion rate, a mean) is a fixed, unknown constant. Only the data are random. Probability statements describe the procedure: "this interval-building method captures $\theta$ in 95% of repeated samples."
• Bayesian. Uncertainty about $\theta$ is described by a probability distribution. Before the data, it is the prior $p(\theta)$. After observing data $y$, Bayes' theorem gives the posterior $p(\theta \mid y)$, and every conclusion is a statement about that distribution: "given the data, $\theta$ lies in this interval with probability 0.95."

Analogy. A frequentist is like a factory inspector who certifies a machine: "this machine produces correct parts 95% of the time", without saying whether a given part is correct. A Bayesian is like a detective who reasons about this case: starting from what was plausible beforehand, the detective weighs the evidence and states how likely each suspect is now.

Aspect	Frequentist	Bayesian
What is random	the data (over hypothetical repetitions)	the parameter, in the sense of uncertainty about it
Information used	the data and the sampling design	the data plus a prior distribution
Typical output	point estimate, standard error, CI, p-value	the posterior distribution and summaries of it
95% interval means	95% of intervals from repeated samples cover $\theta$	$P(a \le \theta \le b \mid y) = 0.95$ for this interval
"Probability $H_0$ is true"	not defined	$P(H_0 \mid y)$, given a prior on the hypotheses
Depends on the stopping rule	yes (p-values, error rates)	no: only the likelihood of the observed data matters
Small samples, many parameters	asymptotic approximations can fail	exact given the model; the prior stabilizes estimates
Main criticism	answers a question that was not asked	results depend on a subjective prior
Computation	usually optimization (maximize a likelihood)	usually integration (often by simulation)

The likelihood principle in action. Two analysts see the same 12 trials: 9 successes and 3 failures. Analyst A planned in advance to run 12 trials (binomial design). Analyst B planned to stop at the 3rd failure (negative binomial design). Both likelihoods are proportional to $\theta^9(1 - \theta)^3$, so by the likelihood principle (Section 8.2) the data carry the same evidence about $\theta$. The two p-values for $H_0: \theta = 0.5$ differ, because each is computed from outcomes that could have occurred under its stopping rule. The Bayesian posterior uses only the observed likelihood, so it is identical for both analysts.

Where the schools agree. With a lot of data and a prior that is not extreme, the posterior is approximately normal, centered at the MLE, with variance equal to the inverse Fisher information (Section 8.2). This result is the Bernstein-von Mises theorem. In that regime Bayesian credible intervals and frequentist confidence intervals nearly coincide numerically, even though their interpretations differ. The schools disagree in practice mainly with small samples, many parameters, sequential designs and hypothesis tests of point nulls (Bayes Factors and Model Comparison, Section 11.2).

import numpy as np
from math import comb
from scipy import stats

# Same data, two stopping rules: 9 successes and 3 failures
k, f = 9, 3
n = k + f

# Frequentist one-sided p-values for H0: p = 0.5 against p > 0.5
p_binom = stats.binom.sf(k - 1, n, 0.5)                    # design A: n = 12 fixed, P(X >= 9)
p_negbin = sum(comb(11, j) for j in range(3)) / 2 ** 11    # design B: stop at 3rd failure, P(N >= 12)
print(f"design A (n fixed at 12):        p-value = {p_binom:.4f}")
print(f"design B (stop at 3rd failure):  p-value = {p_negbin:.4f}")

# Bayesian analysis with a uniform prior: both likelihoods are proportional to p^9 (1-p)^3
post = stats.beta(1 + k, 1 + f)
print(f"posterior in both designs: Beta({1 + k}, {1 + f}), P(p > 0.5 | data) = {post.sf(0.5):.4f}")

# With lots of data, the flat-prior credible interval and the Wald confidence interval agree
for k, n in [(9, 12), (750, 1000)]:
    phat = k / n
    se = np.sqrt(phat * (1 - phat) / n)
    lo, hi = stats.beta(1 + k, 1 + n - k).ppf([0.025, 0.975])
    print(f"k = {k:3d}, n = {n:4d}:  95% Wald CI ({phat - 1.96 * se:.3f}, {phat + 1.96 * se:.3f})   "
          f"95% credible interval ({lo:.3f}, {hi:.3f})")

design A (n fixed at 12):        p-value = 0.0730
design B (stop at 3rd failure):  p-value = 0.0327
posterior in both designs: Beta(10, 4), P(p > 0.5 | data) = 0.9539
k =   9, n =   12:  95% Wald CI (0.505, 0.995)   95% credible interval (0.462, 0.909)
k = 750, n = 1000:  95% Wald CI (0.723, 0.777)   95% credible interval (0.722, 0.776)

The same data give $p = 0.073$ ("not significant at 5%") under one stopping rule and $p = 0.033$ ("significant") under the other, although nothing about the observations differs. The posterior is Beta(10, 4) for both analysts, with $P(\theta > 0.5 \mid \text{data}) = 0.954$. With only 12 trials, the flat-prior credible interval (0.462, 0.909) differs noticeably from the Wald interval (0.505, 0.995), which even runs up against the upper boundary of 1. With 1,000 trials the two intervals agree to the third decimal, as the Bernstein-von Mises theorem predicts.

Key takeaways

• Frequentists treat parameters as fixed and describe the long-run behavior of procedures; Bayesians describe uncertainty about parameters with probability distributions.
• Bayesian inference needs a prior, and in return gives direct probability statements about parameters and hypotheses.
• Bayesian inference obeys the likelihood principle: the stopping rule does not change the posterior, though it changes p-values.
• With large samples and reasonable priors, both approaches give nearly the same numbers (Bernstein-von Mises); they differ most with small samples, many parameters and point-null tests.

🧪 Practice

1. For each statement, say whether it is frequentist, Bayesian, or a misinterpretation: (a) "There is a 95% probability the true mean is in [4.1, 5.3]", computed from a $t$-interval; (b) "If $H_0$ were true, data this extreme would occur 3% of the time"; (c) "Given the data and our prior, the probability that the drug works is 0.88".
2. In the code, change the data to 10 successes and 2 failures. Recompute both p-values and the posterior probability $P(\theta > 0.5 \mid \text{data})$. Does the conclusion of either analyst change at the 5% level?
3. Show that the binomial likelihood $\binom{12}{9}\theta^9(1-\theta)^3$ and the negative binomial likelihood $\binom{11}{2}\theta^9(1-\theta)^3$ give exactly the same posterior for any prior. What happens to the constant factors?
4. Interview: A product manager asks: "Should we switch our A/B testing to Bayesian? I heard it lets us peek at results whenever we want." How do you respond? Hint: the posterior is valid at any time given the prior and model, but if decisions are made by stopping as soon as $P(B > A)$ crosses a threshold, the frequentist error rates of the decision rule still inflate; explain what stays valid and what does not, and mention the role of the prior.

Prior, Likelihood, and Posterior

Why it matters. Section 3.4 applied Bayes' theorem to a few discrete hypotheses (three kinds of coin). In real problems the unknown is a continuous quantity: a resolution rate, a mean, a regression coefficient. The same recipe still works, with densities in place of probabilities: multiply the prior by the likelihood and normalize. Every Bayesian method in this chapter, from conjugate formulas to MCMC, is a way of carrying out this one computation.

Bayes' theorem for parameters. For a parameter $\theta$ and data $y$:

$$\underbrace{p(\theta \mid y)}{\text{posterior}} = \frac{\overbrace{p(y \mid \theta)}^{\text{likelihood}};\overbrace{p(\theta)}^{\text{prior}}}{\underbrace{p(y)}{\text{evidence}}}, \qquad p(y) = \int p(y \mid \theta),p(\theta),d\theta.$$

Piece	What it is	As a function of $\theta$
Prior $p(\theta)$	belief about $\theta$ before seeing $y$	a probability density
Likelihood $p(y \mid \theta)$	probability of the observed $y$ for each value of $\theta$ (Section 8.2)	not a density in $\theta$; need not integrate to 1
Evidence $p(y)$	the probability of the data averaged over the prior (the marginal likelihood)	a constant: does not depend on $\theta$
Posterior $p(\theta \mid y)$	belief about $\theta$ after seeing $y$	a probability density

Because $p(y)$ does not depend on $\theta$, the shape of the posterior is determined by the numerator alone:

$$p(\theta \mid y) ;\propto; p(y \mid \theta),p(\theta).$$

The evidence only rescales the product so that it integrates to 1. This is fortunate, because the integral defining $p(y)$ is usually the hardest part of the computation; most algorithms in Section 11.3 work with the unnormalized posterior. (The evidence does matter for comparing models: Bayes Factors and Model Comparison, Section 11.2.)

The posterior is a compromise. The posterior concentrates where the prior and the likelihood both give reasonable support. A sharp likelihood (much data) dominates a broad prior; a sharp prior (strong previous knowledge) dominates a broad likelihood (little data).

  prior p(theta)             likelihood p(y | theta)        posterior (prior x likelihood, normalized)

        .--.                               .-.                           .-.
      .'    '.                            /   \                         /   \
     /        \                          /     \                      /     \
  __/          \__                    __/       \__                __/       \__
  -----+-------------> theta    ---------+-----------> theta    ----------+--------> theta
      0.70                              0.55                            0.625
   "older teams                  "11 of 20 tickets             between the two, and
    resolve ~70%"                 resolved"                     narrower than either

A worked example on a grid. A new support team resolved 11 of its first 20 tickets on first contact. Older teams resolve about 70%, and that experience is encoded in a Beta(14, 6) prior (mean 0.70, worth about 20 tickets; Section 5.2). The simplest way to compute a posterior is a grid approximation: evaluate prior × likelihood at many values of $\theta$, then divide by the sum (or integral). It works for any prior and likelihood in one or two dimensions.

Posterior mean as a weighted average. For this Beta prior and binomial data the posterior is Beta$(a + k, b + n - k)$ (Conjugate Priors, below), and its mean can be rewritten as

$$E[\theta \mid y] = \frac{a + k}{a + b + n} = (1 - w),\underbrace{\frac{a}{a + b}}{\text{prior mean}} + w,\underbrace{\frac{k}{n}}{\text{MLE}}, \qquad w = \frac{n}{a + b + n}.$$

The weight on the data grows from 0 to 1 as $n$ grows relative to the prior's "sample size" $a + b$.

import numpy as np
from scipy import stats

# First-contact resolution rate of a new support team
a, b = 14, 6                     # prior from older teams: Beta(14, 6), mean 0.70, worth 20 tickets
k, n = 11, 20                    # new team: 11 of 20 tickets resolved on first contact

# Grid approximation: evaluate prior x likelihood on 1,001 points, then normalize
theta = np.linspace(0, 1, 1001)
prior = stats.beta.pdf(theta, a, b)
like = stats.binom.pmf(k, n, theta)                        # likelihood: data fixed, theta varies
unnorm = like * prior
post = unnorm / np.trapezoid(unnorm, theta)                # dividing by p(y) makes it integrate to 1

exact = stats.beta(a + k, b + n - k)                       # conjugate result (next topics)
print(f"grid posterior mean  = {np.trapezoid(theta * post, theta):.4f}")
print(f"exact Beta({a + k}, {b + n - k}) mean = {exact.mean():.4f}")
print(f"max |grid - exact| density difference = {np.max(np.abs(post - exact.pdf(theta))):.1e}")
print(f"evidence p(y) = {np.trapezoid(unnorm, theta):.4f}")

# The posterior mean is a weighted average of the prior mean and the sample proportion
w = n / (n + a + b)
print(f"prior mean {a / (a + b):.3f}, MLE {k / n:.3f}, weight on data {w:.3f} -> "
      f"{(1 - w) * a / (a + b) + w * k / n:.4f}")

# As data accumulate at the same 55% rate, the likelihood takes over
print("\n    n   post. mean   95% credible interval   weight on data")
for n in [0, 20, 100, 1000, 10000]:
    k = round(0.55 * n)
    post_n = stats.beta(a + k, b + n - k)
    lo, hi = post_n.ppf([0.025, 0.975])
    print(f"{n:5d}   {post_n.mean():9.3f}      ({lo:.3f}, {hi:.3f})         {n / (n + a + b):.3f}")

grid posterior mean  = 0.6250
exact Beta(25, 15) mean = 0.6250
max |grid - exact| density difference = 2.0e-14
evidence p(y) = 0.0725
prior mean 0.700, MLE 0.550, weight on data 0.500 -> 0.6250

    n   post. mean   95% credible interval   weight on data
    0       0.700      (0.488, 0.874)         0.000
   20       0.625      (0.472, 0.766)         0.500
  100       0.575      (0.486, 0.662)         0.833
 1000       0.553      (0.522, 0.583)         0.980
10000       0.550      (0.541, 0.560)         0.998

The grid posterior matches the exact Beta(25, 15) to machine precision. Its mean 0.625 lies halfway between the prior mean 0.70 and the observed 0.55, because the prior and the data are each worth 20 tickets. The table shows the likelihood taking over: after 100 tickets the data carry 83% of the weight, and after 10,000 tickets the posterior is concentrated tightly around the observed 55%, whatever the older teams did.

Sequential updating. As in Section 3.4, today's posterior is tomorrow's prior. Updating ticket by ticket, or with all 20 tickets at once, gives the same posterior when observations are conditionally independent given $\theta$. Bayesian analysis is naturally online: there is no separate procedure for "interim" and "final" analyses.

Key takeaways

• The posterior is proportional to likelihood × prior; the evidence $p(y)$ only normalizes it.
• The likelihood is a function of $\theta$ for fixed data and is not a probability density for $\theta$.
• The posterior is a compromise between prior and data; for Beta-binomial, the posterior mean is a weighted average of the prior mean and the MLE with data weight $n/(a + b + n)$.
• Grid approximation computes any one- or two-parameter posterior by brute force; it does not scale to many parameters (Section 11.3).

🧪 Practice

1. With a Beta(2, 2) prior and 7 successes in 10 trials, write the posterior distribution, its mean, and the weight $w$ on the data.
2. Modify the grid code to use a non-conjugate prior: a triangular density peaking at 0.7, $p(\theta) \propto 1 - \lvert\theta - 0.7\rvert / 0.7$ on $[0, 1]$ (set it to 0 where negative). How much does the posterior mean change compared with the Beta(14, 6) prior?
3. Show that updating the Beta(14, 6) prior first with 6 resolved of 10, then with 5 resolved of 10, gives the same posterior as updating with 11 of 20 at once.
4. Interview: Explain to a non-statistician why "the posterior is proportional to the likelihood times the prior" and why we can usually ignore the denominator. Hint: use a picture of two curves multiplied point by point; the denominator only rescales the result to total probability 1.

Choosing Priors

Why it matters. The prior is what makes an analysis Bayesian, and it is the part critics attack first. A good prior encodes real knowledge and stabilizes estimates from small samples. A careless one can dominate the data or put weight on absurd values. Choosing a prior is a modeling decision like choosing a likelihood, and it should be justified, checked and reported in the same way.

Where prior information comes from.

Source	Example	How to encode it
Previous studies or historical data	past A/B tests on similar pages	fit a distribution to past estimates; consider discounting their weight
Physical or logical constraints	a rate is in $[0, 1]$; a variance is positive	choose the support of the prior accordingly
Expert judgment	"defect rates are almost always between 2% and 8%"	elicit quantiles, then fit a distribution
Structure of the problem	many similar groups (stores, schools)	a hierarchical prior learned from the data (Hierarchical Models, Section 11.3)
Regularization goals	effects are usually small	a prior centered at 0 with moderate spread (Noninformative and Weakly Informative Priors, below)

Types of priors.

Type	Idea	Example
Informative	encodes substantial knowledge	Beta(14, 6) from older support teams
Weakly informative	rules out absurd values, otherwise lets the data speak	$N(0, 2.5^2)$ on a standardized logistic coefficient
Noninformative / reference	tries to add as little information as possible	flat prior, Jeffreys prior
Conjugate	chosen so that the posterior has the same family, for convenience	Beta for binomial, Gamma for Poisson
Skeptical	concentrated near "no effect", to demand strong evidence	a narrow prior around a zero treatment effect

Eliciting a prior from an expert. People cannot state $\alpha$ and $\beta$, but they can answer questions about quantities they understand: "what is your best guess?", "what range are you 90% sure of?". Choose a convenient family and find the parameters whose quantiles match the answers. Then show the expert the consequences before using the prior.

Prior predictive checks. A prior on parameters implies a distribution for data: draw $\theta$ from the prior, then draw data given $\theta$. If the simulated data look implausible to someone who knows the domain (negative weights, 90% of users converting), the prior is wrong even if each parameter's prior looked reasonable on paper. This check matters most in models with many parameters, where the joint effect of individually vague priors can be extreme.

Sensitivity analysis. Report how conclusions change under a few reasonable alternative priors. If they hardly change, the data dominate and the prior choice is not critical. If they change materially, say so: the data alone are not enough to settle the question, and the conclusion rests partly on the prior.

import numpy as np
from scipy import stats, optimize

rng = np.random.default_rng(1101)

# Elicitation: an expert is "90% sure the defect rate is between 2% and 8%"
def mismatch(log_ab):
    a, b = np.exp(log_ab)                                   # optimize on the log scale so a, b > 0
    q05, q95 = stats.beta.ppf([0.05, 0.95], a, b)
    return (q05 - 0.02) ** 2 + (q95 - 0.08) ** 2

a, b = np.exp(optimize.minimize(mismatch, x0=np.log([2, 40]), method="Nelder-Mead").x)
print(f"elicited prior: Beta({a:.2f}, {b:.1f}), mean {a / (a + b):.3f}, "
      f"90% interval ({stats.beta.ppf(0.05, a, b):.3f}, {stats.beta.ppf(0.95, a, b):.3f}), worth ~{a + b:.0f} units")

# Prior predictive check: what does this prior say about defects in a batch of 500 units?
p = rng.beta(a, b, 100_000)                                 # 1. draw a rate from the prior
defects = rng.binomial(500, p)                              # 2. draw a batch given that rate
print(f"prior predictive for 500 units: median {np.median(defects):.0f} defects, "
      f"90% range {np.quantile(defects, 0.05):.0f} to {np.quantile(defects, 0.95):.0f}, "
      f"P(0 defects) = {np.mean(defects == 0):.4f}")

# Sensitivity analysis: the same data under three different priors
priors = {"flat Beta(1, 1)": (1, 1), f"expert Beta({a:.1f}, {b:.0f})": (a, b),
          "skeptical Beta(1, 199)": (1, 199)}
for k, n in [(4, 50), (90, 2000)]:
    print(f"\ndata: {k} defects in {n} units (sample rate {k / n:.3f})")
    for name, (a0, b0) in priors.items():
        post = stats.beta(a0 + k, b0 + n - k)
        lo, hi = post.ppf([0.025, 0.975])
        print(f"  {name:22s} posterior mean {post.mean():.4f}   95% CrI ({lo:.4f}, {hi:.4f})")

elicited prior: Beta(5.80, 120.5), mean 0.046, 90% interval (0.020, 0.080), worth ~126 units
prior predictive for 500 units: median 22 defects, 90% range 8 to 42, P(0 defects) = 0.0001

data: 4 defects in 50 units (sample rate 0.080)
  flat Beta(1, 1)        posterior mean 0.0962   95% CrI (0.0326, 0.1888)
  expert Beta(5.8, 121)  posterior mean 0.0556   95% CrI (0.0269, 0.0938)
  skeptical Beta(1, 199) posterior mean 0.0200   95% CrI (0.0066, 0.0406)

data: 90 defects in 2000 units (sample rate 0.045)
  flat Beta(1, 1)        posterior mean 0.0455   95% CrI (0.0368, 0.0550)
  expert Beta(5.8, 121)  posterior mean 0.0451   95% CrI (0.0367, 0.0543)
  skeptical Beta(1, 199) posterior mean 0.0414   95% CrI (0.0334, 0.0501)

The elicited prior Beta(5.8, 120.5) reproduces the expert's 90% range of 2% to 8% exactly and is "worth" about 126 inspected units. Its prior predictive distribution says a batch of 500 should contain about 8 to 42 defects, which the expert can confirm or reject. With 4 defects in 50 units, the three priors give clearly different answers (posterior means 9.6%, 5.6% and 2.0%): 50 units are not enough to overrule an informed prior. With 90 defects in 2,000 units the posteriors nearly agree, and the remaining gap for the skeptical prior signals that it conflicts with the data.

Rules for choosing priors.

• Choose the prior before seeing the data (or from separate data). Tuning a prior until the posterior looks right is the Bayesian version of p-hacking.
• Never give zero prior probability to a possible value (Cromwell's rule, Section 3.4): no amount of data can revive it.
• Think on a scale you understand. A prior that looks vague on the log-odds or log scale can be very informative on the probability or original scale (next topics).
• Report the prior and a sensitivity analysis, so readers with different beliefs can judge the conclusions.

Key takeaways

• Priors come from previous data, constraints, expert judgment, problem structure or regularization goals, and should be justified like any modeling choice.
• Elicit priors through quantities people understand (quantiles, plausible ranges), then fit a distribution.
• Prior predictive checks simulate data from the prior to reveal implausible implications.
• Sensitivity analysis shows how much conclusions depend on the prior; with enough data, reasonable priors converge.

🧪 Practice

1. An expert says a click-through rate is "most likely about 3%, and almost certainly (95%) below 6%". Find a Beta prior that matches both statements approximately (modify the elicitation code).
2. Run a prior predictive check for a regression of adult height (cm) on age with priors $\beta_0 \sim N(0, 100^2)$ and $\beta_1 \sim N(0, 10^2)$ per year, for ages 20 to 70. What fraction of simulated heights are negative or above 250 cm? Propose better priors.
3. Add a fourth prior to the sensitivity analysis: an "optimistic" Beta(1, 999). For which sample size do all four posterior means fall within 0.5 percentage points of each other?
4. Interview: A reviewer says "your Bayesian analysis is subjective because you chose the prior". How do you respond? Hint: likelihoods, test statistics and significance levels are also choices; priors are explicit and checkable; show prior predictive checks and a sensitivity analysis; with enough data reasonable priors agree.

Conjugate Priors

Why it matters. For most combinations of prior and likelihood, the posterior has no closed form, and computing it needs numerical integration or simulation. For a few important models, choosing the prior from the right family makes the posterior land in the same family, with parameters updated by simple arithmetic. These conjugate pairs give exact answers instantly, make the role of the prior transparent (it acts like extra data), and are the building blocks of larger models and of Gibbs sampling (Section 11.3).

Definition. A family of priors is conjugate to a likelihood if, for any data, the posterior belongs to the same family as the prior. Updating then reduces to updating the parameters.

Derivation: Beta-binomial. With $k$ successes in $n$ trials and a Beta$(a, b)$ prior:

$$p(\theta \mid y) \propto \underbrace{\theta^{k}(1 - \theta)^{n - k}}{\text{likelihood}} \cdot \underbrace{\theta^{a - 1}(1 - \theta)^{b - 1}}{\text{prior}} = \theta^{a + k - 1}(1 - \theta)^{b + n - k - 1}.$$

This is the kernel of a Beta$(a + k,; b + n - k)$ density, so that is the posterior; the normalizing constant follows without computing any integral. This is the result quoted in Section 5.2.

Derivation: Gamma-Poisson. With counts $y_1, \dots, y_n \sim$ Poisson$(\lambda)$ and a Gamma$(\alpha, \beta)$ prior (shape $\alpha$, rate $\beta$; Section 5.2):

$$p(\lambda \mid y) \propto \underbrace{\lambda^{\sum y_i} e^{-n\lambda}}{\text{likelihood}} \cdot \underbrace{\lambda^{\alpha - 1} e^{-\beta\lambda}}{\text{prior}} = \lambda^{\alpha + \sum y_i - 1} e^{-(\beta + n)\lambda} ;\Rightarrow; \lambda \mid y \sim \text{Gamma}\Big(\alpha + \sum y_i,; \beta + n\Big).$$

The prior acts like $\beta$ earlier periods of observation containing $\alpha$ events in total.

Normal-normal (known variance). With $y_1, \dots, y_n \sim N(\mu, \sigma^2)$, $\sigma$ known, and prior $\mu \sim N(\mu_0, \tau_0^2)$, the posterior is normal. It is easiest to state with precisions (inverse variances), which add:

$$\frac{1}{\tau_n^2} = \frac{1}{\tau_0^2} + \frac{n}{\sigma^2}, \qquad \mu_n = \tau_n^2\left(\frac{\mu_0}{\tau_0^2} + \frac{n\bar y}{\sigma^2}\right).$$

The posterior mean is a precision-weighted average of the prior mean and the sample mean, and the posterior is always more precise than either source alone.

Common conjugate pairs.

Likelihood	Parameter	Conjugate prior	Posterior	Prior acts like
Binomial / Bernoulli	success probability $\theta$	Beta$(a, b)$	Beta$(a + k, b + n - k)$	$a$ successes and $b$ failures
Poisson	rate $\lambda$	Gamma$(\alpha, \beta)$	Gamma$(\alpha + \sum y_i, \beta + n)$	$\alpha$ events in $\beta$ periods
Exponential	rate $\lambda$	Gamma$(\alpha, \beta)$	Gamma$(\alpha + n, \beta + \sum y_i)$	$\alpha$ lifetimes totaling $\beta$
Normal, $\sigma$ known	mean $\mu$	$N(\mu_0, \tau_0^2)$	$N(\mu_n, \tau_n^2)$ (above)	$\sigma^2/\tau_0^2$ observations at $\mu_0$
Normal, $\mu$ known	variance $\sigma^2$	Inverse-Gamma$(\alpha, \beta)$	Inv-Gamma$(\alpha + n/2, \beta + \frac12\sum(y_i - \mu)^2)$	$2\alpha$ observations
Normal, both unknown	$(\mu, \sigma^2)$	Normal-Inverse-Gamma	Normal-Inverse-Gamma	prior sample size and prior sum of squares
Multinomial	category probabilities $\boldsymbol\pi$	Dirichlet$(\boldsymbol\alpha)$	Dirichlet$(\boldsymbol\alpha + \mathbf{counts})$	$\alpha_j$ observations in category $j$
Linear regression, normal errors	coefficients $\boldsymbol\beta$ ($\sigma$ known)	$N(\mathbf m_0, \mathbf V_0)$	normal	ridge regression when $\mathbf m_0 = 0$ (Section 10.4)

The last row is the Bayesian view of ridge regression mentioned in Section 10.4: with a $N(0, \sigma^2/\lambda)$ prior on each standardized coefficient, the posterior mean (and mode) is the ridge estimate with penalty $\lambda$.

The code applies three conjugate updates and checks the first two against brute-force grid posteriors.

import numpy as np
from scipy import stats

# Gamma-Poisson: calls per hour at a new help desk
a0, b0 = 24, 2                          # prior Gamma(shape 24, rate 2): mean 12 calls/h, "worth" 2 hours
calls = np.array([17, 14, 19, 15, 20])  # 5 observed hours
a1, b1 = a0 + calls.sum(), b0 + len(calls)
post = stats.gamma(a1, scale=1 / b1)    # scipy uses scale = 1 / rate
print(f"Gamma-Poisson: posterior Gamma({a1}, {b1}), mean {post.mean():.2f} calls/h "
      f"(prior mean {a0 / b0:.1f}, sample mean {calls.mean():.1f}), "
      f"95% CrI ({post.ppf(0.025):.2f}, {post.ppf(0.975):.2f})")

# Same posterior by brute force on a grid: prior x likelihood, normalized
lam = np.linspace(0.01, 40, 4000)
unnorm = stats.gamma.pdf(lam, a0, scale=1 / b0) * np.prod(stats.poisson.pmf(calls[:, None], lam), axis=0)
grid = unnorm / np.trapezoid(unnorm, lam)
print(f"  grid check: mean {np.trapezoid(lam * grid, lam):.2f}, "
      f"max density error {np.max(np.abs(grid - post.pdf(lam))):.1e}")

# Normal-Normal (known sigma): caffeine content of a batch, in mg per can
mu0, tau0 = 95.0, 5.0                   # prior from the supplier's specification
sigma = 4.0                             # known measurement error of the lab device
y = np.array([99.1, 103.8, 98.6, 103.3])
prec0, prec_data = 1 / tau0 ** 2, len(y) / sigma ** 2   # precisions add
mu1 = (prec0 * mu0 + prec_data * y.mean()) / (prec0 + prec_data)
sd1 = (prec0 + prec_data) ** -0.5
print(f"Normal-Normal: prior N({mu0}, {tau0}^2), sample mean {y.mean():.2f} (SE {sigma / np.sqrt(len(y)):.1f}) "
      f"-> posterior N({mu1:.2f}, {sd1:.2f}^2), weight on data {prec_data / (prec0 + prec_data):.3f}")

mu = np.linspace(70, 130, 6001)
unnorm = stats.norm.pdf(mu, mu0, tau0) * np.prod(stats.norm.pdf(y[:, None], mu, sigma), axis=0)
grid = unnorm / np.trapezoid(unnorm, mu)
print(f"  grid check: mean {np.trapezoid(mu * grid, mu):.2f}, "
      f"max density error {np.max(np.abs(grid - stats.norm.pdf(mu, mu1, sd1))):.1e}")

# Dirichlet-multinomial: market shares of three plans, uniform Dirichlet(1, 1, 1) prior
counts = np.array([46, 31, 23])
alpha_post = 1 + counts
print(f"Dirichlet-multinomial: posterior Dirichlet{tuple(int(x) for x in alpha_post)}, "
      f"means {np.round(alpha_post / alpha_post.sum(), 3)}")

Gamma-Poisson: posterior Gamma(109, 7), mean 15.57 calls/h (prior mean 12.0, sample mean 17.0), 95% CrI (12.79, 18.63)
  grid check: mean 15.57, max density error 3.8e-14
Normal-Normal: prior N(95.0, 5.0^2), sample mean 101.20 (SE 2.0) -> posterior N(100.34, 1.86^2), weight on data 0.862
  grid check: mean 100.34, max density error 8.3e-16
Dirichlet-multinomial: posterior Dirichlet(47, 32, 24), means [0.456 0.311 0.233]

The help desk prior (12 calls per hour, worth 2 hours of data) and 5 observed hours (17 per hour) combine into a posterior mean of 15.57, which is $(24 + 85)/(2 + 5)$: prior events plus observed events, divided by prior hours plus observed hours. For the caffeine measurements, the data are 6.25 times as precise as the prior ($4/16$ against $1/25$), so they receive 86% of the weight. Both grid checks agree with the formulas to machine precision.

Limits of conjugacy. Conjugate priors exist mainly for exponential-family likelihoods (Section 10.5) with simple structure. Logistic regression, mixture models and most hierarchical models have none. Conjugacy is also a convenience, not a justification: if a Beta prior cannot express what you know (for example, two plausible regions), use a mixture of Betas (still conjugate) or a non-conjugate prior and compute numerically. Many models that are not conjugate overall are conditionally conjugate, piece by piece, which is what Gibbs sampling exploits.

Key takeaways

• A conjugate prior gives a posterior in the same family; updating is arithmetic on the parameters.
• Beta-binomial: add successes to $a$ and failures to $b$. Gamma-Poisson: add events to $\alpha$ and observation periods to $\beta$.
• Normal-normal: precisions add, and the posterior mean is the precision-weighted average of the prior mean and the sample mean.
• Conjugate parameters read as pseudo-data, which makes the strength of a prior easy to judge; ridge regression is the posterior mean under a normal prior.

🧪 Practice

1. Light bulbs have exponential lifetimes with failure rate $\lambda$ per hour. With a Gamma(2, 1000) prior (rate in hours) and five observed lifetimes totaling 4,200 hours, find the posterior of $\lambda$ and its mean. Is the posterior mean lifetime $E[1/\lambda \mid y]$ equal to $1/E[\lambda \mid y]$? (For a Gamma$(\alpha, \beta)$ distribution, $E[1/\lambda] = \beta/(\alpha - 1)$.)
2. Derive the normal-normal posterior by completing the square in $\mu$ in the exponent of prior × likelihood.
3. In the caffeine example, how many measurements would be needed for the data to receive 99% of the weight? What does this say about the "worth" of the prior?
4. Interview: You are building a system that estimates the click-through rate of millions of ads in real time, updating with every impression. Why is a Beta-binomial model attractive, and what is its main weakness? Hint: constant-time updates of two counters, natural shrinkage for new ads; but no sharing of information across similar ads and no drift over time unless you add discounting or a hierarchical structure.

Noninformative and Weakly Informative Priors

Why it matters. Often there is little prior knowledge, or an analysis must look neutral to skeptical readers (regulators, journals). The natural reaction is to use a "flat" or "vague" prior that lets the data speak. This turns out to be subtler than it looks: flatness depends on the scale of the parameter, some flat priors are not proper probability distributions, and very vague priors can produce absurd implications. Weakly informative priors are the modern compromise: they encode only the knowledge that extreme values are implausible.

Flat priors depend on the parameterization. A uniform prior on a success probability $p$ is not uniform on the log-odds $\log\frac{p}{1-p}$: by the change-of-variables formula (Section 4.4), it implies a standard logistic distribution on the log-odds, concentrated near 0. Conversely, a "vague" $N(0, 10^2)$ prior on the log-odds puts most of its mass near $p = 0$ and $p = 1$. There is no single "uninformative" prior: being flat on one scale means being informative on another.

Improper priors. A prior whose integral is infinite, such as $p(\mu) \propto 1$ on the whole real line, is improper. It can still be used if the resulting posterior is proper (integrable). For a normal mean with at least one observation it is, and the posterior gives the familiar $t$-interval (Posterior Summaries, Section 11.2). It is not guaranteed: the Haldane prior Beta(0, 0), $p(\theta) \propto \theta^{-1}(1 - \theta)^{-1}$, gives an improper posterior when no successes (or no failures) are observed. Always check propriety when using an improper prior.

Jeffreys prior. Jeffreys proposed a rule that gives the same answer whatever parameterization you start from:

$$p(\theta) \propto \sqrt{I(\theta)},$$

where $I(\theta)$ is the Fisher information (Section 8.2). The prior is larger where the data are more informative about $\theta$, so that it spreads mass evenly in the sense of "distinguishable" parameter values.

Model	Jeffreys prior	Notes
Binomial $\theta$	Beta$(\frac12, \frac12)$, $\propto \theta^{-1/2}(1-\theta)^{-1/2}$	proper; credible intervals with good frequentist coverage (Credible Intervals, Section 11.2)
Poisson $\lambda$	$\propto \lambda^{-1/2}$	improper; posterior proper after one observation
Normal mean $\mu$	$\propto 1$	flat
Normal SD $\sigma$	$\propto 1/\sigma$	flat on $\log\sigma$; improper

For multi-parameter models the Jeffreys rule can behave badly, and reference priors (Bernardo) refine it. In practice, most analysts now use weakly informative priors instead.

Weakly informative priors. The idea is to use the information you do have about orders of magnitude, without committing to specific values. Typical defaults:

Parameter	Weakly informative prior	Rules out
Logistic regression coefficient, standardized predictor	$N(0, 2.5^2)$ or Cauchy$(0, 2.5)$	odds ratios of $e^{10}$ per SD of the predictor
Linear regression coefficient, standardized $x$ and $y$	$N(0, 1)$	effects much larger than one SD of $y$ per SD of $x$
Standard deviation or group-level scale $\tau$	half-normal or half-Student-$t$ on a sensible scale	absurdly large variability
Correlation matrix	LKJ$(\eta = 2)$	extreme correlations near $\pm 1$

These priors regularize: they behave like ridge (normal prior) or a robust shrinkage (Cauchy, Student-$t$) penalty. A Laplace (double-exponential) prior gives the lasso as its posterior mode (Section 10.4). Their practical payoff is stability in exactly the situations where maximum likelihood fails, such as complete separation in logistic regression (Section 10.5), where the MLE does not exist.

import numpy as np
from scipy import stats, optimize

rng = np.random.default_rng(1105)

# 1. "Flat" depends on the scale: uniform on p is not uniform on the log-odds, and vice versa
p = rng.uniform(0, 1, 1_000_000)
logit = np.log(p / (1 - p))
print(f"p ~ Uniform(0, 1):        P(|log-odds| < 2) = {np.mean(np.abs(logit) < 2):.3f}  (logistic, not flat)")
beta = rng.normal(0, 10, 1_000_000)                        # a "vague" N(0, 10^2) prior on the log-odds
p_implied = 1 / (1 + np.exp(-beta))
print(f"log-odds ~ N(0, 10^2):    P(p < 0.01 or p > 0.99) = {np.mean((p_implied < 0.01) | (p_implied > 0.99)):.3f}")
beta = rng.normal(0, 1.5, 1_000_000)                       # weakly informative alternative
p_implied = 1 / (1 + np.exp(-beta))
print(f"log-odds ~ N(0, 1.5^2):   P(p < 0.01 or p > 0.99) = {np.mean((p_implied < 0.01) | (p_implied > 0.99)):.3f}")

# 2. Zero events in 20 trials under three "noninformative" priors
k, n = 0, 20
for name, (a, b) in {"flat Beta(1, 1)": (1, 1), "Jeffreys Beta(1/2, 1/2)": (0.5, 0.5)}.items():
    post = stats.beta(a + k, b + n - k)
    print(f"{name:24s} posterior mean {post.mean():.4f}, 95% upper bound {post.ppf(0.95):.4f}")
print("Haldane Beta(0, 0)       posterior Beta(0, 20) is improper: the density ~ 1/p cannot be normalized")

# 3. Logistic regression with complete separation: every y = 1 has a larger x than every y = 0
x = np.array([-1.6, -1.1, -0.7, -0.4, -0.1, 0.2, 0.5, 0.9, 1.3, 1.8])
y = np.array([0, 0, 0, 0, 0, 1, 1, 1, 1, 1])

def neg_log_post(b, prior_sd=None):
    eta = b[0] + b[1] * x
    nll = np.sum(np.logaddexp(0, eta) - y * eta)            # Bernoulli negative log-likelihood
    if prior_sd is not None:
        nll += 0.5 * np.sum((b / prior_sd) ** 2)            # normal prior = ridge penalty
    return nll

# No MLE exists: steeper and steeper curves through the gap at x = 0.05 always fit better
for slope in [5, 20, 100, 1000]:
    b = np.array([-0.05 * slope, slope])
    print(f"slope {slope:5d}: log-likelihood = {-neg_log_post(b):.3g}")
fit = optimize.minimize(neg_log_post, [0, 0], args=(2.5,), method="BFGS")
print(f"posterior mode with N(0, 2.5^2) priors: intercept = {fit.x[0]:.2f}, slope = {fit.x[1]:.2f}")

p ~ Uniform(0, 1):        P(|log-odds| < 2) = 0.762  (logistic, not flat)
log-odds ~ N(0, 10^2):    P(p < 0.01 or p > 0.99) = 0.647
log-odds ~ N(0, 1.5^2):   P(p < 0.01 or p > 0.99) = 0.002
flat Beta(1, 1)          posterior mean 0.0455, 95% upper bound 0.1329
Jeffreys Beta(1/2, 1/2)  posterior mean 0.0238, 95% upper bound 0.0905
Haldane Beta(0, 0)       posterior Beta(0, 20) is improper: the density ~ 1/p cannot be normalized
slope     5: log-likelihood = -1.02
slope    20: log-likelihood = -0.0974
slope   100: log-likelihood = -6.12e-07
slope  1000: log-likelihood = -7.18e-66
posterior mode with N(0, 2.5^2) priors: intercept = -0.16, slope = 3.09

A uniform prior on $p$ puts 76% of its mass on log-odds between $-2$ and 2, so it is far from flat on that scale. The "vague" $N(0, 10^2)$ prior on the log-odds says that 65% of the time the probability is below 1% or above 99%, a strong and usually absurd claim. The $N(0, 1.5^2)$ prior is much more neutral on the probability scale. For zero events in 20 trials, the flat and Jeffreys priors give noticeably different answers (95% upper bounds of 13% and 9%), a reminder that "noninformative" priors still matter when data are scarce, and the Haldane prior fails completely. In the separated logistic regression, the log-likelihood keeps approaching 0 as the slope grows, so there is no finite MLE. The weakly informative prior yields a finite, sensible posterior mode (slope 3.09 per SD).

Key takeaways

• A flat prior on one scale is informative on another; no prior is uninformative in every parameterization.
• Improper priors are acceptable only if the posterior is proper; check this, especially with zero counts or many parameters.
• The Jeffreys prior $\propto \sqrt{I(\theta)}$ is invariant to reparameterization; for a binomial proportion it is Beta(1/2, 1/2).
• Weakly informative priors rule out implausible magnitudes, regularize estimates like ridge or lasso penalties, and fix failures of maximum likelihood such as separation.

🧪 Practice

1. Use the change-of-variables formula to derive the density of $\eta = \log\frac{p}{1-p}$ when $p \sim$ Uniform(0, 1), and show it is the standard logistic density $e^{\eta}/(1 + e^{\eta})^2$.
2. Compute the Fisher information for one Poisson observation and verify that the Jeffreys prior is $\propto \lambda^{-1/2}$. What is the posterior after observing counts $y_1, \dots, y_n$?
3. In the separation example, refit the posterior mode with prior SDs of 1, 2.5, 10 and 100. How does the slope estimate depend on the prior SD, and what happens as the SD grows?
4. Interview: A colleague uses $N(0, 1000^2)$ priors on every coefficient of a logistic regression "to be safe and objective". What could go wrong? Hint: the implied prior on probabilities piles up at 0 and 1; separation still produces huge, unstable estimates; and if the model is compared with Bayes factors, very wide priors penalize the larger model heavily (Section 11.2).

<a id="112-bayesian-inference"></a>

11.2 Bayesian Inference

Once the posterior is available, inference means summarizing it, turning it into decisions and predictions, and checking and comparing the models that produced it.

Posterior Summaries

Why it matters. The posterior distribution is the complete answer to a Bayesian analysis, but nobody reads a density in a report. People need a few numbers: a best guess, a measure of uncertainty, and the probability of the events that matter for the decision ("is the mean battery life above the advertised 10.5 hours?"). Choosing and computing these summaries correctly, especially for quantities that are functions of several parameters, is the everyday work of Bayesian inference.

Standard summaries.

Summary	Definition	Use
Posterior mean $E[\theta \mid y]$	center of mass	best estimate under squared error (Bayesian Point Estimation and Loss Functions, below)
Posterior median	50% quantile	best estimate under absolute error; robust for skewed posteriors
Posterior mode (MAP)	highest density point	equals a penalized MLE; available from optimization alone
Posterior SD	spread of the posterior	the Bayesian analogue of a standard error
Credible interval	interval with a stated posterior probability	uncertainty range (Credible Intervals, below)
$P(\theta \in A \mid y)$	posterior probability of a region	direct answer to "how likely is it that...?"

For a symmetric, unimodal posterior, mean, median and mode coincide. For skewed posteriors (variances, rates near 0, ratios) they differ, and you should state which one you report.

Summaries from draws. In practice the posterior is usually represented by a large sample of draws $\theta^{(1)}, \dots, \theta^{(S)}$, from exact sampling or from MCMC (Section 11.3). Every summary is then a simple computation on the draws:

• a mean is the average of the draws; a quantile is a sample quantile;
• a probability $P(\theta > c \mid y)$ is the fraction of draws with $\theta^{(s)} > c$;
• a function of the parameters, such as a ratio, a difference or a percentile $\mu - 1.645\sigma$, is summarized by computing it for each draw and summarizing the results. No delta method (Section 6.5) is needed, and the result is exact up to simulation error.

The last point is one of the biggest practical advantages of the Bayesian approach: uncertainty propagates automatically to any derived quantity.

Marginalization. With several parameters, the summary of one of them integrates over the others: $p(\mu \mid y) = \int p(\mu, \sigma \mid y),d\sigma$ (Section 6.1). With draws this is trivial: just look at the $\mu$ column and ignore the $\sigma$ column. Uncertainty about the nuisance parameter $\sigma$ is then automatically included in the summary of $\mu$.

Monte Carlo error. Summaries computed from $S$ draws carry simulation noise. For independent draws, the Monte Carlo standard error (MCSE) of a posterior mean is (posterior SD)/$\sqrt S$. For MCMC draws, replace $S$ by the effective sample size (Convergence Diagnostics, Section 11.3). Report enough digits to reflect the posterior uncertainty, not the MCSE.

Example. Battery life (hours) of 15 prototype phones, modeled as $N(\mu, \sigma^2)$ with the standard noninformative prior $p(\mu, \sigma^2) \propto 1/\sigma^2$. The joint posterior can be sampled exactly in two steps: $\sigma^2 \mid y$ is a scaled inverse chi-square, and $\mu \mid \sigma^2, y \sim N(\bar y, \sigma^2/n)$. Integrating out $\sigma^2$, the marginal posterior of $\mu$ is a $t_{n-1}$ distribution centered at $\bar y$ with scale $s/\sqrt n$, so the 95% credible interval for $\mu$ equals the classical $t$-interval (Section 8.3).

import numpy as np
from scipy import stats

rng = np.random.default_rng(1106)

# Battery life (hours) of 15 prototype phones; model y ~ N(mu, sigma^2)
y = np.array([11.2, 9.8, 12.5, 10.9, 8.7, 11.8, 10.4, 12.9, 9.5, 11.1, 10.2, 13.4, 9.9, 10.8, 11.6])
n, ybar, s = len(y), y.mean(), y.std(ddof=1)

# Posterior under the standard noninformative prior p(mu, sigma^2) ~ 1/sigma^2:
#   sigma^2 | y ~ (n-1) s^2 / chi^2_{n-1},   mu | sigma^2, y ~ N(ybar, sigma^2 / n)
S = 100_000
sigma2 = (n - 1) * s ** 2 / rng.chisquare(n - 1, S)        # draw sigma^2 first
mu = rng.normal(ybar, np.sqrt(sigma2 / n))                 # then mu given that sigma^2
sigma = np.sqrt(sigma2)

def summarize(name, d):
    lo, hi = np.quantile(d, [0.025, 0.975])
    print(f"{name:26s} mean {d.mean():6.2f}  median {np.median(d):6.2f}  SD {d.std():5.2f}  "
          f"95% CrI ({lo:5.2f}, {hi:5.2f})")

summarize("mu (mean life, h)", mu)
summarize("sigma (SD of life, h)", sigma)
summarize("5th percentile of life, h", mu - 1.645 * sigma)  # any function of the parameters: transform the draws

# Probabilities of statements are just fractions of draws
print(f"P(mu > 10.5 h | data)                 = {np.mean(mu > 10.5):.3f}")
print(f"P(5% of phones die before 9 h | data) = {np.mean(mu - 1.645 * sigma < 9):.3f}")

# Monte Carlo standard error of the posterior mean of mu (draws here are independent)
print(f"MCSE of E[mu | y] = {mu.std() / np.sqrt(S):.4f} h")

mu (mean life, h)          mean  10.98  median  10.98  SD  0.36  95% CrI (10.26, 11.71)
sigma (SD of life, h)      mean   1.38  median   1.34  SD  0.28  95% CrI ( 0.95,  2.06)
5th percentile of life, h  mean   8.71  median   8.79  SD  0.59  95% CrI ( 7.34,  9.63)
P(mu > 10.5 h | data)                 = 0.912
P(5% of phones die before 9 h | data) = 0.662
MCSE of E[mu | y] = 0.0012 h

The posterior mean battery life is 10.98 h with a 95% credible interval of (10.26, 11.71), the same as the $t$-interval (10.26, 11.70) up to simulation error. The posterior of $\sigma$ is right-skewed (mean 1.38 > median 1.34), as is typical for scale parameters. The 5th percentile of battery life, $\mu - 1.645\sigma$, is a function of both parameters; its posterior comes for free from the draws, with a 95% credible interval of (7.34, 9.63) h. The two probability statements answer business questions directly: there is a 91% posterior probability that the mean exceeds 10.5 h, and a 66% probability that at least 5% of phones die before 9 hours.

Key takeaways

• The posterior is the full answer; report a point summary, a measure of spread, credible intervals and the probabilities that matter for the decision.
• Mean, median and mode differ for skewed posteriors; state which one you use.
• With posterior draws, every summary is a sample statistic, and functions of parameters are summarized by transforming each draw.
• Simulation adds Monte Carlo error of about posterior SD / $\sqrt{\text{ESS}}$; use enough draws that it is negligible.

🧪 Practice

1. A posterior for a conversion rate is Beta(30, 270). Compute its mean, median, mode and SD, and $P(\theta > 0.12 \mid y)$, with scipy.stats, and check them with 100,000 draws.
2. In the battery example, compute the posterior distribution of the coefficient of variation $\sigma/\mu$ and give its mean and 95% credible interval.
3. Using the draws, compute the posterior probability that a single new phone lasts more than 12 hours. How does it differ from $P(\mu > 12 \mid y)$? (This anticipates the posterior predictive distribution, below.)
4. Interview: An analyst reports the posterior mean of the ratio of two conversion rates as the ratio of their posterior means. Is that correct? Hint: $E[X/Y] \ne E[X]/E[Y]$ in general (Jensen's inequality, Section 6.5); compute the ratio for each draw and summarize the results.

Credible Intervals

Why it matters. A confidence interval (Section 8.3) does not allow the statement "the parameter lies in this interval with 95% probability", although that is how most people read it. A credible interval does allow it, conditional on the model and the prior. It is also easy to compute for any quantity, including skewed and bounded parameters, where simple confidence intervals can misbehave.

Definition. A $100(1 - \alpha)%$ credible interval (CrI) is any interval $[a, b]$ with

$$P(a \le \theta \le b \mid y) = 1 - \alpha.$$

Many intervals satisfy this; two choices are standard:

	Equal-tailed interval	Highest posterior density (HPD) interval
Definition	$\alpha/2$ posterior probability in each tail	every point inside has higher density than every point outside
Computation	two quantiles	the shortest interval with probability $1 - \alpha$
Width	wider for skewed posteriors	the shortest possible
Invariant to monotone transformations	yes: transform the endpoints	no
Can contain values of low density	yes, in a skewed tail	no
Multimodal posterior	one interval	possibly a union of intervals (the honest summary)

  density                                     equal-tailed: 2.5% cut from each tail,
    |      .-.                                so the low-density right tail is included
    |     /   '.                              and the high-density region near 0 is cut
    |    /      '.
    |   /         '-.                         HPD: a horizontal line cuts the density;
    |--+--------------+--- density level      everything above it (95% of probability)
    |  |/               '-._                  is the interval; the two ends have
    |  |                    '--.__            equal density
    +--+--------------+-------------- theta
       HPD from here  to here

Credible vs confidence intervals.

	Confidence interval	Credible interval
Random quantity	the interval	the parameter (in the sense of uncertainty)
Probability statement	about the procedure over repeated samples	about $\theta$, given these data
Requires	a sampling distribution	a prior
Valid guarantee	coverage over repetitions	coherent given the model and prior

The two often agree numerically, and good credible intervals usually have good frequentist coverage too. The Jeffreys interval for a binomial proportion, the equal-tailed interval from a Beta$(k + \frac12, n - k + \frac12)$ posterior, is a standard recommendation even among frequentists, because the usual Wald interval $\hat p \pm 1.96\sqrt{\hat p(1 - \hat p)/n}$ covers poorly for small $n$ or $p$ near 0 or 1.

Example. A trial observes 2 serious adverse events in 30 patients. With the Jeffreys prior the posterior is Beta(2.5, 28.5), skewed to the right.

import numpy as np
from scipy import stats, optimize

rng = np.random.default_rng(1107)

# 2 serious adverse events in 30 patients, Jeffreys prior Beta(1/2, 1/2)
post = stats.beta(0.5 + 2, 0.5 + 28)

# Equal-tailed interval: 2.5% of posterior probability in each tail
et = post.ppf([0.025, 0.975])

# HPD interval: the shortest interval with 95% probability (slide the lower tail mass from 0 to 5%)
width = lambda lower_tail: post.ppf(lower_tail + 0.95) - post.ppf(lower_tail)
best = optimize.minimize_scalar(width, bounds=(0, 0.05), method="bounded").x
hpd = post.ppf([best, best + 0.95])

# HPD from posterior draws: the shortest window containing 95% of the sorted draws
d = np.sort(post.rvs(200_000, random_state=rng))
m = int(0.95 * len(d))
i = np.argmin(d[m:] - d[:-m])
print(f"equal-tailed 95%: ({et[0]:.4f}, {et[1]:.4f})  width {et[1] - et[0]:.4f}")
print(f"HPD 95% (exact):  ({hpd[0]:.4f}, {hpd[1]:.4f})  width {hpd[1] - hpd[0]:.4f}")
print(f"HPD 95% (draws):  ({d[i]:.4f}, {d[i + m]:.4f})")
print(f"density at the HPD ends: {post.pdf(hpd[0]):.3f} and {post.pdf(hpd[1]):.3f} (equal, as they should be)")

# Frequentist coverage of two 95% intervals for a binomial proportion, n = 40
n, reps = 40, 200_000
print("\ntrue p   Wald CI   Jeffreys CrI")
for p_true in [0.02, 0.05, 0.1, 0.3, 0.5]:
    k = rng.binomial(n, p_true, reps)
    phat = k / n
    half = 1.96 * np.sqrt(phat * (1 - phat) / n)
    wald = np.mean((phat - half <= p_true) & (p_true <= phat + half))
    lo, hi = stats.beta.ppf([[0.025], [0.975]], k + 0.5, n - k + 0.5)
    jeff = np.mean((lo <= p_true) & (p_true <= hi))
    print(f" {p_true:5.2f}    {wald:.3f}      {jeff:.3f}")

equal-tailed 95%: (0.0141, 0.1971)  width 0.1830
HPD 95% (exact):  (0.0057, 0.1748)  width 0.1691
HPD 95% (draws):  (0.0059, 0.1750)
density at the HPD ends: 1.289 and 1.289 (equal, as they should be)

true p   Wald CI   Jeffreys CrI
  0.02    0.553      0.954
  0.05    0.868      0.986
  0.10    0.914      0.970
  0.30    0.929      0.944
  0.50    0.920      0.962

The HPD interval (0.0057, 0.1748) is shorter than the equal-tailed interval (0.0141, 0.1971) and shifted toward 0, where the posterior density is highest. The density is equal at its two ends, as the definition requires, and the draw-based version agrees with the exact one. The coverage table shows that the Jeffreys credible interval keeps close to (or above) 95% across all true rates, while the Wald confidence interval covers only 55% of the time when $p = 0.02$, because with few events it often collapses to a width of 0.

Reporting. State the interval type (equal-tailed or HPD), the probability level, and the prior. A 90% or 89% interval is common in Bayesian practice because tail quantiles of 95% intervals need many draws to be estimated stably.

Key takeaways

• A $1 - \alpha$ credible interval contains $\theta$ with posterior probability $1 - \alpha$, given the model and prior.
• Equal-tailed intervals use quantiles and transform cleanly; HPD intervals are the shortest and contain only the most plausible values.
• For skewed posteriors the two differ; for multimodal posteriors an HPD region can be a union of intervals.
• Well-chosen credible intervals often have good frequentist coverage; the Jeffreys interval beats the Wald interval for proportions.

🧪 Practice

1. A posterior is $N(4.2, 0.8^2)$. Write its 95% equal-tailed and HPD intervals. Why are they the same?
2. For the Beta(2.5, 28.5) posterior, transform the equal-tailed and HPD intervals to the odds scale $\theta/(1 - \theta)$ by transforming their endpoints. Then compute the HPD interval directly on the odds scale from transformed draws. Which interval type is invariant?
3. Extend the coverage simulation to $n = 10$ and $n = 200$. At which sample size does the Wald interval become acceptable for $p = 0.1$?
4. Interview: A dashboard shows a 95% credible interval for a conversion lift of (0.2%, 3.1%) and a colleague says "so it's the same as a confidence interval". When is that roughly true, and when can they differ a lot? Hint: with large samples and weak priors they nearly coincide; with small samples, informative priors, bounded parameters or many comparisons (where hierarchical priors shrink estimates), they can differ substantially.

Bayesian Point Estimation and Loss Functions

Why it matters. Sometimes a single number must be chosen: the failure rate to plan maintenance around, the demand to stock for, the effect size to report. Which summary of the posterior is "best" depends on the consequences of being wrong. Decision theory makes this precise: define a loss function, then choose the estimate that minimizes the expected loss under the posterior. The familiar summaries (mean, median, mode) turn out to be optimal for specific losses, and asymmetric costs call for other quantiles.

The Bayes estimator. A loss function $L(\theta, a)$ measures the cost of reporting $a$ when the truth is $\theta$. The posterior expected loss of $a$ is

$$\rho(a \mid y) = E[L(\theta, a) \mid y] = \int L(\theta, a),p(\theta \mid y),d\theta,$$

and the Bayes estimator is the $a$ that minimizes it. Because the posterior is available, this optimization is done once, for the data at hand.

Loss $L(\theta, a)$	Name	Bayes estimator
$(\theta - a)^2$	squared error	posterior mean
$\lvert\theta - a\rvert$	absolute error	posterior median
$1$ if $\lvert\theta - a\rvert > \delta$, else $0$, with $\delta \to 0$	0-1 loss	posterior mode (MAP)
$c_u(\theta - a)$ if $\theta > a$; $c_o(a - \theta)$ if $\theta \le a$	asymmetric linear (pinball)	posterior quantile at level $c_u/(c_u + c_o)$

Proof for squared error. $E[(\theta - a)^2 \mid y] = \operatorname{Var}(\theta \mid y) + (E[\theta \mid y] - a)^2$, minimized at $a = E[\theta \mid y]$. The asymmetric case gives the same critical-ratio rule as the newsvendor problem of Section 4.2: if under-estimating costs 4 times as much as over-estimating, report the $4/5 = 80$th percentile.

MAP and penalized likelihood. The posterior mode maximizes $\log p(y \mid \theta) + \log p(\theta)$: the log-likelihood plus a penalty. Ridge regression is the MAP estimate under a normal prior, and the lasso under a Laplace prior (Section 10.4). MAP is convenient (only optimization is needed), but it is not invariant to reparameterization, and in high dimensions the mode can be a poor summary because little probability mass lies near it.

Frequentist properties. Bayes estimators are usually biased, because they shrink toward the prior. But they often have lower mean squared error (Section 8.1) than unbiased estimators, especially with small samples, because the reduction in variance outweighs the bias. The price is a higher MSE when the true value is far from where the prior puts its mass.

import numpy as np
from scipy import stats

rng = np.random.default_rng(1108)

# Posterior for a skewed parameter: weekly failure rate of a machine, Gamma(shape 4, rate 2)
post = stats.gamma(4, scale=1 / 2)
theta = post.rvs(400_000, random_state=rng)
cands = np.linspace(0.5, 4, 3501)                           # candidate point estimates

losses = {
    "squared (a - t)^2":        lambda a: ((a - theta) ** 2).mean(),
    "absolute |a - t|":         lambda a: np.abs(a - theta).mean(),
    "0-1 (miss by > 0.05)":     lambda a: 1 - (post.cdf(a + 0.05) - post.cdf(a - 0.05)),  # exact: draws are too noisy
    "under-estimate costs 4x":  lambda a: np.where(theta > a, 4 * (theta - a), a - theta).mean(),
}
theory = {
    "squared (a - t)^2":        ("posterior mean", post.mean()),
    "absolute |a - t|":         ("posterior median", post.median()),
    "0-1 (miss by > 0.05)":     ("posterior mode", (4 - 1) / 2),
    "under-estimate costs 4x":  ("80th percentile", post.ppf(4 / (4 + 1))),
}
for name, loss in losses.items():
    risk = [loss(a) for a in cands[::10]]                   # posterior expected loss of each candidate
    best = cands[::10][np.argmin(risk)]
    label, value = theory[name]
    print(f"{name:25s} best estimate {best:.2f}   {label} = {value:.2f}")

# Frequentist risk (MSE) of the MLE k/n and of the posterior mean under a Beta(2, 2) prior, n = 10
n = 10
k = np.arange(n + 1)
print("\ntrue p   MSE of k/n   MSE of (k+2)/(n+4)")
for p in [0.05, 0.2, 0.5, 0.8]:
    w = stats.binom.pmf(k, n, p)
    mse_mle = np.sum(w * (k / n - p) ** 2)
    mse_bayes = np.sum(w * ((k + 2) / (n + 4) - p) ** 2)
    print(f" {p:4.2f}     {mse_mle:.4f}       {mse_bayes:.4f}")

squared (a - t)^2         best estimate 2.00   posterior mean = 2.00
absolute |a - t|          best estimate 1.84   posterior median = 1.84
0-1 (miss by > 0.05)      best estimate 1.50   posterior mode = 1.50
under-estimate costs 4x   best estimate 2.76   80th percentile = 2.76

true p   MSE of k/n   MSE of (k+2)/(n+4)
 0.05     0.0048       0.0190
 0.20     0.0160       0.0155
 0.50     0.0250       0.0128
 0.80     0.0160       0.0155

For a right-skewed Gamma(4, 2) posterior, the four losses give four different answers: the mean 2.00, the median 1.84, the mode 1.50, and the 80th percentile 2.76 when under-estimating is four times as costly. In each case the brute-force minimizer of the expected loss agrees with the theory. The MSE table shows the trade-off of shrinkage: the posterior mean $(k + 2)/(n + 4)$ beats the MLE $k/n$ for true values of $p$ between about 0.2 and 0.8 (halving the MSE at 0.5), and loses badly at $p = 0.05$, far from the prior's center.

Key takeaways

• A Bayes estimator minimizes the posterior expected loss; the loss function, not habit, should determine which summary to report.
• Squared loss gives the posterior mean, absolute loss the median, 0-1 loss the mode, and asymmetric linear loss a quantile at $c_u/(c_u + c_o)$.
• MAP estimates are penalized maximum likelihood estimates (ridge, lasso), but are not invariant to reparameterization.
• Bayes estimators trade bias for lower variance and often have lower MSE than unbiased estimators, unless the prior is badly placed.

🧪 Practice

1. Show that the posterior median minimizes $E[\lvert\theta - a\rvert \mid y]$ by differentiating with respect to $a$.
2. A Beta(3, 17) posterior describes a defect rate. Compute the Bayes estimates under squared, absolute and 0-1 loss, and under a loss where over-estimating costs twice as much as under-estimating.
3. Extend the MSE table to a grid of $p$ values from 0.01 to 0.99, and find the range of $p$ where the posterior mean beats the MLE for $n = 10$ and for $n = 100$.
4. Interview: A capacity planner asks you for "the" estimate of next month's peak traffic from a Bayesian model. What do you ask before answering? Hint: the relative cost of under-provisioning (outages) and over-provisioning (wasted servers) sets the quantile to report; peak traffic is a prediction, so use the posterior predictive distribution, not the posterior of a parameter.

Posterior Predictive Distribution

Why it matters. Most practical questions are about future observations, not parameters: how many calls will arrive in tomorrow's shift, will the next batch contain a defect, what rent will this apartment fetch? Plugging a point estimate into the model ignores the uncertainty about the parameter and gives prediction intervals that are too narrow. The posterior predictive distribution averages over that uncertainty. It also gives a powerful way to check whether the model fits the data at all.

Definition. For a future observation $\tilde y$ that is conditionally independent of the data given $\theta$:

$$p(\tilde y \mid y) = \int p(\tilde y \mid \theta),p(\theta \mid y),d\theta.$$

It is the model's distribution of new data, averaged over the posterior. By the law of total variance (Section 6.2), its variance has two parts:

$$\operatorname{Var}(\tilde y \mid y) = \underbrace{E[\operatorname{Var}(\tilde y \mid \theta) \mid y]}{\text{noise in the outcome}} + \underbrace{\operatorname{Var}(E[\tilde y \mid \theta] \mid y)}{\text{uncertainty about } \theta}.$$

A plug-in prediction $p(\tilde y \mid \hat\theta)$ keeps only the first term. The second term shrinks as data accumulate, but the first never does.

Simulation. To draw from the posterior predictive distribution, repeat: draw $\theta^{(s)}$ from the posterior, then draw $\tilde y^{(s)}$ from $p(\tilde y \mid \theta^{(s)})$. This works for any model where you can simulate data, and it is how predictions are made from MCMC output.

Closed forms for conjugate models.

Model	Posterior	Posterior predictive for new data
Beta-binomial	Beta$(a', b')$	Beta-binomial$(m, a', b')$ for $m$ new trials
Gamma-Poisson	Gamma$(\alpha', \beta')$	negative binomial with $r = \alpha'$, $p = \beta'/(\beta' + m)$ for the total over $m$ periods
Normal (noninformative prior)	$t$ posterior for $\mu$	$t_{n-1}$ centered at $\bar y$ with scale $s\sqrt{1 + 1/n}$: the classical prediction interval $\bar y \pm t^* s\sqrt{1 + 1/n}$ (regression version in Section 10.1)

Posterior predictive checks. If the model is adequate, data simulated from it should look like the observed data. Choose a test statistic $T$ that captures an aspect the model might get wrong (the variance-to-mean ratio, the number of zeros, the maximum, the autocorrelation), simulate replicated datasets $y^{\text{rep}}$ from the posterior predictive, and compare:

$$p_B = P\big(T(y^{\text{rep}}) \ge T(y) \mid y\big).$$

Values of $p_B$ very close to 0 or 1 mean that the model systematically fails to reproduce that feature. Unlike a significance test, the goal is not to reject the model, but to find how it fails and fix it. Graphical checks, such as overlaying histograms of several replicated datasets on the data, are often more informative than a single $p_B$.

import numpy as np
from scipy import stats

rng = np.random.default_rng(1109)

# Help desk from the conjugate example: posterior for the hourly rate is Gamma(109, rate 7)
a1, b1 = 109, 7
lam_hat = a1 / b1

# Predict the total calls in the next 8-hour shift
plug = stats.poisson(8 * lam_hat)                           # plug-in: pretends lambda is known
pred = stats.nbinom(a1, b1 / (b1 + 8))                      # exact posterior predictive (negative binomial)
lam = rng.gamma(a1, 1 / b1, 200_000)                        # same thing by simulation:
sim = rng.poisson(8 * lam)                                  #   draw lambda, then draw data given lambda
print(f"plug-in Poisson:       mean {plug.mean():.1f}, SD {plug.std():5.2f}, 90% interval {plug.ppf(0.05):.0f} to {plug.ppf(0.95):.0f}")
print(f"posterior predictive:  mean {pred.mean():.1f}, SD {pred.std():5.2f}, 90% interval {pred.ppf(0.05):.0f} to {pred.ppf(0.95):.0f}")
print(f"simulated predictive:  mean {sim.mean():.1f}, SD {sim.std():5.2f}, 90% interval {np.quantile(sim, 0.05):.0f} to {np.quantile(sim, 0.95):.0f}")

# Coverage of 90% intervals for the next shift, after only 3 observed hours (true rate 15, Jeffreys prior)
reps, hits_plug, hits_pred = 20_000, 0, 0
for _ in range(reps):
    obs = rng.poisson(15, 3)
    a, b = 0.5 + obs.sum(), 3
    future = rng.poisson(15 * 8)
    lo, hi = stats.poisson.ppf([0.05, 0.95], 8 * a / b)
    hits_plug += lo <= future <= hi
    lo, hi = stats.nbinom.ppf([0.05, 0.95], a, b / (b + 8))
    hits_pred += lo <= future <= hi
print(f"coverage of nominal 90% intervals: plug-in {hits_plug / reps:.3f}, predictive {hits_pred / reps:.3f}")

# Posterior predictive check: is a Poisson model adequate for daily crash counts?
crashes = rng.negative_binomial(2, 2 / (2 + 4), 60)        # 60 days of overdispersed counts (mean 4)
a, b = 0.5 + crashes.sum(), len(crashes)
T_obs = crashes.var(ddof=1) / crashes.mean()                # test statistic: variance-to-mean ratio
lam = rng.gamma(a, 1 / b, 5000)
reps_data = rng.poisson(lam[:, None], (5000, len(crashes)))  # one replicated dataset per posterior draw
T_rep = reps_data.var(axis=1, ddof=1) / reps_data.mean(axis=1)
print(f"\nvariance/mean: observed {T_obs:.2f}, replicated 95% range "
      f"({np.quantile(T_rep, 0.025):.2f}, {np.quantile(T_rep, 0.975):.2f}), "
      f"posterior predictive p-value {np.mean(T_rep >= T_obs):.4f}")
zeros_rep = (reps_data == 0).sum(axis=1)
print(f"days with zero crashes: observed {np.sum(crashes == 0)}, replicated median {np.median(zeros_rep):.0f}, "
      f"p-value {np.mean(zeros_rep >= np.sum(crashes == 0)):.4f}")

plug-in Poisson:       mean 124.6, SD 11.16, 90% interval 107 to 143
posterior predictive:  mean 124.6, SD 16.34, 90% interval 99 to 152
simulated predictive:  mean 124.6, SD 16.35, 90% interval 99 to 152
coverage of nominal 90% intervals: plug-in 0.625, predictive 0.908

variance/mean: observed 3.77, replicated 95% range (0.68, 1.39), posterior predictive p-value 0.0000
days with zero crashes: observed 8, replicated median 0, p-value 0.0000

For the help desk's next 8-hour shift, the plug-in Poisson prediction has an SD of 11.2 calls, while the posterior predictive SD is 16.3. The extra spread is the uncertainty about $\lambda$: $\sqrt{8E[\lambda] + 64\operatorname{Var}(\lambda)} = \sqrt{124.6 + 142.4}$. The closed-form negative binomial and the two-step simulation agree. The coverage experiment shows why this matters: after only 3 observed hours, nominal 90% plug-in intervals cover the actual shift total just 62.5% of the time, while the posterior predictive intervals cover 90.8%. In the model check, the observed crash counts have a variance 3.8 times their mean and 8 days without crashes. Data replicated from the fitted Poisson model never show either feature ($p_B = 0$ for both), so the Poisson model is inadequate: an overdispersed model such as the negative binomial (Section 10.5) is needed.

Key takeaways

• The posterior predictive distribution $p(\tilde y \mid y) = \int p(\tilde y \mid \theta)p(\theta \mid y),d\theta$ includes both outcome noise and parameter uncertainty.
• Plug-in predictions ignore parameter uncertainty and give intervals that are too narrow, badly so with little data.
• Simulate predictions by drawing a parameter from the posterior, then data given that parameter.
• Posterior predictive checks compare replicated data with the observed data on chosen statistics, revealing how a model fails.

🧪 Practice

1. With a Beta(3, 7) posterior for a success probability, compute the posterior predictive probability that the next trial is a success, and the probability of 0 successes in the next 5 trials. Compare with the plug-in values using the posterior mean.
2. Rerun the coverage experiment with 30 observed hours instead of 3. How close does the plug-in interval's coverage get to 90%? Explain using the variance decomposition.
3. Fit a Gamma-Poisson model separately to weekdays and weekends of simulated data that have a weekly pattern, then design a posterior predictive check that would detect the pattern when a single rate is fitted to all days.
4. Interview: A demand forecasting model produces 95% prediction intervals that contain the actual sales only 70% of the time. List possible causes and how a Bayesian workflow would diagnose each. Hint: plug-in parameter estimates (ignored parameter uncertainty), overdispersion or heavy tails in the outcome, drift over time; use posterior predictive intervals and posterior predictive checks on dispersion and on recent periods.

Bayes Factors and Model Comparison

Why it matters. Hypothesis tests and model selection are where Bayesian and frequentist answers differ most. A p-value cannot say how probable the null hypothesis is, and it can reject a null that the data actually support (Lindley's paradox, below). The Bayes factor measures how much the data shift the odds between two models, and it extends the odds form of Bayes' theorem (Section 3.4) from simple hypotheses to entire models with unknown parameters.

Marginal likelihood. A model $M$ with parameters $\theta$ and prior $p(\theta \mid M)$ predicts data through its marginal likelihood (the evidence of Prior, Likelihood, and Posterior, above):

$$p(y \mid M) = \int p(y \mid \theta, M),p(\theta \mid M),d\theta.$$

It is the average of the likelihood over the prior: how well the model predicted the data before seeing them.

Bayes factor.

$$\text{BF}{10} = \frac{p(y \mid M_1)}{p(y \mid M_0)}, \qquad \underbrace{\frac{P(M_1 \mid y)}{P(M_0 \mid y)}}{\text{posterior odds}} = \text{BF}{10} \times \underbrace{\frac{P(M_1)}{P(M_0)}}{\text{prior odds}}.$$

$\text{BF}{10} > 1$ favors $M_1$, and $\text{BF}{01} = 1/\text{BF}_{10}$. A commonly used scale (Kass and Raftery, 1995):

$\text{BF}_{10}$	Evidence for $M_1$
1 to 3	not worth more than a bare mention
3 to 20	positive
20 to 150	strong
> 150	very strong

Automatic Occam's razor. A flexible model spreads its prior predictions over many possible datasets, so it assigns less probability to any particular one. If a simple model predicted the observed data about as well, it wins, because it "bet" more on what actually happened. Bayes factors therefore penalize complexity without any explicit penalty term.

  p(y | M)   simple model M0: concentrated predictions
     |    ##
     |    ##           complex model M1: spreads probability
     |    ##           over many possible datasets
     |    ##    .......................................
     +----##-----------------------------------------> possible datasets
          ^ observed y here: M0 predicted it better, BF favors M0

Example: is a generator fair? $M_0$: $p = 0.5$ exactly. $M_1$: $p \sim$ Beta(1, 1). For $k$ successes in $n$ trials, the marginal likelihoods are

$$p(y \mid M_0) = \binom{n}{k}0.5^n, \qquad p(y \mid M_1) = \binom{n}{k}\frac{B(1 + k, 1 + n - k)}{B(1, 1)} = \frac{1}{n + 1}.$$

Lindley's paradox. With a large $n$, a result can be "significant" ($p < 0.01$) and yet favor $H_0$ by the Bayes factor. The p-value only asks whether the data are unusual under $H_0$. The Bayes factor asks which model predicted them better. A share of 51.3% in 10,000 trials is unusual for exactly 50%, but under $M_1$, which spread its bets over the whole range from 0 to 1, a share so close to 50% is even less expected.

Sensitivity to the prior. Unlike posteriors, Bayes factors stay sensitive to the prior under $M_1$ however much data there is. A very wide prior spreads its predictions too thin and makes $M_0$ win almost automatically; this is why improper or very vague priors cannot be used for Bayes factors. Always report the prior under each model, and check sensitivity.

BIC approximation. For large $n$, the Laplace approximation of the marginal likelihood gives $\log p(y \mid M) \approx \ell(\hat\theta) - \frac{k}{2}\log n = -\frac12\text{BIC}$, so

$$\log \text{BF}_{10} \approx -\tfrac12(\text{BIC}_1 - \text{BIC}_0).$$

This is the Bayesian justification of BIC mentioned in Section 10.4. The approximation drops terms that do not grow with $n$, including the effect of the prior, so it is rough for small samples.

import numpy as np
from math import comb, log
from scipy import stats
from scipy.special import betaln

# Is a coin (or an A/A test, or a random generator) fair?
#   M0: p = 0.5            M1: p ~ Beta(1, 1), i.e. "any p equally plausible"
def log_bf10(k, n, a=1, b=1):
    log_m1 = log(comb(n, k)) + betaln(a + k, b + n - k) - betaln(a, b)   # marginal likelihood: average over the prior
    log_m0 = log(comb(n, k)) + n * log(0.5)
    return log_m1 - log_m0

print("     k /    n   share   p-value     BF10   P(M1 | y) with 50:50 prior odds")
for k, n in [(7, 10), (60, 100), (62, 100), (5_130, 10_000), (50_500, 100_000)]:
    pval = 2 * stats.binom.sf(k - 1, n, 0.5) if k > n / 2 else 1.0
    bf = np.exp(log_bf10(k, n))
    print(f"{k:6d} / {n:6d}   {k / n:.3f}   {pval:7.4f}   {bf:7.3f}   {bf / (1 + bf):.3f}")

# The Bayes factor depends on the prior width under M1 (the posterior for p barely does)
k, n = 62, 100
print(f"\nk = {k}, n = {n}: BF10 under different priors for p in M1")
for a in [1, 10, 100]:
    label = f"Beta({a}, {a})"
    print(f"  {label:15s} prior SD {np.sqrt(1 / (4 * (2 * a + 1))):.3f}   BF10 = {np.exp(log_bf10(k, n, a, a)):.2f}")

# BIC approximation: log BF10 ~ -(BIC1 - BIC0) / 2
for k, n in [(62, 100), (5_130, 10_000)]:
    ll1 = stats.binom.logpmf(k, n, k / n)                    # M1 maximized at p = k/n, one free parameter
    ll0 = stats.binom.logpmf(k, n, 0.5)                      # M0 has no free parameter
    bic1, bic0 = -2 * ll1 + 1 * log(n), -2 * ll0
    print(f"k = {k}, n = {n}: exact log BF10 = {log_bf10(k, n):+.2f},  BIC approximation = {-(bic1 - bic0) / 2:+.2f}")

     k /    n   share   p-value     BF10   P(M1 | y) with 50:50 prior odds
     7 /     10   0.700    0.3438     0.776   0.437
    60 /    100   0.600    0.0569     0.913   0.477
    62 /    100   0.620    0.0210     2.214   0.689
  5130 /  10000   0.513    0.0096     0.368   0.269
 50500 / 100000   0.505    0.0016     0.588   0.370

k = 62, n = 100: BF10 under different priors for p in M1
  Beta(1, 1)      prior SD 0.289   BF10 = 2.21
  Beta(10, 10)    prior SD 0.109   BF10 = 4.62
  Beta(100, 100)  prior SD 0.035   BF10 = 2.14
k = 62, n = 100: exact log BF10 = +0.79,  BIC approximation = +0.61
k = 5130, n = 10000: exact log BF10 = -1.00,  BIC approximation = -1.22

In the first rows, the Bayes factors stay close to 1 even when the p-value is 0.057 or 0.021: a 62% share in 100 trials raises the odds of $M_1$ only by a factor of 2.2. The last two rows show Lindley's paradox: with $p = 0.0096$ and $p = 0.0016$, the Bayes factor favors the fair model (BF$_{10}$ = 0.37 and 0.59). With 62 of 100, the Bayes factor more than doubles when the prior under $M_1$ narrows from Beta(1, 1) to Beta(10, 10), which concentrates on plausible near-fair values, and falls again for Beta(100, 100), which is almost the same as $M_0$. The BIC approximation is in the right range but off by about 0.2 on the log scale.

Alternatives for model comparison.

Method	Measures	Notes
Bayes factor / posterior model probabilities	which model predicted the data better, a priori	sensitive to priors; answers "which model is true" if one of them is
BIC	approximation to the log marginal likelihood	quick; ignores the prior
Leave-one-out cross-validation (LOO, PSIS-LOO)	expected log predictive density for new data	estimates out-of-sample predictive accuracy; computed from posterior draws
WAIC	approximation to LOO	like AIC (Section 10.4) for Bayesian models
Posterior predictive checks	absolute fit of one model	not a comparison, but often reveals why one model is better
Bayesian model averaging	averages predictions weighted by $P(M \mid y)$	accounts for model uncertainty instead of picking one

As with AIC and BIC (Section 10.4), use Bayes factors when a choice between substantively different models is the goal, and predictive criteria (LOO, WAIC) when the goal is prediction.

Key takeaways

• The marginal likelihood $p(y \mid M)$ averages the likelihood over the prior; the Bayes factor is the ratio of two marginal likelihoods.
• Posterior odds = Bayes factor × prior odds, extending the odds form of Bayes' theorem to models.
• Bayes factors penalize complexity automatically, can disagree with p-values (Lindley's paradox), and remain sensitive to the prior; improper priors cannot be used.
• BIC approximates $-2$ times the log marginal likelihood; LOO and WAIC compare models by predictive accuracy instead.

🧪 Practice

1. With $M_0: p = 0.5$ and $M_1: p \sim$ Beta(1, 1), compute BF$_{10}$ by hand for 9 heads in 10 tosses. With prior odds 1:1, what is $P(M_1 \mid y)$?
2. Use the code to find, for $n = 100$, $1{,}000$ and $10{,}000$, the smallest share $k/n$ for which BF$_{10} > 3$. Compare with the share at which the two-sided p-value crosses 0.05. What happens to the gap as $n$ grows?
3. Show that the log marginal likelihood can be written as a sum of one-step-ahead log predictive probabilities: $\log p(y_1, \dots, y_n) = \sum_i \log p(y_i \mid y_1, \dots, y_{i-1})$. Why does this make the Bayes factor a measure of predictive performance?
4. Interview: A team uses Bayes factors to decide whether a new feature has any effect, with a $N(0, 100^2)$ prior on the effect (in percent) under the alternative. The Bayes factor always favors "no effect". What is going on? Hint: the alternative spreads its predictions over absurdly large effects, so it predicts realistic small effects poorly; use a prior that reflects plausible effect sizes, or report the posterior of the effect instead.

<a id="113-computational-bayesian-methods"></a>

11.3 Computational Bayesian Methods

Real models rarely have conjugate posteriors, so Bayesian computation relies on simulation: drawing samples from the posterior and replacing integrals with averages.

Monte Carlo Integration

Why it matters. Almost every Bayesian quantity is an integral over the posterior: a mean, a probability, a predictive distribution, a marginal likelihood. Outside conjugate models these integrals have no closed form, and grid methods (Section 11.1) become hopeless beyond a few parameters, because the number of grid points grows exponentially with the dimension. Monte Carlo integration replaces the integral with an average over random draws, with an error that shrinks like $1/\sqrt N$ whatever the dimension. It is the foundation for every method in this subchapter.

The idea. By LOTUS (Section 4.3), an expectation is an integral:

$$E[g(\theta)] = \int g(\theta),p(\theta),d\theta ;\approx; \hat I_N = \frac{1}{N}\sum_{i=1}^N g(\theta^{(i)}), \qquad \theta^{(i)} \overset{\text{iid}}{\sim} p.$$

• The law of large numbers (Section 6.5) guarantees $\hat I_N \to E[g(\theta)]$.
• The central limit theorem gives its error: $\hat I_N$ is approximately normal with standard deviation $\sigma_g/\sqrt N$, where $\sigma_g$ is the SD of $g(\theta)$. This Monte Carlo standard error (MCSE) is estimated by the sample SD of the $g(\theta^{(i)})$ divided by $\sqrt N$.
• A probability is the expectation of an indicator: $P(\theta \in A) = E[\mathbf 1{\theta \in A}]$, estimated by the fraction of draws in $A$.

Analogy. To estimate the average depth of a lake, you can survey it on a fine grid, which is accurate but requires $m^2$ measurements for $m$ points per side, or drop a weighted line at $N$ random spots and average. In a lake with 100 dimensions, the grid is impossible; random spots still give an error of order $1/\sqrt N$.

Grid vs Monte Carlo. A grid with $m$ points per axis in $d$ dimensions needs $m^d$ evaluations, and its error (for the midpoint rule) shrinks like $m^{-2} = N^{-2/d}$. Monte Carlo error shrinks like $N^{-1/2}$ in any dimension. Grids win for $d \le 3$ or so; Monte Carlo wins beyond.

Getting the draws: rejection sampling. Monte Carlo needs draws from $p$. The inverse-CDF method (Section 4.4) needs the quantile function. Rejection sampling needs only an unnormalized density $\tilde p(\theta) \propto p(\theta)$, which is exactly what Bayes' theorem provides. Choose a proposal $q$ that is easy to sample and a constant $M$ with $\tilde p(\theta) \le M q(\theta)$ everywhere. Then repeat:

1. Draw $\theta \sim q$ and $u \sim$ Uniform(0, 1).
2. Accept $\theta$ if $u \cdot M q(\theta) < \tilde p(\theta)$; otherwise discard it.

  M q(theta) +-----------------------------------+   envelope (here: a constant times a uniform density)
             |   x  rejected   x        x        |   propose a point uniformly under the envelope:
             |          .-''-.       x           |   keep it if it falls under the target curve,
             |   x    .' o  o '.                 |   discard it otherwise
             |      .' o  o  o  '.     x         |
             |   .-' o o  o  o o  '-.            |   o = accepted: their theta values are exact
             +--'-------------------'------------+       draws from the normalized target
                          target p~(theta)

Accepted draws are exact, independent draws from $p$. The acceptance rate is $\int \tilde p / M$, so a loose envelope wastes most proposals, and in high dimensions it becomes exponentially small. That limitation is what motivates importance sampling and MCMC.

import numpy as np
from scipy import stats

rng = np.random.default_rng(1111)

# 1. Posterior expectation with a known answer: E[p / (1 - p)] under Beta(25, 15) equals 25 / 14
exact = 25 / 14
for N in [100, 10_000, 1_000_000]:
    p = rng.beta(25, 15, N)
    g = p / (1 - p)                                          # the odds, a nonlinear function of p
    print(f"N = {N:>9,}: estimate {g.mean():.4f}, error {g.mean() - exact:+.4f}, MCSE {g.std() / np.sqrt(N):.4f}")

# 2. Same budget (~1 million evaluations), grid vs Monte Carlo: E[cos(X1 + ... + Xd)], X ~ Uniform(0,1)^d
print("\n  d   points/axis   grid error   MC error")
for d, m in [(1, 1_048_576), (2, 1024), (5, 16), (10, 4)]:
    exact = (((np.exp(1j) - 1) / 1j) ** d).real
    mid = (np.arange(m) + 0.5) / m                           # midpoint rule on each axis
    grid = (np.mean(np.exp(1j * mid)) ** d).real             # the integrand factorizes, so the m^d grid sum does too
    x = rng.random((m ** d, d))                              # the same number of random points
    mc = np.cos(x.sum(axis=1)).mean()
    print(f"{d:3d}   {m:11,d}   {abs(grid - exact):10.1e}   {abs(mc - exact):8.1e}")

# 3. Rejection sampling from an unnormalized posterior
#    Conversion rate p with a logit-normal prior (logit p ~ N(-1, 1)) and 12 conversions in 40 visits
def unnorm_post(p):
    prior = stats.norm.pdf(np.log(p / (1 - p)), -1, 1) / (p * (1 - p))   # change of variables to the p scale
    return prior * p ** 12 * (1 - p) ** 28

grid = np.linspace(1e-6, 1 - 1e-6, 100_001)
M = 1.01 * unnorm_post(grid).max()                           # envelope: M x Uniform(0, 1) density lies above the target
prop = rng.random(200_000)                                   # 1. propose from Uniform(0, 1)
keep = rng.random(200_000) * M < unnorm_post(prop)           # 2. accept with probability target / envelope
draws = prop[keep]
dens = unnorm_post(grid)
print(f"\nrejection sampling: acceptance rate {keep.mean():.3f}, {keep.sum():,} draws, "
      f"posterior mean {draws.mean():.4f} (grid: {np.sum(grid * dens) / np.sum(dens):.4f})")

N =       100: estimate 1.7842, error -0.0015, MCSE 0.0612
N =    10,000: estimate 1.7998, error +0.0140, MCSE 0.0064
N = 1,000,000: estimate 1.7866, error +0.0009, MCSE 0.0006

  d   points/axis   grid error   MC error
  1     1,048,576      3.2e-14    3.7e-05
  2         1,024      3.9e-08    4.3e-04
  5            16      5.3e-04    3.8e-04
 10             4      4.9e-03    8.1e-04

rejection sampling: acceptance rate 0.169, 33,796 draws, posterior mean 0.2973 (grid: 0.2971)

The odds $p/(1 - p)$ under a Beta(25, 15) posterior has exact mean $25/14 = 1.7857$. The MCSE falls tenfold for each hundredfold increase in $N$, as $1/\sqrt N$ predicts. The actual errors fluctuate around that scale: at $N = 100$ the error happened to be small, and at $N = 10{,}000$ it is 2.2 MCSEs, which happens in about 3% of runs. The MCSE is a standard deviation, not a bound. In the second table, the grid is unbeatable in one or two dimensions, but with the same budget of about a million evaluations it has only 4 points per axis in 10 dimensions, and its error is 6 times that of Monte Carlo. Rejection sampling accepts 17% of uniform proposals and reproduces the grid posterior mean of a non-conjugate logit-normal model.

Key takeaways

• Monte Carlo integration estimates $E[g(\theta)]$ by the average of $g$ over draws; probabilities are averages of indicators.
• The error is about $\sigma_g/\sqrt N$ (the MCSE), independent of the dimension: 100 times more draws for one more correct digit.
• Grids beat Monte Carlo in one to three dimensions and lose badly beyond.
• Rejection sampling turns draws from an envelope $Mq$ into exact draws from an unnormalized target, but its acceptance rate collapses in high dimensions.

🧪 Practice

1. Estimate $\pi$ by drawing points uniformly in the unit square and counting the fraction inside the quarter circle. How many draws do you need for an MCSE of 0.001?
2. Estimate $P(\theta_A > \theta_B)$ where $\theta_A \sim$ Beta(41, 961) and $\theta_B \sim$ Beta(52, 950), and give its MCSE. How many draws would make the MCSE below 0.001?
3. In the rejection sampler, replace the uniform proposal with a Beta(4, 9) proposal and compute a valid $M$ on a grid. How does the acceptance rate change? Why must $M q$ lie above the target everywhere, including the tails?
4. Interview: Why does Monte Carlo error not depend on the dimension, while the cost of a grid explodes with dimension? Is Monte Carlo then always "dimension free" in practice? Hint: the CLT error depends only on $\operatorname{Var}(g)$ and $N$; but obtaining independent draws from a high-dimensional posterior is itself the hard part, which is why MCMC exists, and $\operatorname{Var}(g)$ can grow with the dimension.

Importance Sampling

Why it matters. Plain Monte Carlo needs draws from the target distribution itself, and it is hopeless for rare events: estimating a probability of $10^{-5}$ needs millions of draws to see a handful of hits. Importance sampling draws from a different, convenient distribution and corrects with weights. It turns rare-event estimation from impossible to routine, lets you reuse draws from one posterior to approximate another (for example after changing the prior), and is the basis of modern tools such as PSIS-LOO cross-validation (Section 11.2) and particle filters.

The identity. For any proposal density $q$ that is positive wherever $g(\theta)p(\theta) \ne 0$:

$$E_p[g(\theta)] = \int g(\theta)\frac{p(\theta)}{q(\theta)},q(\theta),d\theta = E_q\big[g(\theta),w(\theta)\big], \qquad w(\theta) = \frac{p(\theta)}{q(\theta)}.$$

Draw $\theta^{(i)} \sim q$ and average $g(\theta^{(i)})w(\theta^{(i)})$. The weights correct for drawing too often where $q > p$ and too rarely where $q < p$.

Self-normalized importance sampling. For a posterior known only up to a constant, $\tilde p = p(y \mid \theta)p(\theta)$, use unnormalized weights $\tilde w_i = \tilde p(\theta^{(i)})/q(\theta^{(i)})$ and

$$\hat E[g] = \frac{\sum_i \tilde w_i,g(\theta^{(i)})}{\sum_i \tilde w_i}.$$

The unknown constant cancels. The estimator has a small bias that vanishes as $N$ grows. As a by-product, $\frac1N\sum_i \tilde w_i$ estimates the marginal likelihood $p(y)$ (Section 11.2).

Choosing the proposal. The variance of the estimate depends on how variable the weights are.

Proposal	What happens
close to the target (or to $\lvert g \rvert p$)	weights nearly equal; very efficient
wider, with heavier tails than the target	some wasted draws, but stable
narrower or lighter-tailed than the target	rare draws in the target's tails get huge weights; the estimate is dominated by a few draws, and its variance can be infinite
for a rare event ${\theta \in A}$	center the proposal on $A$ (the optimal $q$ is proportional to $\lvert g \rvert p$)

Effective sample size. A diagnostic of weight degeneracy:

$$\text{ESS} = \frac{\left(\sum_i w_i\right)^2}{\sum_i w_i^2},$$

which equals $N$ for equal weights and 1 when a single draw carries all the weight. A small ESS, or a single draw with a large share of the total weight, means the estimate is unreliable. In high dimensions, almost any fixed proposal gives an ESS near 1, because the ratio $p/q$ is a product over many coordinates: importance sampling alone does not scale to large models.

import numpy as np
from scipy import stats

rng = np.random.default_rng(1112)

def ess(w):
    return w.sum() ** 2 / np.sum(w ** 2)                     # effective sample size of a weighted sample

# 1. A rare event: P(Z > 4) for Z ~ N(0, 1)
exact = stats.norm.sf(4)
N = 100_000
z = rng.normal(0, 1, N)
print(f"exact P(Z > 4)        = {exact:.3e}")
print(f"plain Monte Carlo     = {np.mean(z > 4):.3e}   ({np.sum(z > 4)} hits in {N:,} draws)")
x = rng.normal(4, 1, N)                                      # proposal q = N(4, 1), centered on the region that matters
w = stats.norm.pdf(x) / stats.norm.pdf(x, 4, 1)              # weights p(x) / q(x)
est = np.mean((x > 4) * w)
print(f"importance sampling   = {est:.3e}   (relative SE {np.std((x > 4) * w) / np.sqrt(N) / est:.2%})")

# 2. Self-normalized importance sampling from an unnormalized posterior
#    Normal data with a Cauchy(0, 1) prior on the mean (non-conjugate), sigma = 1 known
y = np.array([2.1, 1.4, 2.9, 1.8, 2.5])
log_post = lambda mu: stats.cauchy.logpdf(mu) + np.sum(stats.norm.logpdf(y[:, None], mu, 1), axis=0)

grid = np.linspace(-5, 8, 20_001)
g = np.exp(log_post(grid) - log_post(grid).max())
print(f"\ngrid posterior mean   = {np.sum(grid * g) / g.sum():.4f}")

for name, q in [("prior Cauchy(0, 1)", stats.cauchy(0, 1)),
                ("t_4 at MLE, scale 0.5", stats.t(4, y.mean(), 0.5)),
                ("normal at MLE, scale 0.15", stats.norm(y.mean(), 0.15))]:
    mu = q.rvs(10_000, random_state=rng)
    logw = log_post(mu) - q.logpdf(mu)
    w = np.exp(logw - logw.max())                            # subtract the max before exponentiating, for stability
    print(f"proposal {name:26s} mean {np.sum(w * mu) / w.sum():.4f}   ESS {ess(w):7.0f} of 10,000"
          f"   largest weight share {w.max() / w.sum():.4f}")

exact P(Z > 4)        = 3.167e-05
plain Monte Carlo     = 3.000e-05   (3 hits in 100,000 draws)
importance sampling   = 3.129e-05   (relative SE 0.68%)

grid posterior mean   = 1.9786
proposal prior Cauchy(0, 1)         mean 1.9837   ESS    1009 of 10,000   largest weight share 0.0014
proposal t_4 at MLE, scale 0.5      mean 1.9820   ESS    8525 of 10,000   largest weight share 0.0002
proposal normal at MLE, scale 0.15  mean 2.0504   ESS     609 of 10,000   largest weight share 0.0207

Plain Monte Carlo saw only 3 draws above 4 in 100,000, so its estimate is little more than a guess. The proposal $N(4, 1)$ puts half of its draws in the rare region, and the weighted estimate has a relative standard error under 1% with the same number of draws. For the Cauchy-prior posterior, all three proposals give self-normalized estimates, but their quality differs. The prior proposal is valid but wasteful (ESS about 1,000 of 10,000, because most prior draws fall where the likelihood is negligible). A Student-$t$ proposal centered at the MLE, slightly wider than the posterior, gives an ESS of 8,525. A normal proposal narrower than the posterior (SD 0.15 against 0.46) misses its tails: a single draw carries 2% of the total weight, the ESS drops to 609, and the estimate 2.05 is clearly off the true 1.98.

Key takeaways

• Importance sampling estimates $E_p[g]$ by drawing from $q$ and weighting by $p/q$; self-normalized weights work with unnormalized posteriors.
• The proposal must cover the target, with tails at least as heavy; a too-narrow proposal gives unstable, possibly infinite-variance estimates.
• The effective sample size $(\sum w)^2/\sum w^2$ and the largest weight share diagnose weight degeneracy.
• Importance sampling excels at rare events and at reweighting existing draws, but degenerates in high dimensions.

🧪 Practice

1. Use importance sampling with an exponential proposal shifted to start at 5 to estimate $P(Z > 5)$ for $Z \sim N(0, 1)$. Compare the relative standard error with that of the $N(5, 1)$ proposal.
2. You have 10,000 posterior draws of a conversion rate under a Beta(1, 1) prior and the data 30 successes in 200 trials. Without resampling, reweight them to obtain the posterior mean under a Beta(5, 45) prior. What are the weights? Check against the exact conjugate answer.
3. Estimate the marginal likelihood $p(y)$ of the Cauchy-prior model by averaging the unnormalized weights with the $t_4$ proposal, and compare with the grid value.
4. Interview: Your importance-sampling estimate of a risk metric changes a lot each time you rerun it with a new seed, although you use a million draws. What do you check and how do you fix it? Hint: look at the ESS and the largest weights; the proposal probably has lighter tails than the target or is centered in the wrong place; widen it, use a heavier-tailed (Student-$t$) or mixture proposal, or switch to an adaptive method.

Metropolis-Hastings Algorithm

Why it matters. Rejection and importance sampling need a proposal that resembles the whole posterior, which is impossible to find for models with dozens or thousands of parameters. Markov chain Monte Carlo (MCMC) avoids this: it builds a random walk that explores the posterior locally, one step at a time, needs only the unnormalized posterior $p(y \mid \theta)p(\theta)$, and, in the long run, visits each region in proportion to its posterior probability. The Metropolis-Hastings (MH) algorithm, from 1953 and 1970, is the general recipe behind MCMC and was named one of the top 10 algorithms of the 20th century.

From Markov chains to sampling. Section 5.5 showed that an irreducible, aperiodic Markov chain converges to its stationary distribution $\pi$ from any start, and that the long-run fraction of time it spends in each state equals $\pi$. MCMC reverses the question: instead of finding $\pi$ for a given chain, design a chain whose stationary distribution is the posterior. Then simulate it and treat its states as (dependent) draws.

The algorithm. Given the current state $\theta$ and a proposal distribution $q(\theta' \mid \theta)$:

1. Propose $\theta' \sim q(\theta' \mid \theta)$.
2. Compute the acceptance probability $$\alpha(\theta, \theta') = \min\left(1, ; \frac{\tilde p(\theta'),q(\theta \mid \theta')}{\tilde p(\theta),q(\theta' \mid \theta)}\right).$$
3. With probability $\alpha$ move to $\theta'$; otherwise stay at $\theta$ (and record $\theta$ again).

The unknown normalizing constant cancels in the ratio, so only $\tilde p = p(y \mid \theta)p(\theta)$ is needed. For a symmetric proposal, $q(\theta' \mid \theta) = q(\theta \mid \theta')$, such as a normal step around the current point, the $q$ terms cancel too. This is the original random-walk Metropolis algorithm: always accept uphill moves, and accept downhill moves with probability equal to the density ratio. The factor $q(\theta \mid \theta')/q(\theta' \mid \theta)$ is the Hastings correction for asymmetric proposals.

Why it works: detailed balance. For $\theta \ne \theta'$, the probability flow from $\theta$ to $\theta'$ in equilibrium is

$$p(\theta),q(\theta' \mid \theta),\alpha(\theta, \theta') = \min\big(p(\theta),q(\theta' \mid \theta),; p(\theta'),q(\theta \mid \theta')\big),$$

which is symmetric in $\theta$ and $\theta'$. So the chain satisfies detailed balance with respect to $p$, and $p$ is stationary (Section 5.5). If $q$ can reach every region with positive density, the chain is irreducible, and the ergodic theorem guarantees that averages over the chain converge to posterior expectations.

Tuning the step size. For a random-walk proposal $\theta' = \theta + \varepsilon z$ with $z \sim N(0, I)$:

  step too small: accepts ~95%,      step about right: accepts      step too large: accepts ~5%,
  but crawls (random walk)           ~25-50%, explores quickly       gets stuck for long stretches
  theta                              theta                           theta
    |        ..''''..                  | .  ' .  '. ' .  .'  .       |    ____      ______
    |   ..'''        ''..              |' '. ' '. '  ' ' '. ' '      |___|    |____|      |___
    |..'                 '''           |                             |
    +-------------------- iteration    +-------------------- iter    +-------------------- iter

Theory for well-behaved targets gives optimal acceptance rates of about 0.44 in one dimension and 0.234 in high dimensions (Roberts, Gelman and Gilks, 1997), reached with a step size of about $2.38/\sqrt d$ times the posterior SD. Practical samplers adapt the step size during a warmup (burn-in) phase, which is then discarded.

The price: autocorrelation. Consecutive MCMC draws are correlated, so $S$ draws carry less information than $S$ independent ones. The effective sample size

$$\text{ESS} = \frac{S}{1 + 2\sum_{k=1}^{\infty}\rho_k},$$

where $\rho_k$ is the autocorrelation at lag $k$, is the number of independent draws with the same precision for estimating a mean. The MCSE of a posterior mean is (posterior SD)$/\sqrt{\text{ESS}}$.

Example. A normal mean with known $\sigma = 1$, five observations, and a Cauchy(0, 1) prior. A heavy-tailed prior like this is a common robust choice, but it is not conjugate, so the posterior has no closed form. The chains start far away, at $\mu = -5$.

import numpy as np

rng = np.random.default_rng(1113)

def ess(x):
    """Effective sample size: n / (1 + 2 * sum of autocorrelations), summed until they first turn negative."""
    x = x - x.mean()
    n = len(x)
    f = np.fft.rfft(x, 2 * n)
    acf = np.fft.irfft(f * np.conj(f))[:n] / (x @ x)        # autocorrelation at every lag, via FFT
    neg = np.argmax(acf[1:] < 0) + 1                         # first lag with a negative autocorrelation
    return n / (1 + 2 * acf[1:neg].sum())

# Target: posterior of a normal mean with a Cauchy(0, 1) prior (sigma = 1 known), known only up to a constant
y = np.array([2.1, 1.4, 2.9, 1.8, 2.5])
def log_target(mu):
    return -np.log1p(mu ** 2) - 0.5 * np.sum((y - mu) ** 2)   # log prior + log likelihood, constants dropped

def metropolis(step, n_iter=20_000, start=-5.0):
    chain = np.empty(n_iter)
    mu, lp = start, log_target(start)
    accepted = 0
    for t in range(n_iter):
        prop = mu + step * rng.normal()                      # symmetric random-walk proposal
        lp_prop = log_target(prop)
        if np.log(rng.random()) < lp_prop - lp:              # accept with probability min(1, ratio), on the log scale
            mu, lp = prop, lp_prop
            accepted += 1
        chain[t] = mu                                        # a rejection repeats the current value
    return chain, accepted / n_iter

grid = np.linspace(-5, 8, 20_001)
dens = np.exp([log_target(m) for m in grid])
dens /= dens.sum()
m = np.sum(grid * dens)
print(f"grid answer: posterior mean {m:.3f}, SD {np.sqrt(np.sum((grid - m) ** 2 * dens)):.3f}\n")

print("step   acceptance   mean    SD      ESS (of 19,000 kept)")
for step in [0.05, 1.0, 10.0]:
    chain, acc = metropolis(step)
    kept = chain[1000:]                                      # discard the first 1,000 as burn-in
    print(f"{step:5.2f}     {acc:.3f}     {kept.mean():.3f}   {kept.std():.3f}   {ess(kept):8.0f}")

# Asymmetric proposals need the Hastings correction: sample Gamma(3, 1) with multiplicative steps
def mh_positive(correct, n_iter=200_000):
    th, out = 1.0, np.empty(n_iter)
    for t in range(n_iter):
        prop = th * np.exp(0.8 * rng.normal())               # log-normal step: always positive, not symmetric
        log_r = (2 * np.log(prop) - prop) - (2 * np.log(th) - th)   # Gamma(3, 1) log density, up to a constant
        if correct:
            log_r += np.log(prop / th)                       # q(th | prop) / q(prop | th) = prop / th
        if np.log(rng.random()) < log_r:
            th = prop
        out[t] = th
    return out[5000:].mean()

print(f"\nGamma(3, 1) target, true mean 3: with Hastings correction {mh_positive(True):.3f}, "
      f"without {mh_positive(False):.3f}")

grid answer: posterior mean 1.979, SD 0.457

step   acceptance   mean    SD      ESS (of 19,000 kept)
 0.05     0.953     1.911   0.493         39
 1.00     0.472     1.982   0.457       4416
10.00     0.058     1.957   0.459        807

Gamma(3, 1) target, true mean 3: with Hastings correction 2.992, without 2.003

All three chains target the same posterior (mean 1.979, SD 0.457), but their efficiency differs by a factor of 100. With steps of 0.05, 95% of proposals are accepted, but each moves the chain very little. The 19,000 kept draws are worth only 39 independent ones, and the estimated mean (1.911) is still off. With steps of 10, almost every proposal lands in a region of negligible density and is rejected: 94% of the time the chain stands still. A step of 1.0 (about twice the posterior SD) accepts 47% and gives an ESS of 4,416. The last line shows the role of the Hastings correction. Log-normal multiplicative proposals keep a positive parameter positive, but they are not symmetric. With the correction the chain samples Gamma(3, 1) correctly (mean 2.99); without it, it silently samples the wrong distribution, Gamma(2, 1), with mean 2.

Key takeaways

• MCMC builds a Markov chain whose stationary distribution is the posterior; it needs only the unnormalized posterior.
• Metropolis-Hastings proposes a move and accepts it with probability $\min(1, \frac{\tilde p(\theta')q(\theta \mid \theta')}{\tilde p(\theta)q(\theta' \mid \theta)})$; detailed balance makes the posterior stationary.
• Asymmetric proposals need the Hastings correction; omitting it samples the wrong distribution without any visible error.
• Step size controls efficiency: aim for acceptance rates of roughly 0.2 to 0.5, and judge chains by effective sample size, not raw length.

🧪 Practice

1. Starting from $\theta = 0$ with target $\tilde p(\theta) = e^{-\theta^2/2}$, the proposal is $\theta' = 1.5$. What is the acceptance probability? If accepted, what is the acceptance probability of proposing a return to 0?
2. Run the three samplers in the code for 2,000 iterations each and plot the trace (iteration vs $\mu$) of each one. Which chain still shows the effect of the starting value $\mu = -5$, and for how long?
3. Write a random-walk Metropolis sampler for a Bayesian logistic regression with one predictor (intercept and slope, $N(0, 2.5^2)$ priors), using the separated data from Noninformative and Weakly Informative Priors (Section 11.1). Tune a separate step size for each coefficient and report acceptance rates and ESS.
4. Interview: Explain why MCMC can sample from a posterior whose normalizing constant is unknown, and why its draws are not independent. Hint: only ratios of the target enter the acceptance probability; each draw is generated from the previous one, and rejections repeat values, so the effective sample size is smaller than the number of iterations.

Gibbs Sampling

Why it matters. Tuning proposal distributions is a chore, and random-walk steps are inefficient when there are many parameters. In many models (normal models, mixture models, hierarchical models, models with missing data) the distribution of each parameter given all the others is a standard distribution, even though the joint posterior is not. Gibbs sampling exploits this: it updates one parameter (or block) at a time by drawing directly from its full conditional distribution. Every draw is accepted, and there is nothing to tune. It powered the first wave of applied Bayesian statistics (BUGS, JAGS).

The algorithm. For parameters $\theta = (\theta_1, \dots, \theta_p)$, repeat for $t = 1, 2, \dots$:

$$\begin{aligned} \theta_1^{(t)} &\sim p(\theta_1 \mid \theta_2^{(t-1)}, \theta_3^{(t-1)}, \dots, \theta_p^{(t-1)}, y) \ \theta_2^{(t)} &\sim p(\theta_2 \mid \theta_1^{(t)}, \theta_3^{(t-1)}, \dots, \theta_p^{(t-1)}, y) \ &;;\vdots \ \theta_p^{(t)} &\sim p(\theta_p \mid \theta_1^{(t)}, \dots, \theta_{p-1}^{(t)}, y) \end{aligned}$$

Each full conditional $p(\theta_j \mid \theta_{-j}, y)$ is proportional to the joint posterior viewed as a function of $\theta_j$ alone, so it can be found by dropping every factor that does not involve $\theta_j$. Each Gibbs step is a Metropolis-Hastings step whose proposal is the full conditional, and its acceptance probability is exactly 1. The joint posterior is therefore stationary.

Conditional conjugacy. A model can be non-conjugate as a whole and conjugate piece by piece. For normal data with a normal prior on $\mu$ and an inverse-gamma prior on $\sigma^2$ (the semi-conjugate prior), the joint posterior is not a standard family, but:

Full conditional	Distribution
$\mu \mid \sigma^2, y$	normal: the normal-normal update of Conjugate Priors (Section 11.1) with $\sigma^2$ fixed
$\sigma^2 \mid \mu, y$	Inverse-Gamma$\left(a_0 + \frac n2,; b_0 + \frac12\sum(y_i - \mu)^2\right)$

Data augmentation. Gibbs sampling is especially natural when some quantities are unobserved: missing values, the cluster label of each point in a mixture model, or the latent utility in probit regression. Treat them as extra parameters and alternate: draw the missing pieces given the parameters, then the parameters given the completed data. It is a stochastic relative of the EM algorithm.

Weakness: correlated parameters. Gibbs moves parallel to the coordinate axes. When two parameters are strongly correlated in the posterior, the full conditionals are narrow, and the chain zigzags slowly along the diagonal ridge. For a bivariate normal with correlation $\rho$, the lag-1 autocorrelation of each coordinate is $\rho^2$.

  theta2                                    rho = 0.99: each conditional draw can only move
    |                         _/            a short distance along one axis, so the chain
    |                      _|/              needs many steps to travel along the ridge
    |                    _|/
    |                  _|/                  fixes: reparameterize (center predictors),
    |                _|/                    update correlated parameters together (blocking),
    |              _|/                      or use gradient-based samplers (HMC)
    +---------------------------- theta1

import numpy as np

rng = np.random.default_rng(1114)

def ess(x):
    """Effective sample size (same helper as in the Metropolis-Hastings example)."""
    x = x - x.mean()
    n = len(x)
    f = np.fft.rfft(x, 2 * n)
    acf = np.fft.irfft(f * np.conj(f))[:n] / (x @ x)
    neg = np.argmax(acf[1:] < 0) + 1
    return n / (1 + 2 * acf[1:neg].sum())

# Battery-life data with semi-conjugate priors: mu ~ N(10, 2^2), sigma^2 ~ Inverse-Gamma(2, 2)
y = np.array([11.2, 9.8, 12.5, 10.9, 8.7, 11.8, 10.4, 12.9, 9.5, 11.1, 10.2, 13.4, 9.9, 10.8, 11.6])
n, ybar = len(y), y.mean()
mu0, tau0_2, a0, b0 = 10.0, 4.0, 2.0, 2.0

S = 20_000
mu_draws, s2_draws = np.empty(S), np.empty(S)
mu, s2 = 0.0, 1.0                                            # deliberately poor starting values
for t in range(S):
    # mu | sigma^2, y: normal, precision-weighted average of prior mean and sample mean
    prec = 1 / tau0_2 + n / s2
    mu = rng.normal((mu0 / tau0_2 + n * ybar / s2) / prec, np.sqrt(1 / prec))
    # sigma^2 | mu, y: inverse-gamma (draw a gamma for the precision and invert it)
    s2 = 1 / rng.gamma(a0 + n / 2, 1 / (b0 + 0.5 * np.sum((y - mu) ** 2)))
    mu_draws[t], s2_draws[t] = mu, s2

mu_k, sd_k = mu_draws[1000:], np.sqrt(s2_draws[1000:])
for name, d in [("mu", mu_k), ("sigma", sd_k)]:
    lo, hi = np.quantile(d, [0.025, 0.975])
    print(f"{name:6s} mean {d.mean():6.3f}   95% CrI ({lo:6.3f}, {hi:6.3f})   ESS {ess(d):6.0f} of {len(d):,}")
print(f"first draws of mu: {np.round(mu_draws[:4], 2)}  (the chain leaves the bad start at once)")

# Gibbs on a bivariate normal with correlation rho: each step can only move along one axis
print("\n rho   lag-1 autocorrelation of x (theory rho^2)   ESS of x per 10,000 draws")
for rho in [0.5, 0.9, 0.99]:
    x, z, xs = 0.0, 0.0, np.empty(10_000)
    for t in range(10_000):
        x = rng.normal(rho * z, np.sqrt(1 - rho ** 2))       # x | z ~ N(rho z, 1 - rho^2)
        z = rng.normal(rho * x, np.sqrt(1 - rho ** 2))       # z | x ~ N(rho x, 1 - rho^2)
        xs[t] = x
    r1 = np.corrcoef(xs[:-1], xs[1:])[0, 1]
    print(f"{rho:5.2f}   {r1:.3f} ({rho ** 2:.3f})                              {ess(xs):6.0f}")

mu     mean 10.952   95% CrI (10.272, 11.614)   ESS  19000 of 19,000
sigma  mean  1.298   95% CrI ( 0.943,  1.833)   ESS  16934 of 19,000
first draws of mu: [11.09 11.01 11.18 11.08]  (the chain leaves the bad start at once)

 rho   lag-1 autocorrelation of x (theory rho^2)   ESS of x per 10,000 draws
 0.50   0.261 (0.250)                                6081
 0.90   0.799 (0.810)                                1039
 0.99   0.979 (0.980)                                 111

The semi-conjugate Gibbs sampler for the battery data needs no tuning and is extremely efficient: the 19,000 draws of $\mu$ behave like independent ones, and those of $\sigma$ like 16,934. The posterior mean of $\mu$ (10.95) is pulled slightly toward the prior mean of 10, and $\sigma$ is a little smaller than under the noninformative prior of Posterior Summaries (Section 11.2), because the inverse-gamma prior adds information. The bivariate normal shows the weakness: as the correlation rises from 0.5 to 0.99, the lag-1 autocorrelation follows $\rho^2$ and the ESS per 10,000 draws drops from 6,081 to 111.

Key takeaways

• Gibbs sampling updates each parameter by drawing from its full conditional distribution; every draw is accepted and there is nothing to tune.
• Full conditionals are found by keeping only the factors of the joint posterior that involve the parameter; conditional conjugacy makes them standard distributions.
• Data augmentation treats missing data and latent variables as parameters, which makes many complex models easy to sample.
• Gibbs mixes slowly when parameters are strongly correlated; reparameterize, update blocks jointly, or use HMC.

🧪 Practice

1. Derive the full conditional of $\sigma^2$ in the semi-conjugate normal model: write prior × likelihood, drop the factors without $\sigma^2$, and recognize the inverse-gamma kernel.
2. Write a Gibbs sampler for a change-point model: daily counts are Poisson($\lambda_1$) before an unknown day $\tau$ and Poisson($\lambda_2$) after, with Gamma priors on the rates and a uniform prior on $\tau$. (The full conditional of $\tau$ is discrete: evaluate it for every day.) Test it on simulated data with a known change.
3. In the bivariate normal sampler with $\rho = 0.99$, replace the two one-dimensional updates with a single joint draw from the bivariate normal (blocking). What happens to the ESS?
4. Interview: When would you choose Gibbs sampling over a general-purpose sampler such as HMC, and when not? Hint: Gibbs is natural for conditionally conjugate models and for discrete latent variables, which HMC cannot handle directly; it struggles with strong posterior correlations and with non-conjugate full conditionals.

Hamiltonian Monte Carlo

Why it matters. Random-walk Metropolis and Gibbs move by small, undirected steps, so the number of iterations needed to cross the posterior grows quickly with the number of parameters and with their correlations. Hamiltonian Monte Carlo (HMC) uses the gradient of the log-posterior to make long, directed moves that still have a high acceptance probability. It is the default engine of Stan, PyMC and NumPyro, and it is what makes Bayesian models with hundreds or thousands of parameters practical.

The physics analogy. Picture the negative log-posterior $U(\theta) = -\log \tilde p(\theta)$ as a landscape: valleys are high-probability regions. Place a frictionless puck at the current point, give it a random kick (a momentum $r$), and let it slide for a while. It speeds up downhill and slows uphill, so it travels far along valleys and ridges instead of jittering. Where it stops after a fixed time is the proposal. Because the dynamics conserve the total energy, the proposal is almost always accepted.

The mechanics. Augment the parameters with momentum variables $r \sim N(0, M)$ and define the Hamiltonian (total energy)

$$H(\theta, r) = \underbrace{U(\theta)}{\text{potential: } -\log \tilde p(\theta)} + \underbrace{\tfrac12 r^\top M^{-1} r}{\text{kinetic}}.$$

The joint density $\propto e^{-H}$ has the posterior as its $\theta$-marginal. Hamilton's equations, $\frac{d\theta}{dt} = M^{-1}r$ and $\frac{dr}{dt} = -\nabla U(\theta)$, move $(\theta, r)$ along paths of constant $H$. They are simulated with the leapfrog integrator, step size $\varepsilon$, for $L$ steps:

$$r \leftarrow r - \tfrac{\varepsilon}{2}\nabla U(\theta), \qquad \theta \leftarrow \theta + \varepsilon M^{-1} r, \qquad r \leftarrow r - \tfrac{\varepsilon}{2}\nabla U(\theta).$$

Leapfrog is reversible and volume-preserving, which keeps the Metropolis acceptance step simple: accept the endpoint with probability $\min(1, e^{H_{\text{old}} - H_{\text{new}}})$. The small energy error of the discretization is the only reason for rejections.

Tuning.

Setting	Too small	Too large	In practice
Step size $\varepsilon$	many gradients per unit of distance	energy error grows, acceptance collapses; divergences	adapted during warmup to a target acceptance rate (Stan: 0.8)
Path length $\varepsilon L$	behaves like a random walk	the path turns back on itself (U-turn), wasting computation	the No-U-Turn Sampler (NUTS) stops each path automatically when it starts to turn back
Mass matrix $M$			set to an estimate of the inverse posterior covariance during warmup, so all directions have similar scale

Divergences. When the leapfrog path enters a region of high curvature, such as the narrow neck of the "funnel" in hierarchical models (Hierarchical Models, below), the energy error explodes and the path is rejected. Samplers report these divergent transitions. Even a few of them signal that part of the posterior is not being explored and the results may be biased. The fix is usually a reparameterization, not a smaller step.

Scaling. For a $d$-dimensional target, random-walk Metropolis needs a number of density evaluations of order $d$ to produce a nearly independent draw. HMC's step size must shrink like $d^{-1/4}$ to keep the acceptance rate stable, so it needs of order $d^{1/4}$ gradient evaluations per nearly independent draw. The gap grows with dimension. HMC needs gradients, which probabilistic programming languages compute automatically, and it cannot handle discrete parameters directly (they must be summed out, or sampled by Gibbs steps).

import numpy as np

rng = np.random.default_rng(1115)

def ess(x):
    """Effective sample size (same helper as in the Metropolis-Hastings example)."""
    x = x - x.mean()
    n = len(x)
    f = np.fft.rfft(x, 2 * n)
    acf = np.fft.irfft(f * np.conj(f))[:n] / (x @ x)
    neg = np.argmax(acf[1:] < 0) + 1
    return n / (1 + 2 * acf[1:neg].sum())

# Target: a 100-dimensional standard normal, U(q) = -log p(q) = |q|^2 / 2 (+ const)
d = 100
U = lambda q: 0.5 * q @ q
grad_U = lambda q: q

def hmc(n_iter, eps, L):
    q, out, acc = rng.normal(size=d), np.empty((n_iter, d)), 0
    for t in range(n_iter):
        p = rng.normal(size=d)                               # 1. fresh random momentum
        q_new, p_new = q.copy(), p.copy()
        e = eps * rng.uniform(0.8, 1.2)                      # jitter the step size to avoid periodic paths
        p_new -= 0.5 * e * grad_U(q_new)                     # 2. leapfrog: half step for momentum,
        for i in range(L):
            q_new += e * p_new                               #    full steps for position,
            if i < L - 1:
                p_new -= e * grad_U(q_new)                   #    full steps for momentum,
        p_new -= 0.5 * e * grad_U(q_new)                     #    final half step for momentum
        H_old = U(q) + 0.5 * p @ p                           # 3. accept/reject on the change in total energy
        H_new = U(q_new) + 0.5 * p_new @ p_new
        if np.log(rng.random()) < H_old - H_new:
            q, acc = q_new, acc + 1
        out[t] = q
    return out, acc / n_iter

def rwm(n_iter, step):
    q, out, acc = rng.normal(size=d), np.empty((n_iter, d)), 0
    for t in range(n_iter):
        prop = q + step * rng.normal(size=d)
        if np.log(rng.random()) < U(q) - U(prop):
            q, acc = prop, acc + 1
        out[t] = q
    return out, acc / n_iter

n_rwm, n_hmc, L = 50_000, 2_500, 20                          # same budget: 50,000 gradient/density evaluations
for name, (draws, acc), cost in [
        ("random-walk Metropolis", rwm(n_rwm, 2.38 / np.sqrt(d)), n_rwm),
        ("Hamiltonian Monte Carlo", hmc(n_hmc, 0.08, L), n_hmc * L)]:
    e1 = ess(draws[:, 0])                                    # one coordinate
    e2 = ess(np.sum(draws ** 2, axis=1))                     # |q|^2: a global summary of the whole vector
    print(f"{name:24s} acceptance {acc:.2f}   ESS of q1 {e1:6.0f}   ESS of |q|^2 {e2:6.0f}   "
          f"ESS of q1 per 1,000 evaluations {1000 * e1 / cost:5.1f}")

# Tuning: path length eps * L and step size eps both matter
print("\n  eps     L   path length   acceptance   ESS of q1 (2,500 iterations)")
for eps, L in [(0.02, 20), (0.08, 20), (0.3, 20), (0.8, 2)]:
    draws, acc = hmc(n_hmc, eps, L)
    print(f"{eps:5.2f}  {L:4d}   {eps * L:9.2f}       {acc:.2f}        {ess(draws[:, 0]):6.0f}")

random-walk Metropolis   acceptance 0.24   ESS of q1     97   ESS of |q|^2    168   ESS of q1 per 1,000 evaluations   1.9
Hamiltonian Monte Carlo  acceptance 0.99   ESS of q1   2431   ESS of |q|^2   2115   ESS of q1 per 1,000 evaluations  48.6

  eps     L   path length   acceptance   ESS of q1 (2,500 iterations)
 0.02    20        0.40       1.00           117
 0.08    20        1.60       0.99          2447
 0.30    20        6.00       0.96           389
 0.80     2        1.60       0.33           242

On a 100-dimensional normal, with the same budget of 50,000 density or gradient evaluations, optimally tuned random-walk Metropolis accepts 24% of its proposals (the theoretical optimum) but produces only 97 effective draws of a coordinate. HMC accepts 99% of its long trajectories and produces 2,431, a 25-fold improvement per evaluation. The gap widens with dimension. The tuning table shows both failure modes: a path that is too short (0.4) behaves like a random walk (ESS 117), a path that is too long (6.0, about one full oscillation of this target) comes back near its start (ESS 389), and a step size that is too large (0.8) makes the energy error large enough to reject two thirds of the proposals.

Sampler	Needs	Strengths	Weaknesses
Random-walk Metropolis	$\tilde p(\theta)$	simplest; works for any target	slow in high dimensions and with correlations
Gibbs	full conditionals	no tuning; discrete and latent variables	slow with correlations; needs conditional conjugacy
HMC / NUTS	$\tilde p(\theta)$ and $\nabla \log \tilde p(\theta)$	scales to thousands of parameters; diagnoses its own problems (divergences)	continuous parameters only; needs gradients; struggles with multimodality

Key takeaways

• HMC adds momentum variables and simulates Hamiltonian dynamics with the leapfrog integrator, using gradients to make long moves with high acceptance.
• The step size, path length and mass matrix must be tuned; NUTS and warmup adaptation do this automatically.
• Divergent transitions warn that the sampler cannot explore part of the posterior, usually a region of high curvature; reparameterize.
• HMC scales far better with dimension than random-walk methods, but needs continuous parameters and gradients.

🧪 Practice

1. Show that for $U(\theta) = \theta^2/2$ and $M = 1$ the exact Hamiltonian path is $\theta(t) = \theta_0\cos t + r_0 \sin t$. Use this to explain why a path length of about $2\pi$ gives a poor ESS in the tuning table.
2. Run the HMC code with $d = 10$ and $d = 1{,}000$ (adjusting $\varepsilon$ so the acceptance rate stays near 0.9). How does the ESS per evaluation of HMC and of random-walk Metropolis change with $d$?
3. Modify the target to a 2-dimensional normal with correlation 0.99 and compare HMC with an identity mass matrix against HMC with $M^{-1}$ equal to the target covariance.
4. Interview: A Stan model reports "37 divergent transitions after warmup". A colleague suggests ignoring them because R-hat is 1.00. What do you say? Hint: divergences mean the sampler cannot follow the posterior geometry in some region, so the draws can be biased even if the chains agree with each other; increase the target acceptance rate as a check, but the usual fix is reparameterization (non-centered parameterization for hierarchical scales).

Convergence Diagnostics

Why it matters. MCMC output is only guaranteed to represent the posterior "in the long run", and nothing in a single run announces that it has run long enough. A chain may still be drifting from a poor starting point, crawling with high autocorrelation, or stuck in one mode of a multimodal posterior while looking perfectly stable. Convergence diagnostics cannot prove convergence, but they detect the common failures, and running them is a non-negotiable part of every MCMC analysis.

Run several chains. Start 4 or more chains from dispersed starting points. If they have all converged to the posterior, they must agree with each other. Disagreement is the most reliable sign of a problem, and it is invisible with a single chain.

Trace plots. Plot each parameter against the iteration number, with all chains overlaid. Healthy chains look like overlapping "hairy caterpillars"; problems are visible as trends, long flat stretches or chains in separate bands.

  converged: chains overlap,               not converged: chains disagree,
  no trend, rapid mixing                    each stuck in its own region
  theta                                     theta
    | /\/\/\/\/\/\/\/\/\/\/\/\/\/\            | ~~~~~~~~~~~~~~~~~~~~~~~~~ chain 3, 4
    | \/\/\/\/\/\/\/\/\/\/\/\/\/\/            |
    | /\/\/\/\/\/\/\/\/\/\/\/\/\/\            |
    |                                         | ~~~~~~~~~~~~~~~~~~~~~~~~~ chain 1, 2
    +--------------------------- iteration    +--------------------------- iteration

Split R-hat. The potential scale reduction factor (Gelman and Rubin) compares the variance between chains with the variance within chains. Split each chain in half (so that a trend within one chain also shows up as disagreement), giving $m$ sequences of length $n$. With $W$ the average within-sequence variance and $B$ the variance of the sequence means multiplied by $n$:

$$\hat V = \frac{n - 1}{n}W + \frac{1}{n}B, \qquad \hat R = \sqrt{\frac{\hat V}{W}}.$$

If the sequences have converged, all estimate the same distribution and $\hat R \approx 1$. If they disagree, the pooled variance exceeds the within-sequence variance and $\hat R > 1$. Current guidance (Vehtari et al., 2021) is to require $\hat R < 1.01$ for every parameter of interest; modern implementations also rank-normalize the draws so that $\hat R$ works for heavy-tailed posteriors.

Effective sample size and MCSE. Compute the ESS (Metropolis-Hastings Algorithm, above) for every quantity you report. A common rule is ESS > 400 in total (100 per chain for 4 chains) for both the bulk (means, medians) and the tails (5% and 95% quantiles). Then check that the MCSE is small compared with the precision you report.

Diagnostic checklist.

Check	Healthy	Problem indicates
Trace plots	overlapping, stationary, well mixed	drift, stuck chains, separate modes
$\hat R$ (split, rank-normalized)	< 1.01	chains disagree or have not converged
Bulk and tail ESS	> 400	high autocorrelation; run longer or reparameterize
Divergences (HMC)	0	regions the sampler cannot explore; reparameterize
Max tree depth hits (NUTS)	rare	very long paths: poor scaling or strong correlation
Posterior predictive checks	data look typical	the model, not the sampler, is wrong (Section 11.2)

Warmup and thinning. Discard an initial warmup (burn-in) phase, during which the chains move away from their starting points and samplers adapt their tuning parameters; Stan uses half of each chain by default. Thinning (keeping every $k$-th draw) never increases the information in the chain; it is only worth doing to save memory or storage.

import numpy as np

rng = np.random.default_rng(1116)

def ess(x):
    """Effective sample size (same helper as in the Metropolis-Hastings example)."""
    x = x - x.mean()
    n = len(x)
    f = np.fft.rfft(x, 2 * n)
    acf = np.fft.irfft(f * np.conj(f))[:n] / (x @ x)
    neg = np.argmax(acf[1:] < 0) + 1
    return n / (1 + 2 * acf[1:neg].sum())

def split_rhat(chains):
    """Split R-hat: halve each chain, then compare between-chain and within-chain variance."""
    n = chains.shape[1] // 2
    halves = np.vstack([chains[:, :n], chains[:, n:2 * n]])
    W = halves.var(axis=1, ddof=1).mean()                    # average within-chain variance
    B = n * halves.mean(axis=1).var(ddof=1)                  # between-chain variance of the means
    var_plus = (n - 1) / n * W + B / n                       # pooled estimate of the posterior variance
    return np.sqrt(var_plus / W)

def rwm(log_p, start, step, n_iter=5000):
    x, lp, out = start, log_p(start), np.empty(n_iter)
    for t in range(n_iter):
        prop = x + step * rng.normal()
        lp_prop = log_p(prop)
        if np.log(rng.random()) < lp_prop - lp:
            x, lp = prop, lp_prop
        out[t] = x
    return out

normal = lambda x: -0.5 * x ** 2
bimodal = lambda x: np.logaddexp(-0.5 * (x + 5) ** 2, -0.5 * (x - 5) ** 2)   # equal mixture of N(-5,1), N(5,1)
starts = [-8.0, -3.0, 3.0, 8.0]                                               # overdispersed starting points

for name, log_p, step in [("normal, good step", normal, 2.4),
                          ("normal, tiny step", normal, 0.05),
                          ("bimodal, local step", bimodal, 1.0)]:
    chains = np.array([rwm(log_p, s, step)[1000:] for s in starts])          # drop 1,000 warmup iterations
    total_ess = sum(ess(c) for c in chains)
    means = "  ".join(f"{m:+.2f}" for m in chains.mean(axis=1))
    print(f"{name:25s} R-hat {split_rhat(chains):5.2f}   ESS {total_ess:6.0f} of 16,000   chain means {means}")

# The same bimodal run with only the last chain: nothing looks wrong
one = chains[-1:]
print(f"{'bimodal, one chain only':25s} R-hat {split_rhat(one):5.2f}   ESS {ess(one[0]):6.0f} of  4,000   "
      f"chain mean {one.mean():+.2f}  (true mean 0)")

normal, good step         R-hat  1.00   ESS   3623 of 16,000   chain means +0.01  +0.01  -0.03  -0.01
normal, tiny step         R-hat  1.55   ESS     22 of 16,000   chain means -1.08  -0.47  +0.41  +0.93
bimodal, local step       R-hat  5.32   ESS   1741 of 16,000   chain means -5.03  -4.99  +5.04  +4.97
bimodal, one chain only   R-hat  1.00   ESS    391 of  4,000   chain mean +4.97  (true mean 0)

On a standard normal target with a good step size, four chains from widely dispersed starts agree ($\hat R = 1.00$, chain means all within 0.03 of 0). With a tiny step, the chains are still drifting from their starting points after 5,000 iterations: $\hat R = 1.55$ and the ESS is 22. The bimodal case is the dangerous one. Each chain found one of the two modes and stayed there, so two chains report a mean of about $-5$ and two about $+5$; $\hat R = 5.32$ flags the problem immediately. Looked at alone, the last chain passes every check ($\hat R = 1.00$, ESS 391) while reporting a posterior mean of $+4.97$ for a target whose true mean is 0. The per-chain ESS total of 1,741 is also misleading here, which is why ESS should be computed with the multi-chain method of modern software (which incorporates $\hat R$) and never read without $\hat R$.

Key takeaways

• Convergence cannot be proven; diagnostics detect failure. Always run several chains from dispersed starting points.
• Split $\hat R$ compares between-chain and within-chain variance; require $\hat R < 1.01$.
• Check bulk and tail ESS (> 400) and the MCSE of every reported quantity; with HMC, require zero divergences.
• A single chain stuck in one mode can look perfectly healthy; multimodality calls for multiple chains, better parameterization or specialized samplers.

🧪 Practice

1. Four chains have within-chain variances 1.0, 1.1, 0.9, 1.0 and means 0.0, 0.1, -0.1, 2.0, each of length 1,000 (ignore splitting). Compute $\hat R$. Which chain would you investigate?
2. Rerun the "normal, tiny step" case with 20,000, 100,000 and 500,000 iterations. How long until $\hat R < 1.01$ and ESS > 400?
3. Add a trend to a well-mixed chain (add $0.0005 \times$ iteration to each draw). Show that the split $\hat R$ detects it with a single chain, but the unsplit version does not.
4. Interview: A junior analyst ran one MCMC chain of 100,000 iterations, threw away the first half, thinned by 10, and reports results because "the trace plot looks fine". What would you ask them to change? Hint: multiple chains from dispersed starts with $\hat R$; ESS and MCSE rather than raw length; thinning loses information; check divergences and posterior predictive fit.

Hierarchical Models

Why it matters. Data often come in groups: students in schools, patients in hospitals, customers in stores, measurements from many sensors. Each group has its own parameter, but the groups are similar. Estimating each group separately is noisy for small groups; pooling all groups ignores real differences. Hierarchical (multilevel) models let the data decide how much to share, by treating the group parameters as draws from a common distribution whose parameters are also estimated. They are one of the most useful ideas in applied statistics, and the Bayesian framework handles them naturally.

Three ways to estimate group effects.

Approach	Model	Problem
No pooling	a separate $\theta_j$ per group, estimated from its own data	small groups give extreme, noisy estimates
Complete pooling	one common $\theta$ for all groups	ignores real differences between groups
Partial pooling (hierarchical)	$\theta_j \sim N(\mu, \tau^2)$, with $\mu$ and $\tau$ estimated	none of the above: interpolates according to the data

The model. For group estimates $y_j$ with known standard errors $\sigma_j$ (for example, effect estimates from $J$ separate experiments):

$$\begin{aligned} y_j \mid \theta_j &\sim N(\theta_j, \sigma_j^2) && \text{(data level: group } j\text{'s estimate)} \ \theta_j \mid \mu, \tau &\sim N(\mu, \tau^2) && \text{(group level: groups vary around a common mean)} \ \mu, \tau &\sim p(\mu, \tau) && \text{(hyperprior)} \end{aligned}$$

                 (mu, tau)                    hyperparameters: population of groups
              /   /    \    \
        theta_1 theta_2 ... theta_J           group parameters
           |       |           |
          y_1     y_2   ...   y_J             observed data

Shrinkage. Given $\mu$ and $\tau$, the posterior mean of each group effect is a precision-weighted average of its own estimate and the overall mean:

$$E[\theta_j \mid y_j, \mu, \tau] = (1 - B_j),y_j + B_j,\mu, \qquad B_j = \frac{\sigma_j^2}{\sigma_j^2 + \tau^2}.$$

Noisy groups (large $\sigma_j$) are shrunk more. If the groups really differ (large $\tau$), shrinkage is weak. If they are alike ($\tau \to 0$), everything collapses to the common mean. The crucial difference from fixing a prior by hand is that $\tau$ is estimated from the spread of the groups: the data decide how much to pool. This is the Bayesian counterpart of the James-Stein estimator mentioned in Section 8.1, which shows that with three or more groups, shrinking toward a common mean reduces the total squared error compared with estimating each group separately.

Connection to mixed models. The frequentist counterpart is the linear mixed model with random effects (Section 9.4), fitted by (restricted) maximum likelihood. It gives the same shrinkage, with $\tau$ estimated by a point value. The Bayesian version also propagates the uncertainty about $\tau$, which matters when there are few groups and $\tau$ is poorly determined.

Example: eight schools. A classic dataset (Rubin, 1981): separate randomized experiments in eight schools estimated the effect of SAT coaching programs, each with its own standard error. School A's estimate is 28 points, but with a standard error of 15. Should we believe that A's program is much better than the others? With flat priors on $\mu$ and $\tau$, the posterior can be computed without MCMC. The marginal posterior of $\tau$ is evaluated on a grid (with $\mu$ and the $\theta_j$ integrated out analytically), and then $\tau$, $\mu$ and the $\theta_j$ are drawn in sequence.

import numpy as np

rng = np.random.default_rng(1117)

# Eight schools (Rubin, 1981): estimated coaching effect on SAT scores and its standard error
y = np.array([28, 8, -3, 7, -1, 1, 18, 12], dtype=float)
se = np.array([15, 10, 16, 11, 9, 11, 10, 18], dtype=float)
names = "ABCDEFGH"

# Model: y_j ~ N(theta_j, se_j^2),  theta_j ~ N(mu, tau^2),  flat priors on mu and on tau > 0
# 1. Marginal posterior of tau on a grid (mu and the theta_j integrated out analytically)
tau = np.linspace(0.01, 40, 4000)
log_p = np.empty_like(tau)
for i, t in enumerate(tau):
    v = se ** 2 + t ** 2                                     # marginal variance of y_j given mu and tau
    V_mu = 1 / np.sum(1 / v)
    mu_hat = V_mu * np.sum(y / v)
    log_p[i] = 0.5 * np.log(V_mu) - 0.5 * np.sum(np.log(v)) - 0.5 * np.sum((y - mu_hat) ** 2 / v)
p_tau = np.exp(log_p - log_p.max())
p_tau /= p_tau.sum()

# 2. Draw tau, then mu | tau, y, then each theta_j | mu, tau, y
S = 50_000
t = rng.choice(tau, S, p=p_tau)
v = se ** 2 + t[:, None] ** 2
V_mu = 1 / np.sum(1 / v, axis=1)
mu = rng.normal(V_mu * np.sum(y / v, axis=1), np.sqrt(V_mu))
prec = 1 / se ** 2 + 1 / t[:, None] ** 2
theta = rng.normal((y / se ** 2 + mu[:, None] / t[:, None] ** 2) / prec, np.sqrt(1 / prec))

pooled = np.sum(y / se ** 2) / np.sum(1 / se ** 2)          # complete pooling: one common effect
print(f"tau: posterior median {np.median(t):.1f}, 95% upper bound {np.quantile(t, 0.95):.1f};  "
      f"mu: posterior mean {mu.mean():.1f};  complete-pooling estimate {pooled:.1f}\n")
print("school   no pooling (y +- se)   partial pooling: mean   95% CrI          weight on mu")
for j in range(8):
    lo, hi = np.quantile(theta[:, j], [0.025, 0.975])
    shrink = np.mean(se[j] ** 2 / (se[j] ** 2 + t ** 2))     # shrinkage factor B_j, averaged over tau
    print(f"   {names[j]}      {y[j]:5.0f} +- {se[j]:2.0f}                {theta[:, j].mean():5.1f}   "
          f"({lo:6.1f}, {hi:5.1f})        {shrink:.2f}")
print(f"\nP(school A's effect > 28 | data) = {np.mean(theta[:, 0] > 28):.3f}")

tau: posterior median 5.2, 95% upper bound 17.2;  mu: posterior mean 7.9;  complete-pooling estimate 7.7

school   no pooling (y +- se)   partial pooling: mean   95% CrI          weight on mu
   A         28 +- 15                 11.4   (  -2.0,  31.4)        0.83
   B          8 +- 10                  7.9   (  -4.6,  20.8)        0.72
   C         -3 +- 16                  6.1   ( -11.4,  20.3)        0.84
   D          7 +- 11                  7.6   (  -5.8,  20.8)        0.75
   E         -1 +-  9                  5.2   (  -8.7,  16.4)        0.69
   F          1 +- 11                  6.1   (  -8.4,  18.6)        0.75
   G         18 +- 10                 10.6   (  -1.2,  25.9)        0.72
   H         12 +- 18                  8.5   (  -7.0,  25.7)        0.86

P(school A's effect > 28 | data) = 0.043

The data are consistent with modest variation between schools: the posterior median of $\tau$ is 5.2 points, far smaller than most standard errors. As a result, every school is shrunk most of the way (weights of 0.69 to 0.86 on $\mu$) toward the overall mean of about 8. School A's estimate drops from 28 to 11.4, and there is only a 4% posterior probability that its true effect is as large as its raw estimate. School C's $-3$ rises to 6.1. The partial-pooling estimates lie between the no-pooling estimates (each school alone, wildly different) and the complete-pooling estimate (7.7 for all), and the noisiest schools (A, C, H) are shrunk the most.

Computation and the funnel. In larger hierarchical models with MCMC, the posterior of $(\theta_j, \tau)$ has a funnel shape: when $\tau$ is small, the $\theta_j$ are squeezed tightly around $\mu$, and HMC produces divergences there. The standard fix is the non-centered parameterization: sample $z_j \sim N(0, 1)$ and set $\theta_j = \mu + \tau z_j$, which removes the dependence between the group effects and $\tau$ in the prior.

Where hierarchical models are used. Small-area estimation (unemployment by county), meta-analysis (combining studies), A/B tests with many segments or many simultaneous experiments, sports ratings, recommender systems (user and item effects), repeated measures and longitudinal data, and multilevel regression with poststratification (MRP) for polling.

Key takeaways

• Hierarchical models treat group parameters as draws from a common distribution with estimated hyperparameters, giving partial pooling between no pooling and complete pooling.
• Each group estimate is shrunk toward the overall mean by $B_j = \sigma_j^2/(\sigma_j^2 + \tau^2)$: noisy groups and similar groups are shrunk more.
• The between-group SD $\tau$ is learned from the data, so the data decide how much to pool; Bayesian inference also propagates its uncertainty.
• With MCMC, use the non-centered parameterization to avoid the funnel geometry and divergences when $\tau$ may be small.

🧪 Practice

1. With $\mu = 8$ and $\tau = 5$, compute the shrinkage factor and the conditional posterior mean for a school with $y_j = 28$, $\sigma_j = 15$, and for one with $y_j = 28$, $\sigma_j = 3$. Explain the difference.
2. Modify the code to use a half-normal prior with scale 10 on $\tau$ (multiply the grid posterior by the prior density). How do the posterior median of $\tau$ and the estimate for school A change? Repeat with scale 1.
3. Simulate 50 "stores" with true conversion rates drawn from Beta(20, 380) and between 20 and 2,000 visitors each. Compare the mean squared error of the raw rates $k_j/n_j$ with that of a hierarchical estimate (fit $\alpha, \beta$ of the Beta population by maximizing the beta-binomial marginal likelihood, then use the posterior means). Which stores benefit most?
4. Interview: You run 40 A/B tests on different country segments and the largest observed lift is +12% in a small country. Your manager wants to launch only there. How would a hierarchical model change the analysis? Hint: the largest of many noisy estimates is biased upward (winner's curse); a hierarchical model shrinks small segments toward the overall effect according to the estimated between-segment variation, giving a more realistic estimate and interval for that country.

<a id="12-advanced-and-applied-topics"></a>