Maximum a Posteriori Estimation

Maximum a posteriori (MAP) estimation reports the value of an unknown at which its posterior density is highest: the posterior mode. Maximum Likelihood Estimation introduces it as maximum likelihood plus a log-prior and shows how priors become penalties; Bayesian Inference places it among the posterior summaries. This article covers what is particular to the mode: the loss it is optimal for, when it agrees with the posterior mean and median and when it does not, why a sparse prior makes it return exact zeros that the prior itself never produces, and why it changes when the unknown is expressed differently. An example on a synthetic wavelet image measures both effects.

Definition

The MAP estimate of a parameter $\theta$ from data $\mathcal{D}$ is the maximizer of the posterior, $\hat\theta_{\text{MAP}} = \arg\max_\theta p(\theta \mid \mathcal{D})$ (Murphy, 2023, Eq. 3.7), which Murphy calls the most probable value of the parameter.

Intuition

A noisy wavelet coefficient of an image is observed as $y = x + n$. Most true coefficients are small, a few are large. Asked for one number, the mode answers “the single most plausible $x$”. If $y$ is small, the most plausible single value is exactly zero, because a prior peaked at zero outweighs a weak measurement. The posterior mean answers a different question, “the average of all plausible $x$”, and for this symmetric prior it is exactly zero only when $y$ is (our derivation).

Formal Definition

The loss behind the mode

In Bayesian decision theory the posterior mean minimizes the expected squared loss and the posterior median the expected absolute loss (Murphy, 2022, §5.1.5.1–5.1.5.2); for parameter estimation under 0-1 loss “the optimal action is to pick the posterior mode, i.e., the MAP estimate” (Murphy, 2023, §34.1.5.3). Bayesian Inference works through the three losses, including a practical zero–one loss for a continuous quantity.

Gaussian noise and a Laplace prior

Murphy (2022, Eq. 11.74) writes the Laplace density as $\mathrm{Laplace}(x \mid 0, b) = \frac{1}{2b} \exp(-|x|/b)$. For one observation $y = x + n$ with $n \sim \mathcal{N}(0, \sigma^2)$, the negative log posterior is, up to a constant,

$$\frac{(y - x)^2}{2\sigma^2} + \frac{|x|}{b},$$

and by our derivation its minimizer is the soft threshold

$$\hat x_{\text{MAP}} = \operatorname{sign}(y)\, \max\!\big(|y| - \tau, 0\big), \qquad \tau = \sigma^2 / b .$$

Every observation with $|y| \le \tau$ is set to exactly zero. This is the scalar case of the lasso, which Linear Regression derives as MAP estimation under a Laplace prior. The posterior mean and median have no such simple form; the example integrates them numerically.

Properties

• When mode, median and mean coincide. For a Gaussian posterior the three are the same point, which is why the Wiener filter is at once the MAP estimate and the minimum mean-square-error estimate for Gaussian image and noise models. In our reasoning this extends to any posterior that is symmetric and unimodal; skewed posteriors separate them, and Murphy (2023, §3.2.1.5) notes that the mode can be a poor summary, since it picks a single point from the whole distribution.
• No uncertainty and atypical values. A MAP estimate carries no measure of uncertainty, and the mode is often untypical of the posterior: an isolated peak far from most of the probability mass, or, for a skewed distribution such as a variance posterior, at zero while the mean is not (Murphy, 2023, §7.4.1.1–7.4.1.2).
• Sparse priors give sparse MAP estimates but not sparse posterior samples. Murphy (2023, §15.2.6) notes that under a Laplace prior posterior samples are not sparse, even though the MAP estimate is, because the prior puts infinitesimal probability on any coefficient being exactly zero. A spike-and-slab prior, which puts a point mass at zero, gives sparse posterior samples; its MAP estimate corresponds to an $\ell_0$ penalty (Murphy, 2023, §15.2.5).
• Dependence on parameterization. A change of variable multiplies a density by a Jacobian and moves its maximum, so the mode of $f(\theta)$ is in general not $f$ of the mode of $\theta$. Murphy (2023, §7.4.1.3) gives this as a drawback of MAP estimation and notes that maximum likelihood does not suffer from it, since the likelihood is a function and not a density, and neither does Bayesian inference, which integrates. Medians and other quantiles do not move under a monotone transform (our reasoning); the example shows both.
• Labellings. For discrete labellings the MAP labelling and the maximizer of the posterior marginals (MPM) answer different losses; Markov Random Field compares them and the optimizers that compute them.

