The Bayesian Workflow

A Bayesian model describes uncertainty about unknown parameters through a posterior distribution. That uncertainty depends on the prior, the likelihood, and how well the computation approximates the posterior. A narrow interval can come from an unsuitable model; several chains can agree while exploring only part of a distribution. We will separate these problems using models small enough to check directly.

Start with the outcome the model should explain

For a conversion rate, let \(p\) be the unknown probability of success and let \(s\) be the number of successes in \(n\) trials. The binomial likelihood assumes independent trials with the same probability. A Beta prior puts a distribution on \(p\) between 0 and 1. Bayes’ rule combines them: \(\pi(p\mid s)\propto \Pr(S=s\mid p)\pi(p)\). Here \(\pi(p)\) is the prior density, \(\Pr(S=s\mid p)\) is the binomial likelihood evaluated at the observed count, and \(\pi(p\mid s)\) is the posterior density; the proportionality constant makes the posterior integrate to 1.

With a \(\operatorname{Beta}(a,b)\) prior, the posterior is \(\operatorname{Beta}(a+s,b+n-s)\). For one success in five trials, a uniform Beta(1,1) prior becomes Beta(2,5), with mean \(2/7\). This is an exact calculation: no Markov chain Monte Carlo (MCMC) is needed. A posterior mean can be a useful point estimate, while the distribution records the remaining uncertainty.

Before fitting, ask what outcomes the prior and likelihood produce together. A prior predictive check draws a rate from the prior and then successes from the binomial model. The beta-binomial distribution below computes the same mixture analytically for 20 trials. Its mean and probability of at least 18 successes translate abstract prior parameters into observable outcomes.

import numpy as np
from scipy import stats

priors = ((1, 1, "uniform"), (2, 8, "mean 0.20"), (30, 3, "mean 0.91"))
for a, b, name in priors:
    expected = 20 * a / (a + b)
    tail = stats.betabinom.sf(17, 20, a, b)
    print(f"{name:10s} expected successes {expected:.2f}  P(S >= 18) {tail:.4f}")
# uniform    expected successes 10.00  P(S >= 18) 0.1429
# mean 0.20  expected successes 4.00  P(S >= 18) 0.0001
# mean 0.91  expected successes 18.18  P(S >= 18) 0.7168

If near-universal conversion would be surprising before seeing these data, the prior centered near 0.91 needs justification. This check uses domain knowledge about plausible outcomes; it does not prove that a prior is correct. Dependence between trials or varying conversion rates would also require changing the likelihood, even if the prior looked reasonable.

How long does the prior matter?

import numpy as np
from scipy import stats

true_p = 0.30
for n in (5, 20, 100, 1000, 10000):
    g = np.random.default_rng(0)
    y = g.binomial(1, true_p, n)
    print(f"n={n:6d}  data mean {y.mean():.4f}   ", end="")
    for a, b, name in ((1, 1, "flat"), (2, 8, "mean 0.2"), (30, 3, "mean 0.91")):
        pa, pb = a + y.sum(), b + n - y.sum()
        lo, hi = stats.beta.ppf([0.025, 0.975], pa, pb)
        print(f"{name}: {pa / (pa + pb):.4f} [{lo:.3f},{hi:.3f}]  ", end="")
    print()
# n=     5  data mean 0.2000   flat: 0.2857 [0.043,0.641]  mean 0.2: 0.2000 [0.047,0.428]  mean 0.91: 0.8158 [0.680,0.920]  
# n=    20  data mean 0.4000   flat: 0.4091 [0.218,0.616]  mean 0.2: 0.3333 [0.179,0.508]  mean 0.91: 0.7170 [0.590,0.829]  
# n=   100  data mean 0.3700   flat: 0.3725 [0.282,0.468]  mean 0.2: 0.3545 [0.268,0.446]  mean 0.91: 0.5038 [0.419,0.588]  
# n=  1000  data mean 0.3220   flat: 0.3224 [0.294,0.352]  mean 0.2: 0.3208 [0.292,0.350]  mean 0.91: 0.3408 [0.312,0.370]  
# n= 10000  data mean 0.3001   flat: 0.3001 [0.291,0.309]  mean 0.2: 0.3000 [0.291,0.309]  mean 0.91: 0.3021 [0.293,0.311]

