47. Credible intervals

A posterior distribution contains everything you know about \theta. To communicate it you usually want a single interval — and unlike a confidence interval, this one means what people always thought confidence intervals meant.

The definition

A 95% credible interval is any region C with

P(\theta \in C \mid \text{data}) = 0.95

That is a direct probability statement about \theta, made possible because Bayesian inference treats \theta as a random variable.

Compare carefully with §4:

	Confidence interval	Credible interval
\theta is	fixed, unknown	random
The interval is	random	fixed (given the data)
95% refers to	the long-run procedure	this actual interval
Legitimate claim	"95% of such intervals cover \theta"	"P(\theta \in C) = 0.95"
Needs a prior	no	yes

The frequentist statement is about what would happen in repeated sampling. The Bayesian statement is about this dataset. Neither is wrong; they answer different questions.

Two kinds of interval

The definition doesn't pick a unique region — many intervals contain 95% of the posterior. Two conventions dominate.

Equal-tailed interval. Cut 2.5% from each tail: the interval runs from the 2.5th to the 97.5th percentile.

• Easy to compute (just quantiles), and invariant under monotone transformations: the interval for \log\theta is the log of the interval for \theta.
• But for a skewed posterior it can be wider than necessary, and can even include points of lower density than points it excludes.

Highest posterior density (HPD) interval. The shortest region containing 95% of the mass. Equivalently, every point inside has higher posterior density than every point outside.

• Always the narrowest, and always contains the posterior mode.
• Not transformation-invariant, and can be disjoint if the posterior is multimodal — which is arguably a feature: a bimodal posterior genuinely shouldn't be summarised by one connected range.

For a symmetric posterior the two coincide. They differ when it's skewed, which is exactly when the choice matters.

Point estimates from a posterior

Which summary is "best" depends on the loss function you'd pay for being wrong — this is where Bayesian inference connects to decision theory:

Loss function	Optimal estimate
Squared error (\hat\theta - \theta)^2	posterior mean
Absolute error \lvert\hat\theta - \theta\rvert	posterior median
0–1 loss	posterior mode (MAP)

Squared error is the reason the mean is the default, not any deep principle. If overestimating costs more than underestimating, an asymmetric loss gives a different optimum entirely — and the posterior lets you compute it.

The MAP estimate is the closest Bayesian analogue of the MLE. In fact with a flat prior they coincide exactly, since maximising the posterior reduces to maximising the likelihood. Ridge regression (§6) is the MAP estimate under a Normal prior on the coefficients, and lasso is the MAP under a Laplace prior — the penalties from that lesson are priors in disguise.

Beyond intervals

The real advantage of a posterior is that you're not restricted to intervals at all. Any question about \theta has an answer:

P(\theta > 0 \mid \text{data}), \qquad P(0.4 < \theta < 0.6 \mid \text{data}), \qquad P(\theta_B > \theta_A \mid \text{data})

And you can propagate uncertainty through any function. Want a distribution for e^\theta, or for \theta_1/\theta_2? Draw from the posterior, transform each draw, and look at the result. No delta method, no asymptotics — this works in small samples and for arbitrary transformations, which is a genuine advantage over §3's machinery.

Worked example

Posterior p \mid \text{data} \sim \text{Beta}(9, 5) (from the coin example two lessons ago).

Equal-tailed 95% interval: the 2.5th and 97.5th percentiles of Beta(9,5):

(0.386,\; 0.861)

HPD 95% interval: (0.401, 0.874) — shifted right and slightly narrower, because the posterior is left-skewed and the HPD trims more from the thin left tail.

Point estimates:

\text{mean} = \frac{9}{14} = 0.643, \qquad \text{median} \approx 0.649, \qquad \text{mode} = \frac{9-1}{14-2} = 0.667

All three differ, and the ordering (mean < median < mode) is the left-skew signature from §0, running the opposite way to the salary data.

Direct probability statements, which no confidence interval could give:

P(p > 0.5 \mid \text{data}) \approx 0.867, \qquad P(p > 0.8 \mid \text{data}) \approx 0.099

So: 87% confident the coin favours heads, but only 10% confident it's heavily biased. That's a far more useful summary for a decision than "reject H_0 at \alpha = 0.05".

Doing it in Python

Both interval types, and the point estimates:

import numpy as np
from scipy.stats import beta
from scipy.optimize import minimize_scalar

