11 - Bayesian Estimation: Prior, Posterior, and MAP

You showed a new button to 20 visitors and 8 clicked. The observed fraction is 0.40, but what you really want to know is something else: how far can you trust that number, and what changes if similar buttons converted at 10% in the past? Frequentist statistics answers with an estimator and a confidence interval; Bayesian statistics answers with a probability distribution over the true value, built by combining what you knew before with what you measured. In this article, the eleventh in the Probability and Statistics series, you compute that distribution in full on a case small enough to work by hand, compare it with the maximum likelihood estimate (MLE, the value that makes the data most probable, seen in the article on estimators), and then let a sampler compute it when no closed form exists. Every result is checked with an exact calculation and with a fixed-seed simulation.

From Bayes' theorem to the posterior

Call θ (theta) the parameter you do not know, here the true probability that a visitor clicks, and D the data. The prior p(θ) describes how plausible each value of θ is before you look at the data. The likelihood p(D | θ) says how probable those data would be if θ were exactly that value. Bayes' theorem combines them into the posterior: p(θ | D) = p(D | θ) · p(θ) / p(D). The denominator p(D), called the evidence, is the sum (an integral) of the numerator over all possible θ and only serves to make the area equal 1: to get the shape of the posterior it is therefore enough to multiply likelihood and prior.

With k successes in n independent trials the likelihood is θ^k (1 − θ)^n−k. If you choose as prior a Beta(a, b) distribution, that is a density proportional to θ^a−1 (1 − θ)^b−1 on the interval from 0 to 1, the product has the same shape with the exponents added: the posterior is Beta(a + k, b + n − k). A prior with this property, whereby the posterior stays in the same family, is called conjugate. Read a and b as «virtual» successes and failures observed before the test: Beta(1, 1) is the uniform density, equivalent to having seen nothing; Beta(10, 10) is equivalent to twenty past trials with half of them successes.

A posterior computed in full: 8 successes in 20

The program below starts from Beta(1, 1), arrives at Beta(9, 13), and computes the point estimates and the intervals. First, though, it runs a check that does not depend on the formula: it builds the posterior by brute force, multiplying likelihood and prior on a grid of a thousand values of θ and normalising. If the conjugate formula were wrong, the two curves would not coincide; instead they differ by less than 10^-14. The HPD (highest posterior density) is the shortest interval that contains 95% of the probability; the hpd function finds it by trying all the intervals with that mass and keeping the narrowest. The equal-tail interval (ETI) instead cuts off 2.5% on each side. The Wilson interval is the frequentist interval for a proportion, computed with Wilson's formula (1927).

import numpy as np
from scipy import stats

n, k = 20, 8                      # 8 successes in 20 trials
a0, b0 = 1, 1                     # uniform prior Beta(1, 1)
a1, b1 = a0 + k, b0 + (n - k)     # conjugate posterior: Beta(9, 13)
post = stats.beta(a1, b1)

# Check: posterior computed by brute force on a grid (likelihood x prior, then normalise)
theta = np.linspace(0.0005, 0.9995, 1000)
dx = theta[1] - theta[0]
raw = theta**k * (1 - theta)**(n - k) * stats.beta.pdf(theta, a0, b0)
grid = raw / (raw.sum() * dx)
print(f"max |grid - Beta(9,13)| = {np.abs(grid - post.pdf(theta)).max():.1e}")

# Point estimates and intervals: closed form
mle = k / n
mean = post.mean()
median = post.ppf(0.5)
mode = (a1 - 1) / (a1 + b1 - 2)           # MAP with a uniform prior = MLE
eti = (post.ppf(0.025), post.ppf(0.975))

def hpd(dist, mass=0.95, steps=20000):
    # among all intervals with probability 'mass', take the shortest
    p_low = np.linspace(0, 1 - mass, steps)
    lo = dist.ppf(p_low)
    hi = dist.ppf(p_low + mass)
    j = np.argmin(hi - lo)
    return lo[j], hi[j]