The means get closer as the sample grows, but they do not agree to three decimal places at 10,000 observations: they are 0.3001, 0.3000, and 0.3021. At 1,000, the concentrated Beta(30,3) prior still gives 0.3408 versus 0.3224 under the uniform prior. Each sample is a prefix of the same simulated Bernoulli sequence; these rows illustrate prior sensitivity in one realization.

The posterior mean is a weighted average: \[\frac{a+s}{a+b+n}=\frac{a+b}{a+b+n}\frac{a}{a+b}+\frac{n}{a+b+n}\frac{s}{n}.\] The prior contributes its mean with weight \((a+b)/(a+b+n)\). For Beta(30,3), that weight is about 3.2% at 1,000 trials and 0.33% at 10,000. This explains the remaining pull toward 0.91. It is a property of this fixed-prior binomial model, not a guarantee that data overcome every prior or resolve an unidentified parameter.

At five observations, the Beta(30,3) posterior interval [0.680, 0.920] excludes the generating rate 0.30. Beta(2,8) is less concentrated and centered closer to that rate, but it is not a uniquely correct prior for a fixed parameter of 0.30. Prior information can reduce uncertainty when relevant and distort an answer when poorly matched. Report sensitivity across substantively plausible choices, including whether they change the decision the estimate is intended to support.

Check a random-walk sampler against an exact answer

For larger models, integrating the posterior directly may be impractical. MCMC approximates posterior expectations using a sequence of dependent parameter draws, called a chain. The transitions are designed to preserve the target distribution, but a finite run may not explore it adequately. In this example, \(y_i=\beta x_i+\varepsilon_i\), with independent \(\varepsilon_i\sim N(0,1)\), no intercept, and prior \(\beta\sim N(0,25)\). A random-walk Metropolis proposal adds zero-mean Gaussian noise to the current slope. Because the proposal is symmetric, the acceptance probability is the smaller of 1 and the ratio of the new posterior density to the current density. The code evaluates that ratio on the log scale.

The exact posterior is normal with variance \(v=(\sum_i x_i^2+1/25)^{-1}\) and mean \(m=v\sum_i x_i y_i\). We use it to check both the sampled mean and spread. Proposal scale controls this particular sampler; other samplers have different settings. The first 2,000 iterations are discarded as warmup, without assuming that discarding them establishes convergence.

def rhat(chains):
    c = np.asarray(chains)
    half = c.shape[1] // 2
    c = np.concatenate((c[:, :half], c[:, -half:]), axis=0)
    W = c.var(axis=1, ddof=1).mean()
    B = half * c.mean(axis=1).var(ddof=1)
    return np.sqrt(((half - 1) / half * W + B / half) / W)

g = np.random.default_rng(1)
x = g.normal(size=200)
y = 2 * x + g.normal(size=200)
post_var = 1 / (x @ x + 1 / 25)
post_mean = post_var * (x @ y)
def lp(b):
    return -0.5 * np.sum((y - b * x) ** 2) - 0.5 * (b / 5) ** 2

def chain(seed, step, n_iter=8000):
    gg = np.random.default_rng(seed)
    b = gg.normal(0, 3)
    current = lp(b)
    out = np.empty(n_iter)
    accepted = 0
    for i in range(n_iter):
        proposal = b + gg.normal(0, step)
        proposed = lp(proposal)
        if np.log(gg.random()) < proposed - current:
            b, current = proposal, proposed
            accepted += int(i >= 2000)
        out[i] = b
    return out[2000:], accepted / (n_iter - 2000)