post = beta(9, 5)

lo, hi = post.ppf([0.025, 0.975])
print(f"equal-tailed 95%: ({lo:.4f}, {hi:.4f})   width {hi-lo:.4f}")

def hpd(dist, mass=0.95):
    """Shortest interval containing `mass` of the distribution."""
    def width(lower_tail):
        a = dist.ppf(lower_tail)
        b = dist.ppf(lower_tail + mass)
        return b - a
    res = minimize_scalar(width, bounds=(1e-9, 1 - mass - 1e-9), method="bounded")
    return dist.ppf(res.x), dist.ppf(res.x + mass)

h_lo, h_hi = hpd(post)
print(f"HPD 95%         : ({h_lo:.4f}, {h_hi:.4f})   width {h_hi-h_lo:.4f}")

print(f"\nmean   {post.mean():.4f}")
print(f"median {post.median():.4f}")
print(f"mode   {(9-1)/(9+5-2):.4f}")

Questions a confidence interval can't answer:

from scipy.stats import beta

post = beta(9, 5)

for threshold in (0.4, 0.5, 0.6, 0.7, 0.8):
    print(f"P(p > {threshold}) = {1 - post.cdf(threshold):.4f}")

print(f"\nP(0.4 < p < 0.6) = {post.cdf(0.6) - post.cdf(0.4):.4f}")
print("\nThese are probabilities about the PARAMETER, which is only")
print("meaningful because the Bayesian framework treats it as random.")

Propagating uncertainty through an arbitrary function — no asymptotics needed:

import numpy as np
from scipy.stats import beta

rng = np.random.default_rng(0)
draws = beta(9, 5).rvs(200_000, random_state=rng)

# The odds, p/(1-p) -- a nonlinear, skewed transformation
odds = draws / (1 - draws)

print(f"posterior mean of p     : {draws.mean():.4f}")
print(f"posterior mean of odds  : {odds.mean():.4f}")
print(f"odds at the mean of p   : {draws.mean()/(1-draws.mean()):.4f}  <- NOT the same")
print(f"\n95% interval for odds: ({np.quantile(odds, 0.025):.4f}, "
      f"{np.quantile(odds, 0.975):.4f})")
print(f"P(odds > 2)          : {(odds > 2).mean():.4f}")

That mismatch is Jensen's inequality again (§4) — E[g(\theta)] \ne g(E[\theta]) for nonlinear g. Sampling handles it automatically.

When equal-tailed and HPD really diverge:

import numpy as np
from scipy.stats import beta, gamma
from scipy.optimize import minimize_scalar

def hpd(dist, mass=0.95):
    def width(t):
        return dist.ppf(t + mass) - dist.ppf(t)
    res = minimize_scalar(width, bounds=(1e-9, 1 - mass - 1e-9), method="bounded")
    return dist.ppf(res.x), dist.ppf(res.x + mass)

cases = {
    "symmetric  Beta(10,10)": beta(10, 10),
    "mild skew  Beta(9,5)":   beta(9, 5),
    "heavy skew Beta(1.5,12)": beta(1.5, 12),
    "very skew  Gamma(1.2)":  gamma(1.2),
}
print(f"{'posterior':>24} {'equal-tailed width':>20} {'HPD width':>12} {'saving':>8}")
for name, d in cases.items():
    lo, hi = d.ppf([0.025, 0.975])
    a, b = hpd(d)
    print(f"{name:>24} {hi-lo:>20.4f} {b-a:>12.4f} {100*(1-(b-a)/(hi-lo)):>7.1f}%")

print("\nThe more skewed the posterior, the more the HPD saves.")

Coverage: the credible interval really does contain the parameter 95% of the time when the prior is right:

import numpy as np
from scipy.stats import beta

rng = np.random.default_rng(1)
trials, n = 40_000, 20
a0, b0 = 3, 3                       # the prior we will actually sample from

# Draw the TRUE p from the prior, then data from it -- the Bayesian setup
true_p = rng.beta(a0, b0, trials)
k = rng.binomial(n, true_p)

lo = beta.ppf(0.025, a0 + k, b0 + n - k)
hi = beta.ppf(0.975, a0 + k, b0 + n - k)
print(f"coverage when the prior is correct: {((lo < true_p) & (true_p < hi)).mean():.4f}")

