Bayesian Inference

Bayesian inference treats every unknown quantity, a model parameter as much as a hidden label, as a random variable, describes what is known about it before the data by a prior distribution, and uses Bayes’ rule to turn prior and likelihood into a posterior distribution. The posterior is the answer: a point estimate, an interval, a prediction or a decision are all read from it afterwards, and which point estimate is “best” depends on the loss attached to errors. Probability and Statistics introduces Bayes’ rule and its vocabulary, and Maximum Likelihood Estimation introduces maximum a posteriori (MAP) estimation as maximum likelihood plus a log-prior, ending with the remark that keeping the whole posterior instead of its peak is Bayesian inference. This article starts there. It covers conjugate priors, including the gamma prior for photon counts, posterior summaries and credible intervals, Bayesian decision theory, the posterior predictive distribution, the marginal likelihood and model comparison, how posteriors are computed when no closed form exists, and the parts of vision built on the idea. An example on a photon-limited image patch computes the Bayes risk of five estimators under three losses exactly and by simulation.

Definition

Bayesian inference is the computation of the posterior distribution $p(\theta \mid \mathcal{D})$ of unknown quantities $\theta$ given observed data $\mathcal{D}$ , from a prior $p(\theta)$ and a likelihood $p(\mathcal{D} \mid \theta)$ by Bayes’ rule, followed by answering questions about $\theta$ , or about future data, with expectations under that posterior. In statistics, representing uncertainty about parameters by a probability distribution rather than a point estimate is what Murphy calls inference, and the approach that uses the posterior for this is Bayesian statistics (Murphy, 2022, §4.6).

Intuition

Point a camera at a dim, uniform patch and count photons for a tenth of the usual exposure. Sixty-four pixels collect a few dozen photons between them. The average count per pixel, the maximum likelihood estimate of the brightness, is a single number with no statement of how far off it may be. Before the exposure, though, something was known: patches in this scene have a typical brightness and a typical spread. The Bayesian answer combines the two into a distribution over the brightness. With few photons the prior pulls the answer toward the typical value and the distribution is wide; with a full exposure the counts dominate, the distribution narrows, and the prior hardly matters.

The distribution then serves every later question. Which single brightness should be reported depends on what a mistake costs: if large errors are disproportionately bad, report the posterior mean; if every error costs in proportion to its size, the median; if only exact hits count, the most probable value. An interval that holds the brightness with 90% probability, the probability that the next frame exceeds a threshold, or whether a second, darker model explains the counts better, all come from the same posterior.

Formal Definition

Posterior, marginal likelihood and predictive distribution

For data $\mathcal{D}$ and parameters $\theta$ , Bayes’ rule gives (Murphy, 2022, Eq. 4.107)

p(\theta \mid \mathcal{D}) = \frac{p(\theta)\, p(\mathcal{D} \mid \theta)}{p(\mathcal{D})}, \qquad p(\mathcal{D}) = \int p(\theta')\, p(\mathcal{D} \mid \theta')\, d\theta' .

The denominator is the marginal likelihood, or evidence: the probability of the data averaged over the prior. It does not depend on $\theta$ , so it can be ignored when only the shape of the posterior is wanted, but it is central to comparing models (Murphy, 2022, §4.6, §4.6.2.10).

Prediction averages over the parameters instead of plugging in one estimate. For a conditional model with input $\mathbf{x}$ and output $y$ , the posterior predictive distribution is (Murphy, 2022, Eq. 4.108)

p(y \mid \mathbf{x}, \mathcal{D}) = \int p(y \mid \mathbf{x}, \theta)\, p(\theta \mid \mathcal{D})\, d\theta ,

which Murphy describes as Bayes model averaging over every parameter value weighted by its posterior probability. Replacing the posterior by a spike at one estimate gives the plug-in approximation, which Murphy notes can result in overfitting: after three heads in three tosses the maximum likelihood estimate predicts that tails are impossible (Murphy, 2022, §4.6.2.9). Uncertainty Estimation develops the predictive distribution for Bayesian linear regression on image features, where its variance splits into an aleatoric and an epistemic part.

Conjugate priors

A prior is conjugate to a likelihood when the posterior belongs to the same family as the prior, so that updating only changes the family’s parameters (Murphy, 2022, §4.6.1). Murphy (2023, §3.4.5) notes that conjugate priors exist essentially only for the exponential family, where updating adds the data’s sufficient statistics to the prior’s parameters and the sample size to the prior’s strength. Four standard cases:

Beta–Bernoulli. For $N_1$ successes and $N_0$ failures (pixels that are foreground or not, detections that are correct or not) with a $\mathrm{Beta}(\alpha, \beta)$ prior, the posterior is $\mathrm{Beta}(\alpha + N_1, \beta + N_0)$ . The prior parameters act as pseudo-counts, and $\alpha + \beta$ is the prior’s equivalent sample size. The posterior mean $(\alpha + N_1)/(\alpha + \beta + N)$ is a convex combination of the prior mean and the maximum likelihood estimate, with the prior’s weight $(\alpha + \beta)/(\alpha + \beta + N)$ falling as data accumulate (Murphy, 2022, §4.6.2.4, Eq. 4.124).
Dirichlet–categorical. The multi-class version: a Dirichlet prior on histogram bin probabilities adds pseudo-counts to the observed bin counts. The posterior mode adds $\alpha_k - 1$ to bin $k$ and the posterior mean adds $\alpha_k$ , so “add-one smoothing” means a Dirichlet with $\alpha_k = 2$ when it is a MAP estimate (as Maximum Likelihood Estimation uses the term, and as Murphy, 2022, Eq. 4.120 calls the Beta(2, 2) MAP) and $\alpha_k = 1$ when it is a posterior mean or predictive probability, where Murphy (2022, Eq. 4.131) names it Laplace’s rule of succession.
Gamma–Poisson. Photon counts are Poisson (Image Noise). MacKay (2003, §23.3) lists inferring the rate of a Poisson distribution among the problems in which gamma distributions arise. By our derivation, if the $n$ pixels of a uniform patch record counts $k_i \sim \mathrm{Poisson}(t\lambda)$ during an exposure $t$ , and the flux $\lambda$ has the prior $\mathrm{Gamma}(a, b)$ with shape $a$ and rate $b$ , then multiplying $\prod_i (t\lambda)^{k_i} e^{-t\lambda}$ by $\lambda^{a-1} e^{-b\lambda}$ gives