print(f"exact posterior: mean {post_mean:.4f}  sd {np.sqrt(post_var):.4f}")
print("step    accept  sampled mean  sampled sd  split R-hat")
for step in (0.005, 0.02, 0.1, 0.5, 2.0):
    cs, accs = zip(*[chain(seed, step) for seed in range(4)])
    cs = np.asarray(cs)
    print(f"{step:5.3f}   {np.mean(accs):.3f}      {cs.mean():.4f}       "
          f"{cs.std(ddof=1):.4f}      {rhat(cs):.4f}")
# exact posterior: mean 2.0468  sd 0.0762
# step    accept  sampled mean  sampled sd  split R-hat
# 0.005   0.974      2.0866       0.1287      1.1609
# 0.020   0.912      2.0522       0.0769      1.0193
# 0.100   0.627      2.0484       0.0755      1.0020
# 0.500   0.188      2.0453       0.0765      1.0021
# 2.000   0.050      2.0456       0.0773      1.0100

The acceptance rate here counts retained iterations only. At scale 0.005 it is 97.4%, yet the sampled mean is 2.0866 against the exact 2.0468, and the sampled spread is also too large. A tiny proposal can be accepted frequently while moving slowly; an oversized proposal can leave long runs of repeated values after rejection. Compare the resulting mean and standard deviation with the exact answer, and inspect traces for drift and sticking. Agreement in two summaries is useful evidence in this toy example, but does not establish agreement of entire distributions.

The helper implements classical split \(\widehat R\): it splits each chain in half and compares variation between halves with variation within them. This can expose disagreement across chains or over time. It is an educational calculation; modern analysis should use rank-normalized, folded split diagnostics and bulk/tail effective sample sizes from a tested implementation such as ArviZ. A value near 1 cannot certify that an unexplored region has been found.

Effective sample size (ESS) estimates Monte Carlo precision after accounting for dependence between posterior draws. It refers to independent posterior draws, not additional observed data. Precision depends on the quantity being estimated: a posterior mean and a tail quantile need different checks. Inspect Monte Carlo standard errors alongside posterior uncertainty, diagnostics, and computation time. A large ESS is useful only when the draws adequately represent the target distribution; a chain confined to one mode can still give a misleading impression of precision.

When chains explore only one mode

gg = np.random.default_rng(0)
obs = np.r_[gg.normal(-3, 1, 150), gg.normal(3, 1, 150)]

def lp2(mu):
    terms = np.logaddexp(-0.5 * (obs - mu) ** 2,
                        -0.5 * (obs + mu) ** 2) - np.log(2)
    return terms.sum() - 0.5 * (mu / 10) ** 2

def chain2(seed, start, n_iter=8000, step=0.35):
    r = np.random.default_rng(seed); mu = start; l = lp2(mu)
    out = np.empty(n_iter)
    for i in range(n_iter):
        p = mu + r.normal(0, step); lpp = lp2(p)
        if np.log(r.random()) < lpp - l: mu, l = p, lpp
        out[i] = mu
    return out[2000:]

c_same = [chain2(s, 3.0) for s in range(4)]
c_diff = [chain2(s, 3.0 if s % 2 else -3.0) for s in range(4)]
print(f"4 chains all started near +3:   split R-hat {rhat(c_same):.4f}   pooled mean {np.mean(c_same):+.4f}")
print(f"4 chains started at +3 and -3:  split R-hat {rhat(c_diff):.4f}   pooled mean {np.mean(c_diff):+.4f}")
print(f"per-chain means: {[round(float(np.mean(c)), 3) for c in c_diff]}")
# 4 chains all started near +3:   split R-hat 1.0013   pooled mean +2.9004
# 4 chains started at +3 and -3:  split R-hat 54.5727   pooled mean -0.0010
# per-chain means: [-2.902, 2.898, -2.902, 2.902]

The model is an equally weighted mixture of \(N(-\mu,1)\) and \(N(\mu,1)\), with a symmetric \(N(0,100)\) prior on \(\mu\). Both likelihood and prior are unchanged by reversing its sign. The posterior therefore has equally weighted modes near −2.90 and +2.90. The four positive-start chains agree locally and yield a split R-hat near 1 while missing the negative mode. Starting chains on both sides reveals the disagreement, but no finite collection of starting points guarantees discovery of all modes.