h = hpd(post)
z = stats.norm.ppf(0.975)
center = (k + z**2 / 2) / (n + z**2)
half = z * np.sqrt(k * (n - k) / n + z**2 / 4) / (n + z**2)
print(f"MLE {mle:.4f}  MAP {mode:.4f}  median {median:.4f}  mean {mean:.4f}  sd {post.std():.4f}")
print(f"ETI 95%  [{eti[0]:.4f}, {eti[1]:.4f}]  width {eti[1]-eti[0]:.4f}")
print(f"HPD 95%  [{h[0]:.4f}, {h[1]:.4f}]  width {h[1]-h[0]:.4f}")
print(f"Wilson 95% [{center-half:.4f}, {center+half:.4f}]")

# Fixed-seed simulation: it must confirm the exact result within the Monte Carlo error
rng = np.random.default_rng(2026)
N = 200_000
x = rng.beta(a1, b1, N)
print(f"simulation ({N} draws): mean {x.mean():.4f} +- {x.std()/np.sqrt(N):.4f}, "
      f"q2.5% {np.quantile(x, .025):.4f}, q97.5% {np.quantile(x, .975):.4f}")

# Predictive: probability of at least 5 successes in 10 new trials
pp = stats.betabinom(10, a1, b1)
plug = stats.binom(10, mle)
y = rng.binomial(10, x)                    # one scenario for each posterior draw
print(f"P(Y>=5) predictive {1 - pp.cdf(4):.4f}  simulated {np.mean(y >= 5):.4f}  plug-in MLE {1 - plug.cdf(4):.4f}")
print(f"variance of Y: predictive {pp.var():.3f}  plug-in {plug.var():.3f}")

# Output:
# max |grid - Beta(9,13)| = 9.1e-15
# MLE 0.4000  MAP 0.4000  median 0.4063  mean 0.4091  sd 0.1025
# ETI 95%  [0.2182, 0.6156]  width 0.3974
# HPD 95%  [0.2129, 0.6097]  width 0.3968
# Wilson 95% [0.2188, 0.6134]
# simulation (200000 draws): mean 0.4093 +- 0.0002, q2.5% 0.2180, q97.5% 0.6156
# P(Y>=5) predictive 0.4033  simulated 0.4043  plug-in MLE 0.3669
# variance of Y: predictive 3.363  plug-in 2.400

The MAP (maximum a posteriori, the highest point of the posterior) is 0.4000 and coincides with the MLE: with a uniform prior the posterior is proportional to the likelihood, and the maximum does not move. Mean (0.4091) and median (0.4063) are higher because the Beta(9, 13) has the longer right tail. The equal-tail interval [0.2182, 0.6156] and the Wilson interval [0.2188, 0.6134] are almost identical: with 20 trials and a flat prior the two paradigms say nearly the same thing. The simulation with 200,000 draws confirms the exact calculation: mean 0.4093 against 0.4091 (standard error 0.0002) and quantiles equal to the third digit.

The last part of the program computes the predictive: the probability of getting at least 5 clicks from another 10 visitors. It is 0.4033 when you integrate over the uncertainty in θ, against 0.3669 when you «trust» 0.40 as if it were the true value. The variance of the number of clicks is 3.36 against 2.40: the plug-in underestimates the spread because it ignores that θ is not known.

Which point estimate to choose

The posterior is the complete answer, but often you need a single number. Which one depends on what a mistake costs you, that is, on the loss function:

• squared loss (the error weighs as its square): the best estimate is the posterior mean;
• absolute loss (the error weighs in proportion): the median;
• 0-1 loss (all that matters is guessing exactly): the MAP, the mode.

There is therefore no single «Bayesian» estimate: there are three answers to three different questions, and they coincide only when the posterior is symmetric and unimodal. For the Beta(9, 13) they differ by a few thousandths, but for a strongly asymmetric posterior, or one with two peaks, the choice changes the number you report. One property of the MAP is worth remembering: it depends on how you parametrise the problem. If you estimate θ or its log-odds you get two modes that do not correspond to each other, whereas the median stays the same under every monotone transformation (the mean does not).

MAP and regularisation: the same formula seen from two sides