# Now the same intervals when the truth is always 0.5 -- prior is wrong-ish
true_fixed = np.full(trials, 0.9)
k2 = rng.binomial(n, true_fixed)
lo2 = beta.ppf(0.025, a0 + k2, b0 + n - k2)
hi2 = beta.ppf(0.975, a0 + k2, b0 + n - k2)
print(f"coverage when p is always 0.9      : {((lo2 < 0.9) & (0.9 < hi2)).mean():.4f}")
print("\nBayesian intervals are calibrated ON AVERAGE OVER THE PRIOR.")
print("For a particular parameter value, coverage can be worse -- that's the")
print("price of using prior information, and why the prior should be defensible.")

Your turn

1. A 95% credible interval is (2, 8). State what that means.

2. When do equal-tailed and HPD intervals coincide?

3. Your posterior is N(10, 2^2). Find the 95% credible interval and P(\theta > 12).

Solutions

1. Given the prior and the observed data, there is a 95% probability that \theta lies between 2 and 8.

That is a statement about \theta itself, and it is legitimate here precisely because the Bayesian framework gives \theta a distribution.

Contrast with the frequentist 95% CI (2,8), where the only valid statement is about the procedure: 95% of intervals constructed this way would contain the true \theta. For that specific interval, \theta is either in it or not, with no probability attached (§4).

The Bayesian claim comes with a condition attached — given the prior. Someone with a different prior gets a different interval. Whether that's a bug or a feature is the argument of the next lesson but one.

2. When the posterior is symmetric and unimodal.

Under symmetry, cutting 2.5% from each tail lands at points of equal density, which is precisely the HPD condition. Any other 95% region would have to extend further into one tail than it gains from the other, making it wider.

The standard case is a Normal posterior, where both give \mu \pm 1.96\sigma.

They diverge when the posterior is skewed: the equal-tailed interval keeps cutting 2.5% from a long thin tail, while the HPD trims that tail harder and gains a narrower interval — as the code above shows, saving over 10% of the width for strongly skewed posteriors.

3. For a Normal posterior the interval is symmetric:

10 \pm 1.96(2) = (6.08,\; 13.92)

And because the Normal is symmetric and unimodal, this is also the HPD interval.

For the probability:

P(\theta > 12) = P\left(Z > \frac{12 - 10}{2}\right) = P(Z > 1) = 0.1587

About a 16% probability that \theta exceeds 12 — a statement you simply cannot make from a frequentist confidence interval, however similar the numbers look.

Check yourself in code

Compute both interval types for a Beta(9,5) posterior, confirm the HPD is narrower, and answer a direct probability question.

Print exactly this:

equal-tailed (0.3857, 0.8614)
HPD (0.4013, 0.8737)
HPD is narrower: True
P(p > 0.5) 0.8666

Round every endpoint to 4 decimal places and the probability to 4. Find the HPD by minimising the interval width over the position of the lower tail with scipy.optimize.minimize_scalar.

import numpy as np
from scipy.stats import beta
from scipy.optimize import minimize_scalar

post = beta(9, 5)

lo, hi = post.ppf([0.025, 0.975])
print(f"equal-tailed ({lo:.4f}, {hi:.4f})")

# Find the shortest 95% interval by minimising ppf(t+0.95) - ppf(t) over t,
# print it, say whether it is narrower, and give P(p > 0.5).

import numpy as np
from scipy.stats import beta
from scipy.optimize import minimize_scalar

post = beta(9, 5)

lo, hi = post.ppf([0.025, 0.975])
print(f"equal-tailed ({lo:.4f}, {hi:.4f})")

res = minimize_scalar(lambda t: post.ppf(t + 0.95) - post.ppf(t),
                      bounds=(1e-9, 0.05 - 1e-9), method="bounded")
h_lo, h_hi = post.ppf(res.x), post.ppf(res.x + 0.95)
print(f"HPD ({h_lo:.4f}, {h_hi:.4f})")
print("HPD is narrower:", bool((h_hi - h_lo) < (hi - lo)))
print("P(p > 0.5)", round(1 - post.cdf(0.5), 4))

---

A credible interval says what everyone wants a confidence interval to say: the probability that the parameter lies in this range. Equal-tailed intervals are simple and transformation-invariant; HPD intervals are the shortest and always contain the mode. And because you hold the whole posterior, you're never restricted to intervals — any probability question, and any transformation, is a sampling exercise.

Next: the two frameworks side by side, and when the difference actually matters.