This example also separates a parameterization problem from a prediction problem. Reversing the sign leaves the mixture’s predictions unchanged. The exact posterior mean of the signed parameter is zero by symmetry; it is a valid mean, but substituting zero into the likelihood collapses two separated components into one. Posterior prediction averages predictions over parameter draws. It does not generally equal prediction at the posterior mean.

If only the separation matters, parameterize \(|\mu|\geq0\), with the half-normal prior induced by the original symmetric normal. That removes this duplicate labeling. It does not solve unrelated forms of multimodality. Pooling two trapped chains on each side happens to reproduce the known equal mode weights here; arbitrary numbers of chains in different modes would impose arbitrary weights.

Check what the fitted model predicts

Once computation is credible, a posterior predictive check draws \(\theta^{(s)}\) from the posterior and then a replicated dataset \(y^{\mathrm{rep},(s)}\) from the likelihood conditional on that draw. The result includes parameter uncertainty and outcome randomness. A credible interval for a rate describes uncertainty about that rate; a predictive interval for future counts also includes variation in the future trials.

Compare replicated and observed distributions, including features relevant to the problem: spread, extreme counts, zeros, subgroup patterns, or temporal dependence. Exercise 1 shows a fitted Poisson model reproducing the mean level while missing the spread. A failed check suggests which assumptions to revise. A passed check establishes compatibility with that feature of the data, not overall correctness.

Revision may involve a different observation model, missing group structure, or a prior with more defensible implications. Repeat predictive checks and assess sensitivity of the quantities you intend to report. For forecasting, also evaluate predictions on appropriately held-out observations; checking data already used for fitting is not an independent estimate of future performance. Document the model revisions and preserve an untouched evaluation set when selecting among many alternatives.

Exercises

1. A precise posterior for a poorly fitting model. Fit a common-rate Poisson model to heterogeneous counts. Compare the observed variance and maximum with replicated datasets of the same size.

import numpy as np

g = np.random.default_rng(0)
lam_i = g.gamma(2.0, 3.0, 500)          # the rate varies between units
obs = g.poisson(lam_i)                  # the MODEL will assume a single rate
print(f"observed mean {obs.mean():.3f}  variance {obs.var(ddof=1):.3f}"
      f"  (Poisson implies variance == mean)")

post_a, post_b = 1 + obs.sum(), 1 + len(obs)          # Gamma posterior for one lambda
draws = g.gamma(post_a, 1 / post_b, 2000)
print(f"posterior for lambda: mean {draws.mean():.3f}"
      f"  95% [{np.quantile(draws, 0.025):.3f},{np.quantile(draws, 0.975):.3f}]")

replicated = g.poisson(draws[:, None], size=(len(draws), len(obs)))
rep_var = replicated.var(axis=1, ddof=1)
rep_max = replicated.max(axis=1)
print(f"replicated variance: {np.quantile(rep_var, 0.025):.2f}-{np.quantile(rep_var, 0.975):.2f}"
      f"   observed {obs.var(ddof=1):.2f}   ppp {np.mean(rep_var >= obs.var(ddof=1)):.4f}")
print(f"replicated max:      {np.quantile(rep_max, 0.025):.0f}-{np.quantile(rep_max, 0.975):.0f}"
      f"   observed {obs.max()}   ppp {np.mean(rep_max >= obs.max()):.4f}")
# observed mean 5.506  variance 20.222  (Poisson implies variance == mean)
# posterior for lambda: mean 5.501  95% [5.296,5.710]
# replicated variance: 4.74-6.27   observed 20.22   ppp 0.0000
# replicated max:      12-17   observed 28   ppp 0.0000

2. What kind of 95% coverage? Draw a new rate from Beta(2,8) for each experiment, then observe binomial data. Compare simulated coverage with coverage integrated over that generating distribution.

import numpy as np
from scipy import stats