In a linear regression y = Xβ + noise, with Gaussian noise of standard deviation σ, looking for the MAP when each coefficient has a Gaussian prior N(0, τ²) is equivalent to minimising the sum of squares plus a penalty on the coefficients, that is, ridge regression, with parameter λ = σ² / τ². The prior explains why the penalty works: it says that large coefficients are a priori implausible. The program verifies the equality numerically. It computes the closed-form ridge solution and then looks for the maximum of the log-posterior with a generic optimiser (BFGS) that knows nothing of that formula. If the theory is right, the two results coincide. For regularisation from the machine learning side see Overfitting and regularization.

import numpy as np
from scipy import optimize

rng = np.random.default_rng(7)
n, p = 40, 3
X = rng.standard_normal((n, p))
true_beta = np.array([2.0, -1.0, 0.0])
sigma, tau = 1.5, 1.0                       # data noise and standard deviation of the prior N(0, tau^2)
y = X @ true_beta + sigma * rng.standard_normal(n)

lam = sigma**2 / tau**2                     # the ridge <-> prior link
ridge = np.linalg.solve(X.T @ X + lam * np.eye(p), X.T @ y)

def neg_log_posterior(b):                   # up to constants
    return 0.5 * np.sum((y - X @ b) ** 2) / sigma**2 + 0.5 * np.sum(b**2) / tau**2

map_num = optimize.minimize(neg_log_posterior, np.zeros(p), method="BFGS", options={"gtol": 1e-10}).x
ols = np.linalg.lstsq(X, y, rcond=None)[0]
print("OLS (MLE)  ", np.round(ols, 4))
print("ridge      ", np.round(ridge, 4))
print("numeric MAP", np.round(map_num, 4), " max gap", f"{np.abs(ridge - map_num).max():.1e}")

# Output:
# OLS (MLE)   [ 2.0581 -0.9495  0.6028]
# ridge       [ 1.9052 -0.896   0.5908]
# numeric MAP [ 1.9052 -0.896   0.5908]  max gap 1.2e-08

The two vectors coincide up to the fourth digit shown and the maximum gap is 1.2·10^-8, that is, the optimiser's tolerance. The effect of the prior is visible too: all three coefficients are pulled towards zero compared with least squares (2.0581 becomes 1.9052, 0.6028 becomes 0.5908). Choosing τ is choosing the strength of the penalty: a tight prior (small τ) means a large λ and more strongly shrunk coefficients. This equivalence holds for the MAP, not for the whole posterior: ridge returns only the highest point and loses the uncertainty around it.

What a credible interval really promises

A 95% credible interval is a set of values of θ to which the posterior assigns probability 0.95. The reading «θ is in there with probability 0.95» is legitimate, but it is conditional on the chosen prior: the guarantee does not cover all possible values of θ in the same way. The test below shows it. In the first case θ is drawn from the Beta(2, 2) prior itself and then the data are simulated: the interval covers the true value in 95.0% of cases, as promised. In the next three cases θ is fixed and the prior is wrong for that θ: coverage becomes 92.4% if θ = 0.05, 96.8% if θ = 0.20, 95.7% if θ = 0.50. The simulation error is about 0.001, so the departures from 95% are real.

import numpy as np
from scipy import stats

rng = np.random.default_rng(99)
n = 20

# --- 1) What a credible interval promises -----------------------------------
a0, b0 = 2, 2                               # prior Beta(2, 2)
def coverage(thetas, k):
    lo = stats.beta.ppf(0.025, a0 + k, b0 + n - k)
    hi = stats.beta.ppf(0.975, a0 + k, b0 + n - k)
    return np.mean((lo <= thetas) & (thetas <= hi))

R = 100_000
th = rng.beta(a0, b0, R)                    # theta drawn from the prior itself
print(f"coverage if theta ~ prior   : {coverage(th, rng.binomial(n, th)):.4f}  (expected 0.95, error ~0.0007)")
for fixed_t in (0.05, 0.20, 0.50):
    th = np.full(R, fixed_t)
    print(f"coverage if theta = {fixed_t:.2f}   : {coverage(th, rng.binomial(n, th)):.4f}")