Computing it

A MAP estimate is an optimization problem: minimize the negative log posterior. Where the prior is a penalty this is regularization and its solvers. With latent variables, expectation–maximization still applies with the prior added to the M step; Murphy (2022, §8.7.3.4) uses this to avoid the collapsing variances of maximum likelihood for Gaussian mixtures (Expectation–Maximization, Gaussian Mixture Model). The mode is also the centre of the Laplace approximation, a Gaussian fitted at the mode with the curvature of the negative log posterior there (Murphy, 2022, §4.6.8.2), which Uncertainty Estimation describes for networks.

Examples

A sparse wavelet prior: zeros and reparameterization

The code builds a 128 × 128 image from a hand-written inverse orthonormal Haar transform with three levels. The 16 × 16 approximation band is fixed (grey level 100 with a disk at 180, times the transform’s gain of 8). The 16,128 detail coefficients are drawn independently from $\mathrm{Laplace}(0, b = 6)$, and white Gaussian noise of $\sigma = 12$ is added to the image, unclipped. Because the transform is orthonormal, the noise in the coefficient domain is again white with the same $\sigma$. The prior is exact. Each detail coefficient is estimated four ways: maximum likelihood ($\hat x = y$), MAP (the soft threshold at $\tau = 24$), and the posterior mean and median, read from a table computed by numerical integration of the posterior over an $x$-grid of step 0.02 for every $y$ on a grid of step 0.05. The approximation band has a flat prior, so all four estimators use $y$ there. Five image seeds are run.

Wavelet Transform compares soft and hard thresholding by PSNR on a real image with an estimated noise level; here the prior is known and the comparison is between Bayes estimators. Bayesian Inference compares the Bayes risks of mean, median and mode in general, so this example keeps that comparison short and concentrates on two properties of the mode.

The script also integrates each estimator’s posterior expected loss against the marginal density of $y$ to give exact (quadrature) Bayes risks per coefficient, and the exact probability $P(|y| \le \tau)$. A final block estimates a noise level from 16 pure-noise calibration pixels under one prior, $p(\sigma) \propto 1/\sigma$, which is uniform in $\log \sigma$: once with the posterior written as a density in $\sigma$ and once as a density in $\log\sigma$, on fine grids, with closed forms as checks. Checks test the transform (orthonormality; PyWavelets, sign-agnostic per band), the soft threshold (pywt.threshold; a grid argmax of the log posterior, step $10^{-3}$) and the table’s mean and median (4,000,000 draws from the joint). It runs in about 6 seconds (NumPy 2.3.5, SciPy 1.16.3, PyWavelets reporting pywt.__version__ 1.8.0, whose installed metadata folder is named 1.9.0) and printed identical output on two runs.

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 itertools
import numpy as np
import pywt
from scipy import special

SIZE, LEVELS = 128, 3                          # 128 x 128 image, 3 Haar levels
B, SIGMA = 6.0, 12.0                           # Laplace prior scale, noise sd
TAU = SIGMA**2 / B                             # MAP soft threshold, 24
SEEDS = [0, 1, 2, 3, 4]
NAMES = ["ML", "MAP", "mean", "median"]


# ---- orthonormal 2-D Haar transform, written out
def haar_fwd(img, levels=LEVELS):
    a, details = img.astype(float), []
    for _ in range(levels):
        p, q = a[0::2, 0::2], a[0::2, 1::2]
        r, s = a[1::2, 0::2], a[1::2, 1::2]
        details.append(((p + q - r - s) / 2,  # changes between rows
                        (p - q + r - s) / 2,  # changes between columns
                        (p - q - r + s) / 2))  # diagonal
        a = (p + q + r + s) / 2
    return a, details[::-1]                    # coarsest level first


def haar_inv(a, details):
    for h, v, d in details:
        out = np.empty((2 * a.shape[0], 2 * a.shape[1]))
        out[0::2, 0::2] = (a + h + v + d) / 2
        out[0::2, 1::2] = (a + h - v - d) / 2
        out[1::2, 0::2] = (a - h + v - d) / 2
        out[1::2, 1::2] = (a - h - v + d) / 2
        a = out
    return a


def flat(details):
    return np.concatenate([b.ravel() for lev in details for b in lev])