for a, b, name in ((1, 1, "uniform"), (2, 8, "matched"), (30, 3, "mean 0.91")):
    for n in (10, 100, 1000):
        g = np.random.default_rng(n)
        truth = g.beta(2, 8, 3000)
        successes = g.binomial(n, truth)
        lo = stats.beta.ppf(0.025, a + successes, b + n - successes)
        hi = stats.beta.ppf(0.975, a + successes, b + n - successes)
        coverage = np.mean((lo <= truth) & (truth <= hi))
        counts = np.arange(n + 1)
        lower = stats.beta.ppf(0.025, a + counts, b + n - counts)
        upper = stats.beta.ppf(0.975, a + counts, b + n - counts)
        mass = (stats.beta.cdf(upper, 2 + counts, 8 + n - counts)
                - stats.beta.cdf(lower, 2 + counts, 8 + n - counts))
        integrated = stats.betabinom.pmf(counts, n, 2, 8) @ mass
        print(f"{name:9s} n={n:4d} simulated {coverage:.4f} integrated {integrated:.6f}")
# uniform   n=  10 simulated 0.9533 integrated 0.955719
# uniform   n= 100 simulated 0.9440 integrated 0.950885
# uniform   n=1000 simulated 0.9537 integrated 0.950079
# matched   n=  10 simulated 0.9570 integrated 0.950000
# matched   n= 100 simulated 0.9423 integrated 0.950000
# matched   n=1000 simulated 0.9533 integrated 0.950000
# mean 0.91 n=  10 simulated 0.0000 integrated 0.000249
# mean 0.91 n= 100 simulated 0.0280 integrated 0.027251
# mean 0.91 n=1000 simulated 0.4797 integrated 0.477888

3. Forty small groups. Compare separate estimates, one shared estimate, and beta-binomial empirical Bayes estimates. Account explicitly for what is estimated from the same groups.

import numpy as np
from scipy.optimize import minimize
from scipy.special import betaln

g = np.random.default_rng(3)
theta = g.beta(9, 21, 40)
n_j = g.integers(5, 60, 40)
y_j = g.binomial(n_j, theta)
p_hat = y_j / n_j
pooled = y_j.sum() / n_j.sum()

def objective(log_ab):
    a, b = np.exp(log_ab)
    return -np.sum(betaln(a + y_j, b + n_j - y_j) - betaln(a, b))

fits = [minimize(objective, np.log(start), method="L-BFGS-B",
                 bounds=[(-8, 8), (-8, 8)]) for start in ((1, 1), (2, 8), (9, 21))]
if not all(f.success for f in fits):
    raise RuntimeError("Marginal likelihood optimization failed")
fit = min(fits, key=lambda f: f.fun)
if np.any(np.abs(fit.x) > 7.99):
    raise RuntimeError("Hyperparameter estimate reached the search boundary")
a, b = np.exp(fit.x)
w = n_j / (n_j + a + b)
shrunk = (y_j + a) / (n_j + a + b)
print(f"estimated prior: a={a:.3f} b={b:.3f} mean={a/(a+b):.3f}")
for name, estimate in (("no pooling", p_hat),
                        ("complete pooling", np.full(40, pooled)),
                        ("empirical Bayes", shrunk)):
    print(f"{name:17s} MSE vs truth {np.mean((estimate-theta)**2):.6f}")
print(f"data weights: smallest group {w[np.argmin(n_j)]:.3f}, "
      f"largest group {w[np.argmax(n_j)]:.3f}")
print(f"five trials: zero successes {a/(a+b+5):.3f}, "
      f"five successes {(a+5)/(a+b+5):.3f}")
# estimated prior: a=18.306 b=43.347 mean=0.297
# no pooling        MSE vs truth 0.006271
# complete pooling  MSE vs truth 0.007053
# empirical Bayes   MSE vs truth 0.004727
# data weights: smallest group 0.075, largest group 0.489
# five trials: zero successes 0.275, five successes 0.350

Further reading

The Stan guide to prior and posterior predictive checks develops simulation-based model checking. Vehtari and colleagues’ diagnostic examples explain rank-normalized R-hat and bulk/tail ESS; the ArviZ diagnostic implementation provides the corresponding routines.