p(\lambda \mid \mathbf{k}) = \mathrm{Gamma}\!\left(\lambda \;\middle|\; a + S,\; b + n t\right), \qquad S = \sum_{i=1}^{n} k_i .

The prior acts like $a$ photons collected over an exposure of $b$ . The posterior mean $(a + S)/(b + nt)$ moves from the prior mean $a/b$ toward the maximum likelihood estimate $S/(nt)$ as the exposure grows. SciPy’s gamma takes a scale, so this posterior is stats.gamma(a + S, scale=1/(b + n*t)) (SciPy 1.16.3 source). Averaged over the prior, the total count $S$ is negative binomial (our derivation), which the example below uses to compute risks exactly.

Gaussian–Gaussian. For Gaussian measurements of known variance $\sigma^2$ and a Gaussian prior on their mean, the posterior is Gaussian, with precision equal to the prior precision plus $N/\sigma^2$ and a mean that weights the prior mean and the data mean by their precisions (Murphy, 2022, §4.6.4.1).

Most models have no conjugate prior, and even where one exists it may be too restrictive (Murphy, 2022, §4.6.5).

Properties

Posterior summaries and the loss that picks them

A common way to summarize a posterior is a point estimate and a credible interval (Murphy, 2022, §4.6.6). Bayesian decision theory says which point estimate to use. With a loss $\ell(\theta, a)$ for reporting $a$ when the truth is $\theta$ , the posterior expected loss of $a$ is $\rho(a \mid \mathcal{D}) = \mathbb{E}_{p(\theta \mid \mathcal{D})}[\ell(\theta, a)]$ , and the Bayes estimator chooses the $a$ that minimizes it (Murphy, 2022, §5.1.1, Eqs. 5.1–5.2). For a real-valued quantity:

Squared loss $(\theta - a)^2$ is minimized by the posterior mean, the minimum mean squared error (MMSE) estimate (Murphy, 2022, §5.1.5.1).
Absolute loss $|\theta - a|$ is minimized by the posterior median (Murphy, 2022, §5.1.5.2).
Zero–one loss, which costs nothing for the right answer and 1 for any other, is minimized by the posterior mode, the MAP estimate: Murphy states this for class labels (2022, §5.1.2.1) and for parameters (2023, §34.1.5.3). For a continuous quantity, where an exact hit has probability zero, a practical version counts an estimate as right when it falls within $\varepsilon$ of the truth. By our derivation, its Bayes estimator is the centre of the interval of width $2\varepsilon$ that holds the most posterior probability, which tends to the mode as $\varepsilon \to 0$ . For a $\mathrm{Gamma}(\alpha, \beta)$ posterior with $\alpha > 1$ , setting the densities at the two ends equal gives $a^* = \varepsilon \coth\!\big(\beta\varepsilon / (\alpha - 1)\big)$ , slightly above the mode $(\alpha - 1)/\beta$ .

For a skewed posterior the three summaries differ; for a $\mathrm{Gamma}(\alpha, \beta)$ posterior the mode $(\alpha - 1)/\beta$ lies a distance $1/\beta$ below the mean $\alpha/\beta$ . The MAP estimate has a further weakness: it depends on how the parameter is expressed, because a nonlinear change of variable multiplies the density by a Jacobian and moves its maximum, while (in our reasoning) the median of a monotone transform is the transform of the median. MacKay (2003, §22, p. 306) gives this basis dependence as a reason for disliking MAP, and Fisher’s argument for the invariance of maximum likelihood, in Maximum Likelihood Estimation, is the same point seen from the other side.

Averaging the loss over the data as well as the parameter gives the Bayes risk. Murphy (2023, §34.1.4) shows that the rule which minimizes the posterior expected loss separately for each data set also minimizes the Bayes risk, because the two orders of integration in $\mathbb{E}_{p(\theta)}\mathbb{E}_{p(\mathcal{D} \mid \theta)}[\ell] = \mathbb{E}_{p(\mathcal{D})}\mathbb{E}_{p(\theta \mid \mathcal{D})}[\ell]$ agree (his Eq. 34.11). This is the sense in which a Bayes estimator is optimal: on average over data generated from the prior that the analyst assumed. The guarantee is only as good as that prior, which the example below tests by giving the analyst a wrong one.

Losses for classification work the same way. Under the zero–one loss the best label is the one with the largest posterior probability; unequal costs of false positives and false negatives move the decision threshold (Murphy, 2022, §5.1.2.2), and a reject option lets a classifier abstain when no class is probable enough (§5.1.2.3). The minimum-error rule for skin pixels in Probability and Statistics and the logit adjustment for a changed class prior in Class Imbalance are both Bayes decisions.

Credible intervals

A $100(1-\alpha)\%$ credible interval is an interval that holds $1-\alpha$ of the posterior probability. The usual choice puts $\alpha/2$ in each tail, the central or equal-tailed interval, computed from the inverse CDF of the posterior; the highest posterior density region instead keeps the points of highest density and is the narrowest such interval for a unimodal posterior (Murphy, 2022, §4.6.6).