def unflat(vec, template):
    out, i = [], 0
    for lev in template:
        bands = []
        for b in lev:
            bands.append(vec[i:i + b.size].reshape(b.shape))
            i += b.size
        out.append(tuple(bands))
    return out


def soft(y, t=TAU):
    return np.sign(y) * np.maximum(np.abs(y) - t, 0.0)


# ---- posterior mean and median by 1-D numerical integration on a y-grid
yg = np.arange(-4000, 4001) * 0.05             # -200 ... 200, step 0.05
xg = np.arange(-12500, 12501) * 0.02           # -250 ... 250, step 0.02
dx = 0.02
t_mean, t_med, t_Z = np.empty(yg.size), np.empty(yg.size), np.empty(yg.size)
risk_sq = {k: np.empty(yg.size) for k in NAMES}
risk_abs = {k: np.empty(yg.size) for k in NAMES}
logprior = -np.abs(xg) / B - np.log(2 * B)
for i0 in range(0, yg.size, 250):
    y = yg[i0:i0 + 250, None]
    logp = -(y - xg) ** 2 / (2 * SIGMA**2) - np.log(np.sqrt(2 * np.pi) * SIGMA) + logprior
    m = logp.max(axis=1, keepdims=True)
    p = np.exp(logp - m)
    Z = np.trapezoid(p, dx=dx, axis=1)
    t_Z[i0:i0 + 250] = Z * np.exp(m[:, 0])    # marginal density of y
    p /= Z[:, None]
    mu = np.trapezoid(p * xg, dx=dx, axis=1)
    cdf = np.concatenate([np.zeros((p.shape[0], 1)),
                          np.cumsum((p[:, 1:] + p[:, :-1]) / 2 * dx, axis=1)], axis=1)
    med = np.array([np.interp(0.5, c, xg) for c in cdf])
    t_mean[i0:i0 + 250], t_med[i0:i0 + 250] = mu, med
    est = {"ML": y[:, 0], "MAP": soft(y[:, 0]), "mean": mu, "median": med}
    for k in NAMES:
        a = est[k][:, None]
        risk_sq[k][i0:i0 + 250] = np.trapezoid(p * (xg - a) ** 2, dx=dx, axis=1)
        risk_abs[k][i0:i0 + 250] = np.trapezoid(p * np.abs(xg - a), dx=dx, axis=1)

post_mean = lambda y: np.interp(y, yg, t_mean)
post_median = lambda y: np.interp(y, yg, t_med)
inside = np.abs(yg) <= TAU
p_zero = np.trapezoid(t_Z[inside], yg[inside])
i_zero = np.argmin(np.abs(yg))
print(f"C5 table: integral of p(y) {np.trapezoid(t_Z, yg):.10f} (1 within 1e-6); "
      f"mean, median at y=0: {abs(t_mean[i_zero]):.3e}, {abs(t_med[i_zero]):.3e} (1e-9)")
print("exact risk per detail coefficient:")
for k in NAMES:
    print(f"  {k:7s} squared {np.trapezoid(t_Z * risk_sq[k], yg):10.6f}   "
          f"absolute {np.trapezoid(t_Z * risk_abs[k], yg):9.6f}")
print(f"exact P(|y| <= tau = {TAU:g}) = {p_zero:.6f}")

# ---- C4: posterior mean and median against Monte Carlo from the joint
rng = np.random.default_rng(1000)
xs = rng.laplace(0.0, B, 4_000_000)
ys = xs + SIGMA * rng.standard_normal(xs.size)
edges = np.arange(-60, 61, 2.0)
which = np.digitize(ys, edges) - 1
dm, below = xs - post_mean(ys), (xs < post_median(ys)).astype(float)
zm, zq, nb = [], [], 0
for j in range(edges.size - 1):
    sel = which == j
    n = sel.sum()
    if n < 1000:
        continue
    nb += 1
    zm.append(dm[sel].mean() / (dm[sel].std(ddof=1) / np.sqrt(n)))
    zq.append((below[sel].mean() - 0.5) / (0.5 / np.sqrt(n)))
print(f"C4 Monte Carlo, {nb} bins: max |z| mean {np.max(np.abs(zm)):.2f}, "
      f"median {np.max(np.abs(zq)):.2f} (bar 4)")

