Bayesian inference

Classical statistics gives a point estimate with a confidence interval. Bayesian inference gives a distribution over parameter values.

$\underbrace{p(\theta \mid D)}_{\text{the posterior}} \propto \underbrace{p(D \mid \theta)}_{\text{the likelihood}} \cdot \underbrace{p(\theta)}_{\text{the prior}}$

The difference shows most when the data is scarce. An A/B test with 3 conversions out of 10 visits gives the point estimate 30 %. That is not wrong, but it is not very informative either — the Bayesian posterior shows that anything between 10 % and 60 % is plausible.

The last row is what makes Bayesian intervals easier to communicate: they actually mean what people already think confidence intervals mean.

Conjugate priors give the posterior in closed form — no sampling is needed.

Beta–Binomial is the most useful in practice. The interpretation of the prior is concrete: Beta$(\alpha, \beta)$ corresponds to having seen $\alpha - 1$ successes and $\beta - 1$ failures in advance.

A Bayesian A/B test solves several practical problems at once:

1. The question becomes «how likely is it that B is better than A?» — directly interpretable.
2. You can look whenever you like without inflating the error risk, unlike repeated significance tests.
3. You get an expected loss from choosing wrongly, which is what the decision actually turns on.

Point 2 is a genuine advantage, but it is often misunderstood: you do not escape the data being noisy, and stopping as soon as the posterior looks good still gives bad decisions. What you escape is the specific multiplicity correction.

When the data becomes plentiful the prior matters hardly at all — the likelihood dominates. Bayesian methods are therefore worth the most with small datasets, hierarchical structures and when prior knowledge genuinely exists.

Hierarchical models are the strongest argument: if you have 50 schools with differing numbers of pupils you can let the schools share a common distribution. Small schools «borrow strength» from the large ones and are pulled towards the group mean — a partial pooling that neither treats everyone alike nor each one in isolation.

When there is no closed form, MCMC (NUTS in Stan or PyMC) or variational inference is used. Always check the convergence: $\hat{R} < 1.01$ and a sufficient effective sample size.

import numpy as np
from math import lgamma, exp

# Beta–Binomial: the posterior in closed form
def posterior(alpha_prior, beta_prior, successes, trials):
    return alpha_prior + successes, beta_prior + trials - successes

def beta_quantile(a, b, q, n=200_000, seed=0):
    x = np.random.default_rng(seed).beta(a, b, n)
    return float(np.quantile(x, q))

# 3 conversions out of 10 — what do we actually know?
a, b = posterior(1, 1, 3, 10)
print(f"posterior Beta({a}, {b})")
print(f"  mean {a / (a + b):.3f}")
print(f"  95 % interval [{beta_quantile(a, b, 0.025):.3f}, {beta_quantile(a, b, 0.975):.3f}]")
# posterior Beta(4, 8)
#   mean 0.333
#   95 % interval [0.109, 0.612]     ← the point estimate 0.30 hides this uncertainty

# The prior's influence shrinks with the data
for n, k in ((10, 3), (100, 30), (1000, 300)):
    for name, (ap, bp) in (("uniform", (1, 1)), ("strong 10 %", (10, 90))):
        a, b = posterior(ap, bp, k, n)
        print(f"  n={n:>4} {name:<12} posterior mean {a / (a + b):.4f}")
#  ↑ at n=1000 the two priors are practically impossible to tell apart

# A Bayesian A/B test
def ab_test(k_a, n_a, k_b, n_b, prior=(1, 1), n_sim=200_000, seed=0):
    rng = np.random.default_rng(seed)
    pa = rng.beta(prior[0] + k_a, prior[1] + n_a - k_a, n_sim)
    pb = rng.beta(prior[0] + k_b, prior[1] + n_b - k_b, n_sim)
    b_better = float((pb > pa).mean())
    return {
        "P(B > A)": round(b_better, 4),
        "expected_lift": round(float((pb - pa).mean()), 4),
        "95 %_interval_lift": [round(float(np.quantile(pb - pa, q)), 4) for q in (0.025, 0.975)],
        "expected_loss_of_choosing_B": round(float(np.maximum(pa - pb, 0).mean()), 5),
    }

print(ab_test(45, 500, 60, 500))
# {'P(B > A)': 0.9..., 'expected_lift': 0.03, ...}

# A hierarchical model: small groups borrow strength from the large ones
def partial_pooling(k, n, iterations=200):
    """Empirical Bayes: estimate a common prior from the data, pull the groups towards it."""
    k, n = np.asarray(k, float), np.asarray(n, float)
    p_global = k.sum() / n.sum()
    strength = 10.0
    for _ in range(iterations):
        a, b = strength * p_global, strength * (1 - p_global)
        post = (a + k) / (a + b + n)
        p_global = float(post.mean())
    return {"raw": [round(float(x), 3) for x in k / n],
            "pooled": [round(float(x), 3) for x in post],
            "group_mean": round(p_global, 3)}

print(partial_pooling(k=[1, 12, 48], n=[2, 40, 160]))
# {'raw': [0.5, 0.3, 0.3], 'pooled': [0.348, 0.303, 0.301], 'group_mean': 0.317}
#  ↑ the group with 2 observations is pulled hard towards the mean; the one with 160 barely at all