A credible interval is not a confidence interval. A 95% confidence interval is produced by a procedure that covers the fixed true parameter in 95% of repeated data sets; it does not say that the parameter lies in the one interval computed from the data at hand with probability 0.95, which is what a credible interval says (Murphy, 2023, §3.3.5.1). The two connect through coverage. When the parameter really is drawn from the prior the analyst uses, a credible interval’s coverage over repeated draws of parameter and data equals its nominal level exactly, since averaging the posterior probability of the interval over the data gives that probability (our derivation, by the same exchange of integrals as the Bayes risk). The example checks this and shows what a wrong prior does to it. Uncertainty Estimation measures coverage of predictive intervals for regression, where the intervals cover a future label rather than a parameter.

Choosing the prior

Murphy (2023, §3.2.3) names the need to specify a prior as a challenge of the Bayesian approach, one that may be difficult in large models such as neural networks, and separates four approaches:

Conjugate priors are chosen for convenience, as above.
Noninformative priors try to “let the data speak”, but there is no unique way to define one and all of them encode some knowledge, so Murphy prefers the terms diffuse, minimally informative or default prior (2023, §3.5). Two constructions are the maximum-entropy prior, uniform for a Bernoulli rate (§3.5.1), and the Jeffreys prior, proportional to the square root of the Fisher information, which makes the posterior independent of how the parameter is expressed (§3.5.2). MacKay (2003, exercise 27.1) uses the improper prior $1/\lambda$ for a photon counter’s rate, which is uniform in $\log\lambda$ .
Hierarchical priors put a prior on the prior’s own parameters, so that related quantities, such as the brightness of many patches, inform each other (Murphy, 2022, §4.6.5.2).
Empirical Bayes, or type II maximum likelihood, estimates the prior’s parameters by maximizing the marginal likelihood and then computes the posterior with them fixed (Murphy, 2023, §3.7).

In the conjugate cases above the prior’s weight in the posterior mean falls as data accumulate (Murphy, 2022, Eq. 4.124), so the choice matters most when data are few, which in our reading is the regime of photon-limited images and ill-posed restoration.

Model comparison and Occam’s razor

Comparing models is inference one level up. With a prior over models $m$ , the posterior is $p(m \mid \mathcal{D}) \propto p(\mathcal{D} \mid m)\, p(m)$ , and under the zero–one loss the best choice is the most probable model, which for a uniform prior is the one with the largest marginal likelihood $p(\mathcal{D} \mid m) = \int p(\mathcal{D} \mid \theta, m)\, p(\theta \mid m)\, d\theta$ (Murphy, 2022, §5.2.2). The marginal likelihood penalizes complexity without an explicit penalty term. A model with more parameters must spread its prior over a larger space, so less prior mass sits on the parameters that fit the data; equivalently, a model that can predict many data sets assigns each of them less probability. Murphy (2022, §5.2.3) calls this the Bayesian Occam’s razor. His example of Bayesian model selection fits polynomials of degree 1, 2 and 3, and the most probable model is the straight line for five data points and the quadratic for thirty (§5.2.2.1).

Computing the posterior

Exact computation is usually intractable except in special cases, such as conjugate models and models whose unknowns take a few discrete values (Murphy, 2022, §4.6.8). Otherwise the posterior is approximated:

Grid approximation evaluates likelihood times prior on a grid and normalizes. It handles any shape of posterior but the number of grid points grows exponentially with the dimension (§4.6.8.1). The example uses it to check a closed form.
Laplace approximation fits a Gaussian at the MAP estimate with the inverse Hessian of the negative log posterior as covariance (§4.6.8.2). It is poor for skewed posteriors and for bounded parameters, which a change of variable such as the logit can help; MacKay (2003, §27) notes that the approximation itself depends on the parameterization. Logistic Regression uses it for Bayesian logistic regression.
Variational inference chooses the member of a tractable family closest to the posterior in KL divergence, by maximizing a lower bound on the log marginal likelihood (§4.6.8.3); Evidence Lower Bound derives the bound and its use in variational autoencoders.
Markov chain Monte Carlo draws samples from the posterior without its normalizing constant, and Hamiltonian Monte Carlo adds gradient information (§4.6.8.4). Uncertainty Estimation summarizes a study that ran full-batch Hamiltonian Monte Carlo on modern network architectures.

When the unknown is a state that changes over time, the posterior is updated recursively as each frame arrives: the Bayes filter, in its Gaussian form the Kalman filter and in its sampled form the particle filter.

Examples

Photon-limited flux: the loss that picks the estimate

The code simulates 20,000 dim, uniform $8 \times 8$ patches. Each patch’s flux $\lambda$ (expected photons per pixel per unit exposure) is drawn from the true prior $\mathrm{Gamma}(4, 0.8)$ , with mean 5 and standard deviation 2.5, and each of the 64 pixels records a Poisson count at a short exposure, $t = 0.1$ , and a long one, $t = 1$ : about 0.5 and 5 photons per pixel. An analyst estimates each patch’s flux from the posterior $\mathrm{Gamma}(a' + S, b' + 64t)$ , once with the correct prior and once with a misspecified prior $\mathrm{Gamma}(4, 0.4)$ , whose mean of 10 is twice the true mean. Five estimators are compared: the posterior mean, median and mode, the maximum likelihood estimate (MLE) $S/(64t)$ , and the window centre $a^*$ , the Bayes estimator of the zero–one-type loss with $\varepsilon = 0.1$ . The prior’s strength was chosen when the example was planned, so that its four pseudo-counts are about an eighth of the 32 photons a patch collects on average at the short exposure and about an eightieth of the 320 at the long one; the example illustrates how the prior’s weight changes with exposure under that choice rather than testing a general claim about it.