# ---- the image: fixed approximation band, Laplace detail coefficients
yy, xx = np.mgrid[0:16, 0:16]
approx = 8.0 * np.where((yy - 7.5) ** 2 + (xx - 7.5) ** 2 <= 25, 180.0, 100.0)
template = haar_fwd(np.zeros((SIZE, SIZE)))[1]
n_det = sum(b.size for lev in template for b in lev)

results = {}
for seed in SEEDS:
    rng = np.random.default_rng(seed)
    x = rng.laplace(0.0, B, n_det)                         # true detail coefficients
    clean = haar_inv(approx, unflat(x, template))
    noise = SIGMA * rng.standard_normal((SIZE, SIZE))
    noisy = clean + noise
    a_y, d_y = haar_fwd(noisy)
    y = flat(d_y)
    w_noise = flat(haar_fwd(noise)[1])
    z_b = (np.abs(x).mean() - B) / (B / np.sqrt(n_det))
    z_s = (w_noise.std() - SIGMA) / (SIGMA / np.sqrt(2 * n_det))
    print(f"\nseed {seed}: D1 mean|x| {np.abs(x).mean():.4f} (b 6, z {z_b:.2f}), "
          f"prior zeros {int(np.sum(x == 0.0))}; D2 coefficient-domain noise sd "
          f"{w_noise.std():.4f} (sigma 12, z {z_s:.2f})")
    if seed == 0:
        v = noisy
        a0, d0 = haar_fwd(v)
        rec = haar_inv(a0, d0)
        e_rec = np.linalg.norm(rec - v) / np.linalg.norm(v)
        e_en = abs(np.sqrt(np.sum(a0**2) + np.sum(flat(d0) ** 2)) - np.linalg.norm(v)) / np.linalg.norm(v)
        print(f"C1 orthonormality: reconstruction {e_rec:.3e}, energy {e_en:.3e} (1e-12)")
        coeffs = pywt.wavedec2(v, "haar", level=LEVELS)
        worst = np.abs(coeffs[0] - a0).max()
        for lev, mine in zip(coeffs[1:], d0):
            best = min(max(np.abs(np.abs(mine[i]) - np.abs(lev[j])).max()
                           for i, j in zip(range(3), perm))
                       for perm in itertools.permutations(range(3)))
            worst = max(worst, best)
        e_thr = np.abs(pywt.threshold(y, TAU, mode="soft") - soft(y)).max()
        print(f"C2 PyWavelets {pywt.__version__}: Haar {worst:.3e} (1e-9), "
              f"soft threshold {e_thr:.3e} (1e-12)")
        tests = np.concatenate([y[:500], [0, 23.999, -23.999, 24, -24, 24.001, -24.001, 100, -100]])
        err = 0.0
        for t in tests:
            g = np.arange(-int(round((abs(t) + 1) * 1000)), int(round((abs(t) + 1) * 1000)) + 1) * 1e-3
            lp = -(t - g) ** 2 / (2 * SIGMA**2) - np.abs(g) / B
            err = max(err, abs(g[np.argmax(lp)] - soft(t)))
        print(f"C3 MAP vs grid argmax over {tests.size} values: {err:.3e} (5e-4)")
    est = {"ML": y, "MAP": soft(y), "mean": post_mean(y), "median": post_median(y)}
    print(f"{'':8s}{'coef MSE':>10s}{'coef MAE':>10s}{'img MSE':>10s}{'img MAE':>10s}"
          f"{'zeros':>7s}{'fraction':>10s}")
    row = {}
    for k in NAMES:
        img = haar_inv(a_y, unflat(est[k], template))
        e = est[k] - x
        row[k] = ((e**2).mean(), np.abs(e).mean(),
                  ((img - clean) ** 2).mean(), np.abs(img - clean).mean(),
                  int(np.sum(est[k] == 0.0)))
        print(f"{k:8s}{row[k][0]:10.4f}{row[k][1]:10.4f}{row[k][2]:10.4f}{row[k][3]:10.4f}"
              f"{row[k][4]:7d}{row[k][4] / n_det:10.4f}")
    f0 = row["MAP"][4] / n_det
    zz = (f0 - p_zero) / np.sqrt(p_zero * (1 - p_zero) / n_det)
    best = [min(NAMES, key=lambda k: row[k][c]) for c in range(4)]
    print(f"  lowest: coef MSE {best[0]}, coef MAE {best[1]}, img MSE {best[2]}, "
          f"img MAE {best[3]}; MAP zero fraction vs exact z {zz:.2f}")
    results[seed] = row