# --- 2) How much the prior weighs: n=10, k=7 ---------------------------------
print("prior        mean    ETI 95%")
for a, b in [(1, 1), (0.5, 0.5), (10, 10), (1, 10)]:
    d = stats.beta(a + 7, b + 3)
    print(f"Beta({a},{b})".ljust(12), f"{d.mean():.3f}  [{d.ppf(.025):.3f}, {d.ppf(.975):.3f}]")

# Output:
# coverage if theta ~ prior   : 0.9502  (expected 0.95, error ~0.0007)
# coverage if theta = 0.05   : 0.9241
# coverage if theta = 0.20   : 0.9676
# coverage if theta = 0.50   : 0.9571
# prior        mean    ETI 95%
# Beta(1,1)    0.667  [0.390, 0.891]
# Beta(0.5,0.5) 0.682  [0.394, 0.907]
# Beta(10,10)  0.567  [0.389, 0.736]
# Beta(1,10)   0.381  [0.191, 0.592]

This contrast with the frequentist confidence interval is the point not to confuse. The confidence interval guarantees a coverage frequency whatever θ is, but gives no probability about the θ of this experiment. The credible interval gives that probability, but with a frequentist coverage that depends on θ and on the prior. The second part of the program shows how much the prior weighs when the data are few: with 7 successes in 10, the Beta(1, 1) prior brings the mean to 0.667 and the Beta(10, 10), which «remembers» twenty balanced trials, to 0.567; a Beta(1, 10) prior, which held that successes are rare, brings it to 0.381 and the interval [0.191, 0.592] ends up excluding 0.7. With thousands of observations the difference between these priors vanishes, because the likelihood dominates: this is the content of the Bernstein-von Mises theorem, which in intuitive form says that the posterior becomes Gaussian and centred on the MLE.

A worked A/B test, and the trap of peeking

The classic comparison: variant A has 84 clicks out of 1200 visits (7.0%), B has 99 out of 1100 (9.0%). With a uniform prior the posteriors are Beta(85, 1117) and Beta(100, 1002). The question «what is the probability that B is better?» is the integral of P(θ_B > t) weighted by the density of θ_A at t. The program computes it with an accurate numerical integration (quad) and compares it with a simulation of 400,000 pairs of draws. At the end it answers a different question, one that A/B-testing articles usually skip: if A and B are identical and you look at the data every 100 visits per arm, stopping as soon as P(B > A) exceeds 0.95, how often do you declare a win that is not there?

import numpy as np
from scipy import stats, integrate

kA, nA, kB, nB = 84, 1200, 99, 1100
pa, pb = stats.beta(1 + kA, 1 + nA - kA), stats.beta(1 + kB, 1 + nB - kB)
# exact: integrate P(thetaB > t) * density of thetaA
exact = integrate.quad(lambda t: pb.sf(t) * pa.pdf(t), 0, 1, epsabs=1e-12, limit=500)[0]
rng = np.random.default_rng(11)
M = 400_000
a, b = pa.rvs(M, random_state=rng), pb.rvs(M, random_state=rng)
d = b - a
print(f"P(B>A) exact {exact:.4f}  simulated {np.mean(b > a):.4f} +- {np.sqrt(exact*(1-exact)/M):.4f}")
for name, dist in (("A", pa), ("B", pb)):
    print(f"{name}: posterior mean {dist.mean():.4f}  ETI 95% [{dist.ppf(.025):.4f}, {dist.ppf(.975):.4f}]")
print(f"difference: mean {d.mean():.4f}  ETI 95% [{np.quantile(d, .025):.4f}, {np.quantile(d, .975):.4f}]")

# Stop as soon as P(B>A) > 0.95: how often do we "win" by chance? ---
def prob_b_better(kA, nA, kB, nB):          # normal approximation of the two Beta posteriors
    def mv(k, n):
        a_, b_ = 1 + k, 1 + n - k
        return a_ / (a_ + b_), a_ * b_ / ((a_ + b_) ** 2 * (a_ + b_ + 1))
    (ma, va), (mb, vb) = mv(kA, nA), mv(kB, nB)
    return stats.norm.cdf((mb - ma) / np.sqrt(va + vb))