The script computes each estimator’s Bayes risk under the three losses in two ways, the two sides of Murphy’s Eq. 34.11:

Monte Carlo, the average realized loss against each patch’s true flux, with its standard error and, for each loss, the paired standard error of each estimator’s difference from the one expected to win;
exactly, by summing the closed-form posterior expected loss over every possible total count $S$ , weighted by the negative binomial probability of $S$ under the true prior. The expectation is taken under the true posterior, while the estimate comes from the analyst’s posterior, which makes the sum exact under the misspecified prior too.

It also prints the coverage of 90% equal-tailed credible intervals, by simulation and exactly. Before the risks, it checks the data (prior moments, mean counts per pixel, one flux per patch), the closed-form posterior against a grid posterior (likelihood × prior normalized by the trapezoid rule on a grid of step $10^{-4}$ ), the hand-written posterior summaries against scipy.stats and scipy.optimize, the window centre against a grid search, and the hand-written negative binomial against SciPy’s. The run takes under a second (NumPy 2.3.5, SciPy 1.16.3).

import os
for v in ("OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS"):
    os.environ[v] = "1"                       # one BLAS thread, before NumPy loads
import numpy as np
from scipy import optimize, special, stats

N, R, EPS = 64, 20_000, 0.1                   # 8 x 8 patch, patches, 0-1 window
A, B = 4.0, 0.8                               # true prior Gamma(shape A, rate B)
PRIORS = {"correct": (4.0, 0.8), "misspecified": (4.0, 0.4)}  # analyst's prior
EXPOSURES = [0.1, 1.0]
NAMES = ["mean", "median", "mode", "MLE", "window"]
LOSSES = ["squared", "absolute", "0-1"]


def estimates(S, t, a, b):
    """Point estimates of the flux from the posterior Gamma(a + S, b + N t)."""
    al, be = a + S, b + N * t
    return {"mean": al / be,
            "median": special.gammaincinv(al, 0.5) / be,
            "mode": (al - 1) / be,
            "MLE": S / (N * t),
            "window": EPS / np.tanh(be * EPS / (al - 1))}  # centre of best 2*EPS window


def realized_loss(loss, d, lam):
    e = np.abs(d - lam)
    return {"squared": e ** 2, "absolute": e, "0-1": (e >= EPS) * 1.0}[loss]


def posterior_loss(loss, d, al, be):
    """E[loss(lam, d)] for lam ~ Gamma(al, rate be), in closed form."""
    F = lambda x, shape: special.gammainc(shape, be * np.maximum(x, 0.0))
    mu = al / be
    if loss == "squared":
        return al / be ** 2 + (mu - d) ** 2
    if loss == "absolute":
        return mu - d + 2 * (d * F(d, al) - mu * F(d, al + 1))
    return 1 - (F(d + EPS, al) - F(d - EPS, al))


# ---- data: one flux per patch from the true prior, Poisson counts per pixel
rng = np.random.default_rng(0)
lam = rng.gamma(A, 1 / B, R)
counts = {t: rng.poisson(t * lam[:, None], (R, N)) for t in EXPOSURES}
S = {t: counts[t].sum(axis=1) for t in EXPOSURES}

sd_mean, sd_var = np.sqrt(A / B**2 / R), np.sqrt((A / B**2) ** 2 * (2 + 6 / A) / R)
print(f"data: lambda mean {lam.mean():.4f} (prior 5, |z| {abs(lam.mean() - 5) / sd_mean:.2f}), "
      f"variance {lam.var(ddof=1):.4f} (prior 6.25, |z| {abs(lam.var(ddof=1) - 6.25) / sd_var:.2f})")
for t in EXPOSURES:
    se = np.sqrt((t * A / B / N + t**2 * A / B**2) / R)
    expected = t * lam[:, None] * np.ones((1, N))
    print(f"data: t={t}: mean count/pixel {counts[t].mean():.5f} (expected {t * A / B:.1f}, "
          f"|z| {abs(counts[t].mean() - t * A / B) / se:.2f}); S in [{S[t].min()}, {S[t].max()}]; "
          f"one rate per patch: {np.ptp(expected, axis=1).max() == 0}")

# ---- checks on 20 posteriors: first 5 patches, both exposures, both priors
grid = np.arange(1, 300_001) * 1e-4
worst = dict.fromkeys(["dens", "gmean", "gmed", "gmode", "sp", "spmode", "win"], 0.0)
for t in EXPOSURES:
    for a, b in PRIORS.values():
        for s in S[t][:5]:
            al, be = a + s, b + N * t
            logp = (s + a - 1) * np.log(grid) - (N * t + b) * grid  # likelihood x prior
            p = np.exp(logp - logp.max())
            p /= np.trapezoid(p, grid)
            post = stats.gamma(al, scale=1 / be)
            closed = post.pdf(grid)
            est = estimates(s, t, a, b)
            cdf = np.concatenate([[0], np.cumsum((p[1:] + p[:-1]) / 2 * 1e-4)])
            res = optimize.minimize_scalar(lambda x: -post.logpdf(x), method="bounded",
                                           bounds=(est["mode"] - 1, est["mode"] + 1),
                                           options={"xatol": 1e-12})
            ag = np.arange(est["mode"] - 0.01, est["mode"] + 0.01, 1e-6)
            ag_best = ag[np.argmax(post.cdf(ag + EPS) - post.cdf(ag - EPS))]
            for key, val in [("dens", np.abs(p - closed).max() / closed.max()),
                             ("gmean", abs(np.trapezoid(grid * p, grid) - est["mean"])),
                             ("gmed", abs(np.interp(0.5, cdf, grid) - est["median"])),
                             ("gmode", abs(grid[np.argmax(p)] - est["mode"])),
                             ("sp", max(abs(post.mean() - est["mean"]),
                                        abs(post.median() - est["median"]))),
                             ("spmode", abs(res.x - est["mode"])),
                             ("win", abs(ag_best - est["window"]))]:
                worst[key] = max(worst[key], val)