# ---- reparameterization: one prior p(sigma) ~ 1/sigma, two parameterizations
rng = np.random.default_rng(2000)
r = SIGMA * rng.standard_normal(16)              # 4 x 4 flat calibration patch
n, S = r.size, np.sum(r**2)
sg = np.arange(100, 1_000_001) * 1e-4            # sigma grid 0.01 ... 100
lp_sig = -(n + 1) * np.log(sg) - S / (2 * sg**2)
sgrid = np.arange(int(np.log(0.01) / 1e-5), int(np.log(100) / 1e-5) + 1) * 1e-5
lp_s = -n * sgrid - S * np.exp(-2 * sgrid) / 2
def grid_median(lp, g):
    p = np.exp(lp - lp.max())
    c = np.concatenate([[0], np.cumsum((p[1:] + p[:-1]) / 2 * np.diff(g))])
    return np.interp(0.5, c / c[-1], g)
mode_sig, mode_s = sg[np.argmax(lp_sig)], np.exp(sgrid[np.argmax(lp_s)])
med_sig, med_s = grid_median(lp_sig, sg), np.exp(grid_median(lp_s, sgrid))
med_cf = np.sqrt((S / 2) / special.gammaincinv(n / 2, 0.5))
print(f"\nreparameterization: n {n}, S {S:.4f}, ML sigma {np.sqrt(S / n):.6f}")
print(f"  most probable sigma, density in sigma:     {mode_sig:.6f} (closed form {np.sqrt(S / (n + 1)):.6f})")
print(f"  most probable sigma, density in log sigma: {mode_s:.6f} (closed form {np.sqrt(S / n):.6f})")
print(f"  ratio {mode_s / mode_sig:.6f} (sqrt((n+1)/n) = {np.sqrt((n + 1) / n):.6f})")
print(f"  posterior median: in sigma {med_sig:.6f}, in log sigma {med_s:.6f}, closed form {med_cf:.6f}")
print(f"C6 errors: sigma-mode {abs(mode_sig - np.sqrt(S / (n + 1))):.3e} (5e-5), "
      f"log-mode {abs(np.log(mode_s) - np.log(np.sqrt(S / n))):.3e} (5e-6), "
      f"medians rel {abs(med_sig / med_cf - 1):.3e}, {abs(med_s / med_cf - 1):.3e} (1e-6)")

Output:

C5 table: integral of p(y) 1.0000009259 (1 within 1e-6); mean, median at y=0: 8.882e-16, 2.249e-15 (1e-9)
exact risk per detail coefficient:
  ML      squared 144.000133   absolute  9.574623
  MAP     squared  59.073694   absolute  5.601429
  mean    squared  46.903864   absolute  5.194184
  median  squared  47.547151   absolute  5.155700
exact P(|y| <= tau = 24) = 0.899610
C4 Monte Carlo, 50 bins: max |z| mean 1.94, median 2.03 (bar 4)

seed 0: D1 mean|x| 6.0214 (b 6, z 0.45), prior zeros 0; D2 coefficient-domain noise sd 11.9743 (sigma 12, z -0.38)
C1 orthonormality: reconstruction 1.934e-16, energy 0.000e+00 (1e-12)
C2 PyWavelets 1.8.0: Haar 9.095e-13 (1e-9), soft threshold 7.105e-15 (1e-12)
C3 MAP vs grid argmax over 509 values: 4.913e-04 (5e-4)
          coef MSE  coef MAE   img MSE   img MAE  zeros  fraction
ML        143.3860    9.5625  143.0451    9.5311      0    0.0000
MAP        59.6593    5.6016   60.6266    6.1479  14508    0.8996
mean       47.5407    5.2079   48.6974    5.5318      0    0.0000
median     48.2112    5.1659   49.3574    5.5621      0    0.0000
  lowest: coef MSE mean, coef MAE median, img MSE mean, img MAE mean; MAP zero fraction vs exact z -0.02

seed 1: D1 mean|x| 6.0019 (b 6, z 0.04), prior zeros 0; D2 coefficient-domain noise sd 11.8194 (sigma 12, z -2.70)
          coef MSE  coef MAE   img MSE   img MAE  zeros  fraction