rng = np.random.default_rng(2027)
true_p, step, looks, repetitions = 0.07, 100, 20, 4000   # A and B IDENTICAL
false_with_peeking = false_at_end = 0
for _ in range(repetitions):
    ca = np.cumsum(rng.binomial(step, true_p, looks))
    cb = np.cumsum(rng.binomial(step, true_p, looks))
    n_cum = step * np.arange(1, looks + 1)
    pr = prob_b_better(ca, n_cum, cb, n_cum)
    false_with_peeking += np.any(pr > 0.95)
    false_at_end += pr[-1] > 0.95
e1, e2 = false_with_peeking / repetitions, false_at_end / repetitions
print(f"false wins, 20 peeks: {e1:.3f} +- {np.sqrt(e1*(1-e1)/repetitions):.3f}")
print(f"false wins, a single look at 2000 visits per arm: {e2:.3f} +- {np.sqrt(e2*(1-e2)/repetitions):.3f}")

# Output:
# P(B>A) exact 0.9615  simulated 0.9607 +- 0.0003
# A: posterior mean 0.0707  ETI 95% [0.0569, 0.0859]
# B: posterior mean 0.0907  ETI 95% [0.0745, 0.1084]
# difference: mean 0.0200  ETI 95% [-0.0023, 0.0425]
# false wins, 20 peeks: 0.197 +- 0.006
# false wins, a single look at 2000 visits per arm: 0.049 +- 0.003

The probability that B beats A is 0.9615, and the simulation gives 0.9607: the gap is 0.0008, about 2.6 standard errors, so within the 3 we use as a threshold (with three other seeds tried the gap stayed between 0.0000 and 0.0007; the numerical integral is the reference value). The difference between the two rates has mean 0.0200 and a 95% credible interval of [−0.0023, 0.0425]: it includes zero. There is no contradiction with the 96%: P(B > A) = 0.96 means that zero lies in the left-hand 4%, while a two-sided 95% interval leaves out 2.5% on each side. Whoever reads «96%» and announces that the new button wins with certainty is saying more than the data support.

The last line corrects a commonplace. It is true that the posterior, computed on all the data seen up to a given moment, remains a coherent description of the uncertainty even if you decided when to stop by looking at it. But a stopping rule such as «stop when P(B > A) > 0.95» produces false wins at a frequency that depends on how many times you look: in the simulation they are 19.7% (error 0.6%) with twenty peeks, against 4.9% (error 0.3%) if you look only once at 2000 visits per arm. For a decision that costs, fix the threshold and the sample sizes in advance, or check the error frequencies of your rule with the simulator. For experiment design see Test and learn: controlled experiments and A/B tests.

When the posterior cannot be computed: Metropolis

Outside conjugate models the evidence integral is intractable, but to sample from the posterior you do not need to know it: it is enough to be able to evaluate it up to a constant. The Metropolis algorithm (1953) builds a Markov chain, that is, a sequence in which each value depends only on the previous one, whose long-run distribution is the posterior (chains are the subject of the next article). At each step it proposes a new value by adding Gaussian noise of size step, and accepts it with probability min(1, ratio of the posteriors); if it rejects it, it repeats the current value. Here the target is the Beta(9, 13) already known, so you can compare the outcome with the exact answer.

The listing also shows the two diagnostics you always need: R̂ by Gelman and Rubin (1992), which compares the variance within and between four chains started from different points and is about 1 if they have reached the same distribution, and the effective sample size (ESS), the number of independent draws equivalent to the correlated draws produced.

import numpy as np
from scipy import stats

# Random-walk Metropolis on the Beta(9, 13) posterior ---------------------
def log_post(t):                            # log of theta^8 (1-theta)^12, -inf outside (0,1)
    return np.where((t > 0) & (t < 1), 8*np.log(np.clip(t, 1e-300, 1)) + 12*np.log(np.clip(1 - t, 1e-300, 1)), -np.inf)

def metropolis(step, start, T, seed):
    r = np.random.default_rng(seed)
    x = np.empty(T); x[0] = start; acc = 0
    lp = float(log_post(x[0]))
    for i in range(1, T):
        prop = x[i-1] + step * r.standard_normal()
        lpp = float(log_post(prop))
        if np.log(r.random()) < lpp - lp:   # accept with prob. min(1, ratio)
            x[i], lp, acc = prop, lpp, acc + 1
        else:
            x[i] = x[i-1]
    return x, acc / (T - 1)