print(f"C1 grid posterior vs closed form: density {worst['dens']:.3e} (tol 1e-9), "
      f"mean {worst['gmean']:.3e} (1e-9), median {worst['gmed']:.3e} (1e-7), "
      f"mode {worst['gmode']:.3e} (5e-5)")
print(f"C2 scipy.stats: mean/median {worst['sp']:.3e} (1e-12), mode {worst['spmode']:.3e} (1e-5)")
print(f"C3 window centre vs grid argmax {worst['win']:.3e} (1e-5)")

# ---- Bayes risk: Monte Carlo over the patches, and exact by summing over S
zs = []
for t in EXPOSURES:
    p_nb = B / (B + N * t)                    # true marginal of S: NB(A, p_nb)
    smax = int(stats.nbinom.isf(1e-16, A, p_nb)) + 1
    s_all = np.arange(smax + 1)
    w = np.exp(special.gammaln(A + s_all) - special.gammaln(A) - special.gammaln(s_all + 1)
               + A * np.log(p_nb) + s_all * np.log1p(-p_nb))
    print(f"\nC4 t={t}: S_max {smax}, |hand - scipy NB| {np.abs(w - stats.nbinom.pmf(s_all, A, p_nb)).max():.3e}, "
          f"1 - sum {1 - w.sum():.3e}, E[S] - NtA/B {(s_all * w).sum() - N * t * A / B:.3e}")
    al_true, be_true = A + s_all, B + N * t   # true posterior for each S
    for pname, (a, b) in PRIORS.items():
        print(f"\nt={t}, {pname} prior Gamma({a:g}, {b:g})")
        print(f"{'loss':9s}{'estimator':10s}{'MC risk':>10s}{'SE':>10s}{'exact risk':>16s}"
              f"{'z':>7s}{'paired SE':>11s}")
        mc_est, ex_est = estimates(S[t], t, a, b), estimates(s_all, t, a, b)
        for loss in LOSSES:
            mc = {k: realized_loss(loss, mc_est[k], lam) for k in NAMES}
            ex = {k: (w * posterior_loss(loss, ex_est[k], al_true, be_true)).sum() for k in NAMES}
            ref = {"squared": "mean", "absolute": "median", "0-1": "window"}[loss]
            for k in NAMES:
                se = mc[k].std(ddof=1) / np.sqrt(R)
                z = (mc[k].mean() - ex[k]) / se
                zs.append(z)
                pse = "-" if k == ref else f"{(mc[k] - mc[ref]).std(ddof=1) / np.sqrt(R):.2e}"
                print(f"{loss:9s}{k:10s}{mc[k].mean():10.6f}{se:10.2e}{ex[k]:16.7e}{z:7.2f}{pse:>11s}")
            order = sorted(NAMES, key=lambda k: ex[k])
            four = sorted(NAMES[:4], key=lambda k: ex[k])
            print(f"  {loss}: exact best {order[0]} (margin {ex[order[1]] - ex[order[0]]:.3e} to "
                  f"{order[1]}), worst {order[-1]}; best of mean/median/mode/MLE {four[0]}; "
                  f"MC best {min(NAMES, key=lambda k: mc[k].mean())}")
        lo, hi = (special.gammaincinv(a + S[t], q) / (b + N * t) for q in (0.05, 0.95))
        cov_mc = np.mean((lam >= lo) & (lam <= hi))
        lo_s, hi_s = (special.gammaincinv(a + s_all, q) / (b + N * t) for q in (0.05, 0.95))
        cov_ex = (w * (special.gammainc(al_true, be_true * hi_s)
                       - special.gammainc(al_true, be_true * lo_s))).sum()
        z = (cov_mc - cov_ex) / np.sqrt(cov_ex * (1 - cov_ex) / R)
        zs.append(z)
        print(f"  90% interval coverage: MC {cov_mc:.4f}, exact {cov_ex:.10f}, z {z:.2f}, "
              f"MC vs 0.9 z {(cov_mc - 0.9) / np.sqrt(0.09 / R):.2f}")
print(f"\nC5 MC vs exact over {len(zs)} cells: max |z| {np.max(np.abs(zs)):.2f}, "
      f"cells with |z| > 4: {int(np.sum(np.abs(zs) > 4))}")

Output:

data: lambda mean 5.0135 (prior 5, |z| 0.76), variance 6.2365 (prior 6.25, |z| 0.16)
data: t=0.1: mean count/pixel 0.50162 (expected 0.5, |z| 0.87); S in [0, 151]; one rate per patch: True
data: t=1.0: mean count/pixel 5.01144 (expected 5.0, |z| 0.64); S in [13, 1383]; one rate per patch: True
C1 grid posterior vs closed form: density 8.565e-13 (tol 1e-9), mean 3.553e-15 (1e-9), median 3.901e-10 (1e-7), mode 4.706e-05 (5e-5)
C2 scipy.stats: mean/median 1.776e-15 (1e-12), mode 2.734e-07 (1e-5)
C3 window centre vs grid argmax 4.815e-07 (1e-5)

C4 t=0.1: S_max 395, |hand - scipy NB| 5.013e-16, 1 - sum -1.110e-15, E[S] - NtA/B 0.000e+00