ML        139.6998    9.4420  139.7787    9.4239      0    0.0000
MAP        58.8691    5.6167   60.2110    6.0895  14586    0.9044
mean       46.3629    5.1714   47.9002    5.4602      0    0.0000
median     47.0957    5.1464   48.6215    5.4910      0    0.0000
  lowest: coef MSE mean, coef MAE median, img MSE mean, img MAE mean; MAP zero fraction vs exact z 2.02

seed 2: D1 mean|x| 6.0739 (b 6, z 1.56), prior zeros 0; D2 coefficient-domain noise sd 11.9844 (sigma 12, z -0.23)
          coef MSE  coef MAE   img MSE   img MAE  zeros  fraction
ML        143.6573    9.5498  143.7694    9.5681      0    0.0000
MAP        60.3993    5.6776   61.8123    6.1889  14494    0.8987
mean       47.7060    5.2233   49.3173    5.5434      0    0.0000
median     48.4729    5.1986   50.0722    5.5725      0    0.0000
  lowest: coef MSE mean, coef MAE median, img MSE mean, img MAE mean; MAP zero fraction vs exact z -0.39

seed 3: D1 mean|x| 6.0367 (b 6, z 0.78), prior zeros 0; D2 coefficient-domain noise sd 11.9471 (sigma 12, z -0.79)
          coef MSE  coef MAE   img MSE   img MAE  zeros  fraction
ML        142.7385    9.5332  142.5577    9.4944      0    0.0000
MAP        59.9082    5.6342   61.0216    6.1339  14548    0.9020
mean       47.4014    5.2001   48.7102    5.5015      0    0.0000
median     48.1453    5.1644   49.4425    5.5277      0    0.0000
  lowest: coef MSE mean, coef MAE median, img MSE mean, img MAE mean; MAP zero fraction vs exact z 1.02

seed 4: D1 mean|x| 6.0044 (b 6, z 0.09), prior zeros 0; D2 coefficient-domain noise sd 11.9924 (sigma 12, z -0.11)
          coef MSE  coef MAE   img MSE   img MAE  zeros  fraction
ML        143.8180    9.5715  143.7777    9.6062      0    0.0000
MAP        58.8773    5.6212   60.1642    6.1196  14531    0.9010
mean       46.4661    5.1875   47.9469    5.4735      0    0.0000
median     47.1237    5.1525   48.5943    5.5024      0    0.0000
  lowest: coef MSE mean, coef MAE median, img MSE mean, img MAE mean; MAP zero fraction vs exact z 0.58

reparameterization: n 16, S 1867.8993, ML sigma 10.804800
  most probable sigma, density in sigma:     10.482200 (closed form 10.482196)
  most probable sigma, density in log sigma: 10.804795 (closed form 10.804800)
  ratio 1.030775 (sqrt((n+1)/n) = 1.030776)
  posterior median: in sigma 11.035329, in log sigma 11.035329, closed form 11.035329
C6 errors: sigma-mode 4.096e-06 (5e-5), log-mode 5.007e-07 (5e-6), medians rel 2.479e-11, 1.053e-11 (1e-6)

Every check held: the data checks D1 and D2 on all five seeds (largest $|z|$ 2.70, seed 1’s D2), orthonormality C1, the PyWavelets cross-checks C2, the Monte Carlo checks C4 (largest $|z|$ 2.03 over 50 bins) and the reparameterization checks C6. Two held with small margins. C3, the MAP against the grid argmax, printed $4.913 \times 10^{-4}$ against its half-step bound of $5 \times 10^{-4}$, a margin of $8.7 \times 10^{-6}$, as a half-step bound allows. C5, the table’s normalization, printed $1.0000009259$ against a tolerance of $10^{-6}$, a margin of $7.4 \times 10^{-8}$. The same quadrature error shows in the maximum likelihood estimator’s squared risk, which is $\sigma^2 = 144$ exactly by construction and is printed as 144.000133: an error of about $10^{-4}$, far below the differences between the four estimators’ risks.

Mean, median and mode. On every seed the posterior mean has the lowest coefficient MSE and the posterior median the lowest coefficient MAE, and image MSE ranks the four estimators in the same order as coefficient MSE (mean, median, MAP, maximum likelihood): hypotheses H1a and H1b held. The exact risks rank them the same way; that part holds by construction, since each estimator minimizes its own posterior expected loss for every $y$, and checks the implementation. In pixel MAE the posterior mean, not the median, is lowest on every seed (for example 5.5318 against 5.5621 on seed 0): H1c held. Our reasoning for H1c, written before the run, is that each pixel error is a weighted sum of coefficient errors (by the construction, one per detail band at each of the three levels and one from the approximation band), and the median’s optimality for absolute loss holds per coefficient, not per pixel.