def ess(x):                                 # sum of autocorrelations while they stay positive
    x = x - x.mean(); T = len(x)
    f = np.fft.rfft(x, 2*T); ac = np.fft.irfft(f*np.conj(f))[:T]; ac /= ac[0]
    s = 0.0
    for lag in range(1, T - 1, 2):          # pairs of lags, Geyer's rule
        g = ac[lag] + ac[lag+1]
        if g < 0: break
        s += g
    return T / (1 + 2*s)

def rhat(chains):                           # classic Gelman-Rubin, chains of equal length
    m, T = chains.shape
    W = chains.var(axis=1, ddof=1).mean()
    B = T * chains.mean(axis=1).var(ddof=1)
    return np.sqrt(((T-1)/T * W + B/T) / W)

exact = stats.beta(9, 13)
T, burn = 20_000, 2_000
print("step    acceptance    mean     sd      ESS   R-hat")
for step in (0.005, 0.3, 5.0):
    chains = np.array([metropolis(step, s, T, 100 + j)[0][burn:] for j, s in enumerate((0.05, 0.3, 0.7, 0.95))])
    acc = np.mean([metropolis(step, 0.3, T, 100 + j)[1] for j in range(4)])
    print(f"{step:<7} {acc:>10.3f}  {chains.mean():.4f}  {chains.std():.4f}  {sum(ess(c) for c in chains):>6.0f}  {rhat(chains):.3f}")
print(f"exact            {exact.mean():.4f}  {exact.std():.4f}")

# Output:
# step    acceptance    mean     sd      ESS   R-hat
# 0.005        0.985  0.4129  0.1036      66  1.066
# 0.3          0.388  0.4076  0.1019   16934  1.000
# 5.0          0.027  0.4099  0.1019    1245  1.001
# exact            0.4091  0.1025

The step 0.3 works: acceptance 38.8%, mean 0.4076 and standard deviation 0.1019 against the exact 0.4091 and 0.1025, R̂ = 1.000, ESS about 16,900 out of 72,000 samples (four chains of 18,000 after discarding the first 2,000). The step 0.005 accepts almost everything (98.5%) but moves in tiny increments: ESS 66, R̂ 1.066 and a wrong mean of 0.4129, because the chains have not yet explored the distribution. The step 5.0 proposes enormous jumps that almost always land where the posterior is negligible and are rejected (acceptance 2.7%): the chain stays put and the ESS drops to about 1,250. Tuning the step is therefore a compromise between accepting often and moving far. For models with many parameters more efficient versions are used, such as the Hamiltonian Monte Carlo implemented in Stan or PyMC, not run here because the libraries are not on the test machine: the principle, and the diagnostics, are those shown.

A worked exercise and a test to do on your own

A laboratory wants to estimate the fraction of defective parts. With a Beta(2, 18) prior (mean 0.10, equivalent to 20 past parts) it observes 3 defects in 25 parts. The posterior is Beta(2 + 3, 18 + 22) = Beta(5, 40): mean 5/45 = 0.111, MAP (5 − 1)/(5 + 40 − 2) = 4/43 = 0.0930. The MLE would be 3/25 = 0.12. The posterior weighs 25/45 = 56% on the data and 44% on the prior: mean = 0.56 · 0.12 + 0.44 · 0.10 = 0.111, as expected.

What you can do now

You can turn a prior and some data into a conjugate posterior, tell MAP, median and mean apart knowing which loss each one corresponds to, and read a credible interval for what it does and does not promise. You can link the Gaussian MAP to ridge, compute P(B > A) in an A/B test without taking it for a certainty, and you know that stopping at the first threshold crossed inflates false wins. When the closed form is missing, you can write a Metropolis sampler and check its step, R̂ and ESS. The next article, on Markov chains, explains why that sampler converges; the bootstrap shows the frequentist alternative for building intervals without a prior. For the machine learning version see Bayesian statistics: prior, posterior, conjugates. The previous article is linear regression.

The same series also includes Exploratory Data Analysis (EDA) with Python.