t=0.1, correct prior Gamma(4, 0.8)
loss     estimator    MC risk        SE      exact risk      z  paired SE
squared  mean        0.699130  8.54e-03   6.9444444e-01   0.55          -
squared  median      0.701178  8.67e-03   6.9657869e-01   0.53   5.47e-04
squared  mode        0.718161  9.03e-03   7.1373457e-01   0.49   1.64e-03
squared  MLE         0.787331  9.36e-03   7.8125000e-01   0.65   4.32e-03
squared  window      0.717917  9.02e-03   7.1348878e-01   0.49   1.63e-03
  squared: exact best mean (margin 2.134e-03 to median), worst MLE; best of mean/median/mode/MLE mean; MC best mean
absolute mean        0.644879  3.76e-03   6.4533416e-01  -0.12   3.24e-04
absolute median      0.643465  3.79e-03   6.4421535e-01  -0.20          -
absolute mode        0.647725  3.86e-03   6.4876004e-01  -0.27   6.43e-04
absolute MLE         0.684321  3.99e-03   6.8359460e-01   0.18   1.94e-03
absolute window      0.647626  3.86e-03   6.4865941e-01  -0.27   6.37e-04
  absolute: exact best median (margin 1.119e-03 to mean), worst MLE; best of mean/median/mode/MLE median; MC best median
0-1      mean        0.895700  2.16e-03   8.9623634e-01  -0.25   2.69e-03
0-1      median      0.891100  2.20e-03   8.9508187e-01  -1.81   2.21e-03
0-1      mode        0.889200  2.22e-03   8.9413756e-01  -2.22   2.78e-04
0-1      MLE         0.899950  2.12e-03   9.0073329e-01  -0.37   2.67e-03
0-1      window      0.888850  2.22e-03   8.9413718e-01  -2.38          -
  0-1: exact best window (margin 3.774e-07 to mode), worst MLE; best of mean/median/mode/MLE mode; MC best window
  90% interval coverage: MC 0.9000, exact 0.9000000000, z -0.02, MC vs 0.9 z -0.02

t=0.1, misspecified prior Gamma(4, 0.4)
loss     estimator    MC risk        SE      exact risk      z  paired SE
squared  mean        0.806586  9.10e-03   8.0017301e-01   0.70          -
squared  median      0.780018  8.93e-03   7.7377874e-01   0.70   5.87e-04
squared  mode        0.741184  8.69e-03   7.3529412e-01   0.68   1.76e-03
squared  MLE         0.787331  9.36e-03   7.8125000e-01   0.65   3.38e-03
squared  window      0.741323  8.69e-03   7.3543227e-01   0.68   1.76e-03
  squared: exact best mode (margin 1.382e-04 to window), worst mean; best of mean/median/mode/MLE mode; MC best mode
absolute mean        0.702492  3.96e-03   7.0124675e-01   0.31   3.29e-04
absolute median      0.688382  3.91e-03   6.8729443e-01   0.28          -
absolute mode        0.667168  3.85e-03   6.6630931e-01   0.22   6.66e-04
absolute MLE         0.684321  3.99e-03   6.8359460e-01   0.18   1.88e-03
absolute window      0.667270  3.85e-03   6.6640760e-01   0.22   6.60e-04
  absolute: exact best mode (margin 9.829e-05 to window), worst mean; best of mean/median/mode/MLE mode; MC best mode
0-1      mean        0.904400  2.08e-03   9.0698843e-01  -1.24   2.65e-03
0-1      median      0.901300  2.11e-03   9.0409196e-01  -1.32   2.18e-03
0-1      mode        0.900900  2.11e-03   8.9936373e-01   0.73   2.18e-04
0-1      MLE         0.899950  2.12e-03   9.0073329e-01  -0.37   2.81e-03
0-1      window      0.900750  2.11e-03   8.9940281e-01   0.64          -
  0-1: exact best mode (margin 3.907e-05 to window), worst mean; best of mean/median/mode/MLE mode; MC best MLE
  90% interval coverage: MC 0.8843, exact 0.8868220878, z -1.10, MC vs 0.9 z -7.38

C4 t=1.0: S_max 3754, |hand - scipy NB| 1.318e-15, 1 - sum -2.442e-15, E[S] - NtA/B 4.547e-13

t=1.0, correct prior Gamma(4, 0.8)
loss     estimator    MC risk        SE      exact risk      z  paired SE
squared  mean        0.075858  8.81e-04   7.7160494e-02  -1.48          -
squared  median      0.075907  8.82e-04   7.7186942e-02  -1.45   2.00e-05
squared  mode        0.076164  8.86e-04   7.7398643e-02  -1.39   6.01e-05
squared  MLE         0.076939  8.95e-04   7.8125000e-02  -1.33   1.47e-04
squared  window      0.076134  8.86e-04   7.7372384e-02  -1.40   5.75e-05
  squared: exact best mean (margin 2.645e-05 to median), worst MLE; best of mean/median/mode/MLE mean; MC best mean
absolute mean        0.212690  1.24e-03   2.1486429e-01  -1.76   3.63e-05
absolute median      0.212703  1.24e-03   2.1482227e-01  -1.71          -
absolute mode        0.212982  1.24e-03   2.1499066e-01  -1.62   7.23e-05
absolute MLE         0.214106  1.25e-03   2.1617133e-01  -1.66   2.17e-04
absolute window      0.212943  1.24e-03   2.1495879e-01  -1.62   6.63e-05
  absolute: exact best median (margin 4.202e-05 to mean), worst MLE; best of mean/median/mode/MLE median; MC best mean