Exact zeros. On every seed the MAP estimate sets about nine in ten of the 16,128 detail coefficients to exactly zero, more than the planned half, and each seed’s fraction is within 4 binomial standard errors of the exact probability $P(|y| \le \tau)$ (the largest $|z|$ is 2.02, on seed 1), while the maximum likelihood, mean and median estimates and the true coefficients drawn from the prior contain no exact zeros: H2 held. The prior gives probability zero to any coefficient being zero, yet its MAP estimate is mostly zeros. This is the MAP half of Murphy’s (2023, §15.2.6) remark about the Laplace prior; the example draws no posterior samples, so it does not test the other half.

Reparameterization. The same posterior over the noise level has two “most probable” values: 10.482200 as a density in $\sigma$, and 10.804795 as a density in $\log\sigma$, which is the maximum likelihood value 10.804800 to within the planned $10^{-6}$ relative, because the prior is flat in $\log\sigma$. Their ratio, 1.030775, is $\sqrt{(n+1)/n}$ for $n = 16$ (by our derivation, the Jacobian $d\sigma/d\log\sigma = \sigma$ removes one power of $\sigma$ from the density), and the posterior median, 11.035329, is the same in both parameterizations and in closed form: H3 held.

The example was planned in full, with every check bar and hypothesis, before any of its code ran. Most bars follow from the computation (half a grid step for C3 and C6, four standard errors for C4, D1, D2 and the zero fraction). A planning script, the pilot script with the hypothesis measures removed and seed 0 only, then ran once to confirm that each bar was reachable; it printed no estimator error, zero fraction, risk, mode or median, and no bar or hypothesis was changed after it. The pilot and the run shown are the same script with the same seeds, with identical outputs. The plan departed from the batch’s outline of the example in eight ways, all before the pilot: (1) it fixed the parameters (three levels, $b = 6$, $\sigma = 12$, the approximation pattern); (2) it used five seeds; (3) it gave the approximation band a flat prior; (4) it added the exact risks and zero fraction as references; (5) it added the median’s Monte Carlo check, the orthonormality, normalization and data checks, and a cross-check against PyWavelets, which is not among the batch’s planned libraries; (6) it split the MSE and MAE hypothesis into coefficient, image MSE and image MAE parts, the last predicting the mean; (7) it fixed the reparameterization block (16 pixels, prior $1/\sigma$, closed forms); (8) it put the overlap with Bayesian Inference’s example to the editor, who kept the comparison short.

Common Misconceptions

• “The MAP estimate is the Bayesian answer.” It is one summary, optimal for zero–one loss. In the example the posterior mean has the lower coefficient MSE on every seed, 47.5407 against 59.6593 on seed 0.
• “A sparse estimate means the prior believes in sparsity.” The Laplace prior gives zero probability to exact zeros, yet its MAP estimate zeroed 0.8996 of the coefficients on seed 0 of the example. Posterior samples under the same prior are not sparse (Murphy, 2023, §15.2.6).
• “The most probable value is a property of the posterior.” It is a property of the posterior and a choice of coordinates: the example’s two modes of one posterior differ by 3%. The median and maximum likelihood do not change with the coordinates.

Where It Is Used

• Restoration and denoising. Image Restoration treats restoration of linear degradations as MAP estimation with a quadratic prior, Total Variation Denoising reports Getreuer’s derivation of the ROF model as a MAP estimate, and Image Denoising covers wavelet shrinkage. Image Deblurring reads its regularized filters as MAP estimates and reports Levin et al.’s (2009) analysis of MAP estimation in blind deconvolution.
• Labelling. Graph Cuts computes exact MAP labellings of binary Markov random fields; Markov Random Field compares approximate optimizers.
• Geometry and SLAM. In a factor graph with Gaussian measurement models, MAP inference is nonlinear least squares, as in Pose Graph Optimization and Bundle Adjustment.
• Training. Weight decay and other regularization penalties read as log-priors; MAP-EM fits mixtures without collapsing variances (Expectation–Maximization).