0-1      mean        0.686950  3.28e-03   6.9273157e-01  -1.76   1.46e-03
0-1      median      0.686300  3.28e-03   6.9240410e-01  -1.86   1.17e-03
0-1      mode        0.685800  3.28e-03   6.9219344e-01  -1.95   3.94e-04
0-1      MLE         0.686700  3.28e-03   6.9413642e-01  -2.27   1.82e-03
0-1      window      0.685800  3.28e-03   6.9218259e-01  -1.94          -
  0-1: exact best window (margin 1.085e-05 to mode), worst MLE; best of mean/median/mode/MLE mode; MC best mode
  90% interval coverage: MC 0.9030, exact 0.9000000000, z 1.39, MC vs 0.9 z 1.39

t=1.0, misspecified prior Gamma(4, 0.4)
loss     estimator    MC risk        SE      exact risk      z  paired SE
squared  mean        0.076986  8.94e-04   7.8363103e-02  -1.54          -
squared  median      0.076714  8.91e-04   7.8068450e-02  -1.52   2.02e-05
squared  mode        0.076329  8.87e-04   7.7639752e-02  -1.48   6.06e-05
squared  MLE         0.076939  8.95e-04   7.8125000e-02  -1.33   1.13e-04
squared  window      0.076340  8.87e-04   7.7654796e-02  -1.48   5.80e-05
  squared: exact best mode (margin 1.504e-05 to window), worst mean; best of mean/median/mode/MLE mode; MC best mode
absolute mean        0.214348  1.25e-03   2.1681262e-01  -1.98   3.64e-05
absolute median      0.213907  1.24e-03   2.1632503e-01  -1.94          -
absolute mode        0.213293  1.24e-03   2.1560290e-01  -1.86   7.27e-05
absolute MLE         0.214106  1.25e-03   2.1617133e-01  -1.66   2.08e-04
absolute window      0.213321  1.24e-03   2.1563656e-01  -1.86   6.67e-05
  absolute: exact best mode (margin 3.366e-05 to window), worst mean; best of mean/median/mode/MLE mode; MC best mode
0-1      mean        0.690600  3.27e-03   6.9599245e-01  -1.65   1.43e-03
0-1      median      0.689600  3.27e-03   6.9507666e-01  -1.67   1.14e-03
0-1      mode        0.688700  3.27e-03   6.9368130e-01  -1.52   3.74e-04
0-1      MLE         0.686700  3.28e-03   6.9413642e-01  -2.27   1.78e-03
0-1      window      0.688700  3.27e-03   6.9378544e-01  -1.55          -
  0-1: exact best mode (margin 1.041e-04 to window), worst mean; best of mean/median/mode/MLE mode; MC best MLE
  90% interval coverage: MC 0.9013, exact 0.8985964892, z 1.27, MC vs 0.9 z 0.61

C5 MC vs exact over 64 cells: max |z| 2.38, cells with |z| > 4: 0

Every check held: the data match the prior and exposures within one standard error (the one-rate-per-patch line is true by the way the counts are drawn, so it cannot fail), the grid posterior matches the closed form, the summaries match SciPy, and the window centre matches its grid search (checks C1 to C4 and the data checks). C5, the agreement of all 64 Monte Carlo cells with the exact values within four standard errors, held with a largest $|z|$ of 2.38. Exact coverage under the correct prior is 0.9 to the ten printed decimals at both exposures, and the simulated coverage is within two standard errors of it: the coverage hypothesis held.

The loss picks the estimator. Under the correct prior the exact Bayes risk ranks the mean first under squared loss, the median first under absolute loss and the window centre first under the zero–one-type loss, at both exposures: the hypothesis held. Since each of these estimators minimizes its posterior expected loss for every possible count, this outcome follows from the derivation and checks the implementation rather than testing anything new. The more informative comparison is among the four estimators of the original plan, without the window centre: the mode has the lowest exact zero–one-type risk at both exposures, a hypothesis that held, in our reading because the window centre lies only slightly above the mode, by about $\beta\varepsilon^2 / (3(\alpha - 1))$ for small $\varepsilon$ by our series expansion. The exact margin between window centre and mode is $3.774 \times 10^{-7}$ at the short exposure.

Simulation cannot always see it. The Monte Carlo risks name the same winners at the short exposure: that hypothesis held there, and at the long exposure under squared loss. It failed in two cells at the long exposure. Under absolute loss the simulated risk of the mean, 0.212690, is below the median’s, 0.212703, against an exact margin of $4.202 \times 10^{-5}$ in the median’s favour and a paired standard error of $3.63 \times 10^{-5}$ for the difference. Under the zero–one-type loss the mode and the window centre both print 0.685800. Each zero–one risk is a count of misses divided by 20,000, a multiple of $5 \times 10^{-5}$ , so six decimals show that the two estimators missed the same number of patches; the code’s min() returns the first of tied names in its list, which is the mode, so the printed winner is the mode and the hypothesis that the window centre wins failed as printed, by a tie rather than a lower risk. The exact sums decide these cases, which is why the script computes them: in the absolute-loss cell the exact margin is of the same size as the paired standard error, too small for 20,000 patches to rank the two estimators reliably.

Little data, strong prior effect. Under the correct prior the MLE has the highest squared risk at the short exposure, exactly and by simulation: the hypothesis held. Its exact squared risk, 0.78125, is $\mathbb{E}[\lambda]/(nt) = 5/6.4$ , the Poisson variance $\lambda/(nt)$ of the MLE averaged over the prior, while the posterior mean’s, 0.69444, is the expected posterior variance $\mathbb{E}[a + S]/(b + nt)^2 = 36/7.2^2$ (both from the construction). At the long exposure the same formulas give 0.078125 for the MLE, $5/64$ , and 0.077160 for the posterior mean, $324/64.8^2$ : with ten times the photons the prior’s four pseudo-counts barely move the posterior.

A wrong prior. With the misspecified prior, no hypothesis was set. The output shows the posterior mode with the lowest exact risk under all three losses at both exposures and the posterior mean with the highest; in simulation the MLE has the lowest zero–one-type risk at both exposures. At the short exposure the exact coverage of the 90% intervals is 0.8868, and the simulated coverage, 0.8843, is 7.38 standard errors below 0.9; at the long exposure the exact coverage is 0.8986. The guarantees of the preceding paragraphs belong to the prior that generated the data, and with the wrong prior they no longer hold: the posterior mean is no longer the best estimate under squared loss, and coverage falls below 90% at both exposures, less so at the long one.

The planning of this example followed the site’s process for code examples. It was planned before any code of the example ran; a script then printed the accuracy of checks C1 to C4 on the example’s own data (the first five patches for C1 to C3), together with standard errors estimated on a different seed, and the bars for checks C1 (density $10^{-9}$ , mean $10^{-9}$ , median $10^{-7}$ , mode $5 \times 10^{-5}$ ), C2 ( $10^{-12}$ and $10^{-5}$ ), C3 ( $10^{-5}$ ) and C4 ( $10^{-12}$ for the probabilities and their sum, $10^{-8}$ for the mean count) were set after it, from its printed values. The mode bar is half the grid step, but it was chosen after that script had printed a mode error of $4.706 \times 10^{-5}$ on the same patches, and the run above prints the same $4.706 \times 10^{-5}$ , so that part of C1 is an illustration rather than a test. The four-standard-error bars of C5, of the coverage check and of the data checks, the $10^{-9}$ bar on exact coverage and every hypothesis were written in the plan before any code of the example ran, and were not changed afterwards. The plan departed from the batch’s outline of the example in eight ways, all decided before the pilot: exposure enters the posterior as $b + nt$ ; the prior values and the form of the misspecified prior were fixed; the exact Bayes risk was added because arithmetic showed that 20,000 patches could not resolve the smaller margins; the window centre was added as the exact zero–one-type Bayes estimator, splitting the loss hypothesis into an exact, a simulated and a four-estimator form; coverage is also computed exactly; checks of the window centre, the negative binomial, simulation against exact values and the data were added; the mode check uses scipy.optimize, since scipy.stats gives no mode; and the sources for the three Bayes estimators are cited section by section. The pilot run and the run shown are the same script with the same seed, and their outputs are identical.

Common Misconceptions

“The posterior mean is the Bayesian estimate.” Each point summary is optimal for a different loss: mean for squared, median for absolute, mode for zero–one. In the example’s exact risks under the correct prior, the mean wins under squared loss and loses under the other two.
“MAP estimation is Bayesian inference.” A MAP estimate is a single point, the mode, and it changes when the parameter is re-expressed (MacKay, 2003, p. 306). Keeping the posterior, or at least its spread, is what distinguishes Bayesian inference from regularized maximum likelihood (Regularization and Weight Decay cover the latter).
“A 90% credible interval is a 90% confidence interval.” It is a statement of posterior probability given the observed data. Its coverage over repeated draws of parameter and data is exactly 90% when the parameters are drawn from the analyst’s prior; under the misspecified prior in the example it is 0.8868 at the short exposure.
“A flat prior is no prior.” A prior that is flat in one parameterization is not flat in another: MacKay’s prior $1/\lambda$ for a photon rate is flat in $\log\lambda$ (2003, exercise 27.1). Every default prior encodes some knowledge (Murphy, 2023, §3.5).
“A softmax output is a posterior.” A classifier’s output is a distribution over labels for one input, produced by fitted weights; it is not a posterior over the weights, and its calibration has to be checked (Expected Calibration Error).
“The prior washes out, so it does not matter.” It washes out only as data accumulate. In photon-limited imaging and ill-posed restoration the data can be too few for that, in our reading, and the prior shapes the answer, as the short exposure in the example shows.

Where It Is Used

Image restoration. Geman and Geman (1984) took a Bayesian approach to restoring degraded images: a Markov random field prior on the image, expressed as a Gibbs distribution, combined with a model of blur and noise gives a posterior that is again a Markov random field, and their stochastic relaxation with gradually lowered temperature (“annealing”) generates a sequence of images that converges, in their words “in an appropriate sense”, to the MAP estimate. Graph Cuts computes such MAP labellings exactly for binary labels. Image Restoration, Total Variation Denoising and the Wiener filter, which is the posterior mean for a Gaussian image and noise model, read restoration as Bayesian estimation; Image Deblurring covers the Poisson likelihood of photon-limited deconvolution.
Low-light and scientific imaging. Photon counts in fluorescence microscopy and other dim imaging follow the Poisson model of the example, where priors matter most.
Tracking and state estimation. The Bayes filter, Kalman filter and particle filter are recursive Bayesian inference over a moving state.
Classification. Bayes decisions with class priors and costs, from the skin-colour example of Probability and Statistics to prior-shift correction in Class Imbalance and the generative classifiers of Linear Discriminant Analysis.
Uncertainty in deep learning. Bayesian neural networks, Laplace approximations, Monte Carlo dropout and their checks against exact posteriors, in Uncertainty Estimation; variational inference and latent-variable generative models in Evidence Lower Bound.
Mixtures and latent variables. Expectation–maximization for Gaussian mixture models computes posterior responsibilities of the components; MAP estimation and regularization connect priors to penalties.
Choosing hyperparameters. Empirical Bayes picks a prior’s parameters by the marginal likelihood (Murphy, 2023, §3.7); Bayesian optimization, a different use of the name, searches hyperparameters with a probabilistic model of the score (Hyperparameter Tuning).