How old is it? Radiocarbon calibration and Bayesian chronologies

#The question, in one paragraph

Living things take in carbon from the air, including a tiny, slowly decaying radioactive kind. After death the clock starts: measure how much is left and you know how long ago the plant or animal died. In 1988 three laboratories did this with a postage-stamp-sized piece of the Shroud of Turin, the linen cloth many believe wrapped the body of Jesus. They reported that the flax was harvested between AD 1260 and 1390. How sure can we be - and could the cloth still be from the first century? The catch is that the radioactive clock does not run at a steady speed through history: the amount of radiocarbon in the air has wobbled, so a lab measurement has to be translated into a calendar date through a wiggly "exchange-rate" curve. Sometimes that translation is sharp, sometimes it smears one measurement over centuries. In the second half we see how archaeologists and geologists beat the smear: when samples come from layers that must be in order - deeper is older - the order itself sharpens every date.

#The technical plan

A radiocarbon laboratory reports a conventional radiocarbon age (in "14C years BP", before 1950) with a standard error. The true calendar date $\theta$ enters only through the calibration curve, $\mu(\theta) \pm s(\theta)$: the radiocarbon age that material of calendar age $\theta$ would have, measured on tree rings of known age (IntCal20, Reimer et al. 2020, Radiocarbon 62). Calibration is Bayes' theorem through that curve (Bronk Ramsey 2009, Radiocarbon 51; Buck, Cavanagh & Litton 1996):

$$p(\theta \mid y) \propto p(\theta)\, \mathcal N\bigl(y \mid \mu(\theta),\ \sigma^2 + s(\theta)^2\bigr).$$

Because $\mu$ wiggles, the posterior is often multimodal, and the likelihood surface is a gradient-based sampler's nightmare. Our strategy throughout: the calendar date is one- dimensional, so integrate it out exactly on a grid, and let NUTS handle only the smooth parameters around it (laboratory offsets, phase boundaries, outlier rates).

import bz2
import gzip
import json
import logging
import lzma
import struct
from pathlib import Path

import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import plotly.io as pio
import pymc as pm
import pytensor.tensor as pt
from IPython.display import HTML, display
from matplotlib import animation
from matplotlib.patches import Ellipse, Rectangle

from pymc_challenges import data

RANDOM_SEED = 1988
rng = np.random.default_rng(RANDOM_SEED)
az.style.use("arviz-variat")
pio.renderers.default = "plotly_mimetype+notebook_connected"
logging.getLogger("pymc").setLevel(logging.WARNING)

# One hue per role: radiocarbon (lab) scale, calendar scale, the sequence model; grey for data.
C14, CAL, SEQ, GREY, INK = "#7b3294", "#d6452a", "#2a78d6", "#8a8a86", "#222222"
LABS = ["Arizona", "Oxford", "Zurich"]
LAB_COL = {"Arizona": "#e08214", "Oxford": "#2a78d6", "Zurich": "#1b9e77"}


def ad_label(year):
    """Astronomical year (0 = 1 BC) -> 'AD 1260' / '250 BC'."""
    year = int(round(year))
    return f"AD {year}" if year > 0 else f"{1 - year} BC"


print(f"PyMC {pm.__version__}, ArviZ {az.__version__}")

PyMC 6.3.2, ArviZ 1.3.0

#1 · The calibration curve

IntCal20 gives, for every calendar year back to 55,000 years ago, the radiocarbon age that a sample of that age should show, with a 1-sd band. Calendar ages are in cal BP (calendar years before 1950). We convert to AD with $\text{AD} = 1950 - \text{BP}$ and use astronomical year numbering, where year 0 is 1 BC.

data.describe("intcal20")
curve = data.load("intcal20").sort_values("cal_bp").reset_index(drop=True)
print(curve.shape, "rows; calendar step (years) by range:",
      {f"{lo}-{hi} cal BP": sorted(np.unique(np.diff(curve.cal_bp[(curve.cal_bp >= lo) & (curve.cal_bp <= hi)])).tolist())
       for lo, hi in [(0, 14000), (14000, 15000), (15000, 25000), (25000, 55000)]})
curve.head(3)

intcal20: Calibration curve on a calendar grid (5-yr steps in the Holocene, coarser further back): calendar age (cal BP, years before 1950), radiocarbon age (14C yr BP) and its 1-sd uncertainty, and Delta-14C (per mil) with its 1-sd uncertainty. Rows run from 55,000 cal BP down to 0.
  source: Reimer et al. (2020), The IntCal20 Northern Hemisphere radiocarbon age calibration curve (0-55 cal kBP), Radiocarbon 62:725-757, doi:10.1017/RDC.2020.41; distributed by intcal.org.
  url:    https://intcal.org/curves/intcal20.14c
(9501, 5) rows; calendar step (years) by range: {'0-14000 cal BP': [1, 5], '14000-15000 cal BP': [5], '15000-25000 cal BP': [10], '25000-55000 cal BP': [20]}

fig, axes = plt.subplots(1, 2, figsize=(13, 4.4), gridspec_kw={"width_ratios": [1, 1.3]})
ax = axes[0]
ax.plot(curve.cal_bp / 1000, curve.c14_bp / 1000, color=C14, lw=1)
ax.plot([0, 50], [0, 50], color=GREY, ls="--", lw=1, label="if radiocarbon years were calendar years")
ax.set(xlabel="calendar age (thousand years cal BP)", ylabel="radiocarbon age (thousand 14C yr BP)",
       title="IntCal20: the whole curve")
ax.legend(fontsize=9)
ax = axes[1]
z = curve[(curve.cal_bp >= 0) & (curve.cal_bp <= 3200)]
ad = 1950 - z.cal_bp
ax.fill_between(ad, z.c14_bp - 2 * z.c14_sd, z.c14_bp + 2 * z.c14_sd, color=C14, alpha=0.25, lw=0,
                label="curve ± 2 sd")
ax.plot(ad, z.c14_bp, color=C14, lw=1)
ax.axvspan(-800, -400, color=GREY, alpha=0.15)
ax.text(-790, 2150, "Hallstatt\nplateau", fontsize=9, color=INK)
ax.axvspan(1260, 1390, color=CAL, alpha=0.12)
ax.text(1000, 350, "Damon et al.'s\nshroud range", fontsize=9, color=CAL, ha="right")
ax.set(xlabel="calendar year (negative = BC)", ylabel="radiocarbon age (14C yr BP)",
       title="The last 3,200 years: a wiggly exchange rate", xlim=(-1250, 1950))
ax.legend(loc="upper right", fontsize=9);

Two features matter for everything that follows. On the left, radiocarbon years are not calendar years: 50,000 radiocarbon years correspond to 55,000 calendar years, and the relation bends. On the right, the curve wiggles, and in places it is flat: between about 750 and 400 BC (the Hallstatt plateau, named after the Iron Age culture it plagues) the radiocarbon age hardly changes for three centuries, so a measurement there cannot distinguish 750 BC from 400 BC however precise the laboratory is.

In plain words: A lab does not measure a date; it measures how much radiocarbon is left. Turning that into a year needs a conversion table built from tree rings, and the table is bumpy - some centuries all "look the same" to radiocarbon.

#2 · Calibrating one date on a grid

For a single date the posterior lives on one axis, so we simply evaluate it on every calendar year and normalise: exact, up to the 1-year grid. The curve is interpolated linearly between its 5-year nodes. We use a flat prior over 1000 BC to AD 1900.

YEAR = np.arange(-1000, 1901)                 # astronomical years: 0 = 1 BC
MU_Y = np.interp(1950 - YEAR, curve.cal_bp, curve.c14_bp)
S_Y = np.interp(1950 - YEAR, curve.cal_bp, curve.c14_sd)


def calibrate(age, sd, year=YEAR, mu=MU_Y, s=S_Y):
    """Posterior probability of each calendar year (flat prior) for one radiocarbon age."""
    v = sd**2 + s**2
    lp = -0.5 * (age - mu) ** 2 / v - 0.5 * np.log(v)
    p = np.exp(lp - lp.max())
    return p / p.sum()


def hpd_ranges(p, year=YEAR, prob=0.95):
    """Highest-posterior-density region of a grid posterior, as a list of (from, to) years."""
    order = np.argsort(p)[::-1]
    keep = np.zeros(len(p), bool)
    keep[order[: np.searchsorted(np.cumsum(p[order]), prob) + 1]] = True
    edges = np.flatnonzero(np.diff(np.r_[0, keep.astype(int), 0]))
    return [(int(year[a]), int(year[b - 1])) for a, b in zip(edges[::2], edges[1::2])]


def fmt_ranges(ranges):
    return ", ".join(f"{ad_label(a)}-{ad_label(b)[3:] if a > 0 else ad_label(b)}" for a, b in ranges)


examples = {"Shroud, pooled (Damon et al.): 691 ± 31": (691, 31),
            "A date on the Hallstatt plateau: 2450 ± 25": (2450, 25),
            "Cope of St Louis d'Anjou (control): 724 ± 20": (724, 20)}
for name, (a, e) in examples.items():
    p = calibrate(a, e)
    print(f"{name:<46} 95% HPD: {fmt_ranges(hpd_ranges(p))}")

Shroud, pooled (Damon et al.): 691 ± 31        95% HPD: AD 1273-1317, AD 1360-1388
A date on the Hallstatt plateau: 2450 ± 25     95% HPD: 750 BC-684 BC, 667 BC-634 BC, 622 BC-613 BC, 591 BC-414 BC
Cope of St Louis d'Anjou (control): 724 ± 20   95% HPD: AD 1268-1297

The shroud's pooled date, 691 ± 31 radiocarbon years, calibrates to two separate ranges - the same double range Damon et al. reported with the "intercept" method in 1989 (AD 1262-1312 and 1353-1384 on the calibration curve of the day; IntCal20 moves the ends by a few years). The plateau date, with a smaller error, spreads over more than three centuries. The cope's date lands on a steep part of the curve and calibrates to a range narrower than its own error. The figure shows how: the lab's bell curve on the radiocarbon axis "casts a shadow" through the calibration curve onto the calendar axis.

def shadow_panel(ax_main, ax_left, ax_bottom, age, sd, xlim, color=CAL, label=None):
    m = (YEAR >= xlim[0]) & (YEAR <= xlim[1])
    ax_main.fill_between(YEAR[m], MU_Y[m] - 2 * S_Y[m], MU_Y[m] + 2 * S_Y[m], color=C14, alpha=0.3, lw=0)
    ax_main.plot(YEAR[m], MU_Y[m], color=C14, lw=1.2)
    ax_main.axhspan(age - 2 * sd, age + 2 * sd, color=color, alpha=0.12, lw=0)
    yy = np.linspace(age - 4 * sd, age + 4 * sd, 200)
    ax_left.fill_betweenx(yy, 0, np.exp(-0.5 * ((yy - age) / sd) ** 2), color=C14, alpha=0.6, lw=0)
    p = calibrate(age, sd)
    ax_bottom.fill_between(YEAR[m], 0, p[m], color=color, alpha=0.7, lw=0)
    for a, b in hpd_ranges(p):
        ax_bottom.plot([a, b], [-p.max() * 0.08] * 2, color=INK, lw=3, solid_capstyle="butt")
    ax_bottom.set(xlim=xlim, ylim=(-p.max() * 0.15, p.max() * 1.1), yticks=[], xlabel="calendar year AD")
    ax_left.set(xticks=[], xlim=(1.1, 0), ylabel="radiocarbon age (14C yr BP)")
    ax_main.set(xlim=xlim, ylim=ax_left.get_ylim())
    ax_main.tick_params(labelbottom=False, labelleft=False)
    if label:
        ax_main.set_title(label, fontsize=11)


fig = plt.figure(figsize=(13, 5.2), layout="none")
for col, (name, (a, e), xlim) in enumerate([
        ("Shroud, pooled: 691 ± 31", (691, 31), (1150, 1450)),
        ("Hallstatt plateau: 2450 ± 25", (2450, 25), (-850, -300))]):
    x0 = 0.06 + col * 0.49
    ax_left = fig.add_axes([x0, 0.36, 0.06, 0.5])
    ax_main = fig.add_axes([x0 + 0.065, 0.36, 0.36, 0.5])
    ax_bot = fig.add_axes([x0 + 0.065, 0.1, 0.36, 0.24])
    shadow_panel(ax_main, ax_left, ax_bot, a, e, xlim, label=name)
    ax_left.set_ylim(a - 4.2 * e, a + 4.2 * e)
    ax_main.set_ylim(a - 4.2 * e, a + 4.2 * e)
    if col == 1:
        ax_bot.set_xlabel("calendar year (negative = BC)")
fig.suptitle("Calibration: the lab's bell curve (left) seen through the calibration curve (purple band) "
             "gives the calendar distribution (bottom; bars = 95% HPD)", fontsize=11);

In plain words: The shroud measurement points to two stretches of history, one around 1260-1310 and another around 1350-1390, because the conversion curve passes through the measured value twice. A measurement in the Iron Age can be hopeless: the curve is flat, so anywhere in 750-400 BC fits equally well.

#3 · Why not just use NUTS?

The obvious PyMC model makes the calendar date a continuous parameter and interpolates the curve with pytensor. It compiles, it samples in two seconds, it reports zero divergences - and it is wrong.

MU_T, S_T = pt.as_tensor(MU_Y), pt.as_tensor(S_Y)


def curve_at(theta):
    """Linear interpolation of the curve (mean and sd) at a continuous calendar year."""
    x = theta - YEAR[0]
    i0 = pt.clip(pt.cast(pt.floor(x), "int64"), 0, len(YEAR) - 2)
    f = x - i0
    return MU_T[i0] * (1 - f) + MU_T[i0 + 1] * f, S_T[i0] * (1 - f) + S_T[i0 + 1] * f


with pm.Model() as single_model:
    year = pm.Uniform("year", -1000.0, 1900.0)
    mu_t, s_t = curve_at(year)
    pm.Normal("y", mu_t, pt.sqrt(31.0**2 + s_t**2), observed=691.0)
    single_idata = pm.sample(random_seed=RANDOM_SEED, progressbar=False)

draws_single = single_idata.posterior["year"].to_numpy()
p_shroud_pooled = calibrate(691, 31)
print(f"divergences: {int(single_idata.sample_stats['diverging'].sum())}, "
      f"r_hat: {float(az.rhat(single_idata)['year']):.2f}, ESS: {float(az.ess(single_idata)['year']):.0f}")
for c in range(draws_single.shape[0]):
    print(f"  chain {c}: median {ad_label(np.median(draws_single[c]))}, "
          f"90% range {ad_label(np.quantile(draws_single[c], 0.05))} to {ad_label(np.quantile(draws_single[c], 0.95))}")
print(f"exact grid posterior: median {ad_label(YEAR[np.searchsorted(np.cumsum(p_shroud_pooled), 0.5)])}")

divergences: 0, r_hat: 2.85, ESS: 5
  chain 0: median AD 90, 90% range AD 88 to AD 91
  chain 1: median AD 1290, 90% range AD 1277 to AD 1305
  chain 2: median AD 494, 90% range AD 489 to AD 495
  chain 3: median 21 BC, 90% range 22 BC to 19 BC
exact grid posterior: median AD 1297

Each chain sits in a different place, most of them centuries away from the real answer, and none moves. r_hat catches it; divergences do not. What went wrong:

• Local modes everywhere. The likelihood through a wiggly curve has a local maximum at every wiggle that comes closest to the measurement, even where it is 10-20 standard errors away. A chain initialised there climbs to the nearest wiggle and stays: the valleys between wiggles are hundreds of log units deep.
• Gradients from a piecewise-linear curve are piecewise constant and jump at every node, so the leapfrog integrator's energy error is noisy even inside the right mode.

Smoothing the curve or tempering would help a little. For radiocarbon the right move is simpler: the calendar date is one-dimensional, so never sample it - integrate it. On a 1-year grid a sum over 2,900 years is exact to plotting precision and costs nothing. Whenever a model has more structure (laboratory offsets, stratigraphic order), we will sum each date out on the grid inside the model, and let NUTS see only parameters on which the likelihood is smooth. The calendar dates are then recovered afterwards, exactly, from their grid conditional distributions.

fig, ax = plt.subplots(figsize=(11, 3.4))
ax.fill_between(YEAR, 0, p_shroud_pooled, color=CAL, alpha=0.35, lw=0, label="exact grid posterior")
for c in range(draws_single.shape[0]):
    h, edges = np.histogram(draws_single[c], bins=np.arange(-1000, 1901, 5), density=True)
    ax.step(edges[:-1], h / h.max() * p_shroud_pooled.max() * 0.8, color=plt.cm.tab10(c), lw=1,
            where="post", label=f"NUTS chain {c} (rescaled)")
ax.set(xlabel="calendar year (negative = BC)", yticks=[], title="One radiocarbon date: NUTS chains vs the exact grid")
ax.legend(fontsize=8, ncol=5, loc="upper left");

#4 · The Shroud of Turin data

In April 1988 a strip was cut from the shroud and split between three accelerator laboratories, together with three textiles of known age as controls (the laboratories knew their approximate ages): linen from an Islamic-Christian-era tomb at Qasr Ibrîm in Nubia (11th-12th century AD), linen associated with an early 2nd-century-AD mummy from Thebes (radiocarbon-dated earlier by the British Museum to the 1st century BC - 1st century AD), and threads from the cope of St Louis d'Anjou (c. AD 1290-1310). Each laboratory measured several subsamples with different cleaning treatments. We use every individual measurement, transcribed from Table 1 of Damon et al. (1989, Nature 337).

data.describe("shroud_damon1989")
d = data.load("shroud_damon1989")
TEXTILES = ["shroud", "Nubia", "mummy", "cope"]
d["textile"] = [TEXTILES[s - 1] for s in d["sample"]]
lab_i = pd.Index(LABS).get_indexer(d["lab"])
smp_i = d["sample"].to_numpy() - 1
assert (lab_i >= 0).all()
y_obs, sd_obs = d.age_bp.to_numpy(float), d.sd.to_numpy(float)
w_obs = 1 / sd_obs**2


def lab_means(y, sd, smp, lab):
    """Inverse-variance weighted mean and its error for every (sample, lab) cell."""
    w = 1 / sd**2
    num, den = np.zeros((4, 3)), np.zeros((4, 3))
    np.add.at(num, (smp, lab), w * y)
    np.add.at(den, (smp, lab), w)
    return num / den, den ** -0.5


def ward_wilson_T(m, e):
    """Chi-square statistic of the three lab means around their weighted mean (2 d.f.)."""
    w = 1 / e**2
    pooled = (w * m).sum(-1, keepdims=True) / w.sum(-1, keepdims=True)
    return (w * (m - pooled) ** 2).sum(-1)


M_obs, E_obs = lab_means(y_obs, sd_obs, smp_i, lab_i)
T_obs = ward_wilson_T(M_obs, E_obs)
print(pd.crosstab(d.textile, d.lab).loc[TEXTILES])
pd.DataFrame({f"{lab} mean": M_obs[:, j].round() for j, lab in enumerate(LABS)} |
             {f"{lab} error": E_obs[:, j].round() for j, lab in enumerate(LABS)} |
             {"chi-square (2 df)": T_obs.round(1)}, index=TEXTILES)

shroud_damon1989: All 49 individual AMS measurements (14C yr BP, 1 sd) by the three laboratories (Arizona, Oxford, Zurich) on sample 1 (the Shroud of Turin) and three known-age controls: 2 = linen from a tomb at Qasr Ibrim, Nubia (11th-12th century AD), 3 = linen associated with an early 2nd-century AD mummy of Cleopatra from Thebes, 4 = threads from the cope of St Louis d'Anjou, Saint-Maximin (c. AD 1290-1310). meas_id is the paper's code (lab, sample, run, pretreatment, replicate).
  source: Damon et al. (1989), Radiocarbon dating of the Shroud of Turin, Nature 337:611-615, doi:10.1038/337611a0. HAND-TRANSCRIBED from Table 1 ('Basic data, individual measurements') of the paper; the URL is the full text of the paper (HTML) at shroud.com, not a CSV, so keep the cached file in data/. The 12 shroud rows were cross-checked against the `shroud` data of the CRAN package rice (Blaauw), and the per-lab weighted means reproduce Table 2.
  url:    https://www.shroud.com/nature.htm
lab      Arizona  Oxford  Zurich
textile                         
shroud         4       3       5
Nubia          5       3       5
mummy          5       3       3
cope           5       3       5

Our weighted lab means reproduce Damon et al.'s Table 2 (646, 750, 676 for the shroud; the Oxford values differ by a year or two; the paper's footnotes mention an adjusted Oxford replicate and Oxford's rounding of small errors up to 40). One detail matters: for its shroud mean Arizona quoted an error of ±31, not the ±17 the four individual errors imply, and with the paper's errors the shroud's chi-square is 6.4 (p = 0.04), the number behind thirty years of argument about whether the three labs really agree. With the errors implied by the individual measurements, the disagreement looks stronger still (the value above). The controls agree well.

fig, ax = plt.subplots(figsize=(12, 4.2))
xpos = {t: i for i, t in enumerate(TEXTILES)}
for j, lab in enumerate(LABS):
    sub = d[d.lab == lab]
    xs = sub.textile.map(xpos).to_numpy() + (j - 1) * 0.22 + rng.uniform(-0.05, 0.05, len(sub))
    ax.errorbar(xs, sub.age_bp - sub.textile.map(dict(zip(TEXTILES, M_obs.mean(1)))), yerr=sub.sd, fmt="o",
                ms=4, color=LAB_COL[lab], alpha=0.8, lw=1, label=lab)
    ax.errorbar(np.arange(4) + (j - 1) * 0.22 + 0.08, M_obs[:, j] - M_obs.mean(1), yerr=E_obs[:, j],
                fmt="s", ms=8, color=LAB_COL[lab], mec=INK, lw=2)
ax.axhline(0, color=GREY, lw=1)
ax.set_xticks(range(4), [f"{t}\nχ² = {T:.1f}" for t, T in zip(TEXTILES, T_obs)])
ax.set(ylabel="radiocarbon age minus the textile's\naverage of lab means (14C yr)",
       title="Individual measurements (dots) and lab means (squares), ± 1 sd")
ax.legend(ncol=3, fontsize=9);

On the shroud, Arizona's four measurements sit low and Oxford's three high; on the controls the labs overlap. Is that a shroud-specific problem (a heterogeneous or contaminated sample), a laboratory bias that also shows on the controls, or chance? A single chi-square cannot tell these apart; a hierarchical model can.

#5 · A hierarchical model of three laboratories and four textiles

For measurement $k$ of textile $i$ by laboratory $j$:

$$y_k = c_i + \delta_j + \eta_{ij} + \varepsilon_k,\qquad \varepsilon_k \sim \mathcal N(0, \sigma_k^2),$$

• $c_i \sim \mathcal N(\mu(\theta_i), s(\theta_i)^2)$: the textile's true radiocarbon age, the calibration curve at its calendar date $\theta_i$ (shared by every measurement of it);
• $\delta_j \sim \mathcal N(0, \omega^2)$: a persistent laboratory offset;
• $\eta_{ij} \sim \mathcal N(0, \tau^2)$: lab-by-sample heterogeneity (the subsample a lab received, its cleaning, a bad run) - exactly what the chi-square test is about. We give the shroud its own $\tau_{\text{shroud}}$ and the controls a shared $\tau_{\text{ctrl}}$, so the data can say whether the shroud is special;
• calendar priors: the shroud is flat over 1000 BC - AD 1900; the controls get their historical dates as uniform windows - Nubia AD 1001-1200, cope AD 1290-1310, and the mummy linen only loosely, 300 BC - AD 200 (the linen can be older than the burial).

The controls are what make laboratory offsets estimable: with a known calendar age, the curve says what each lab should have measured.

Integrating out $c_i$ and $\theta_i$. Given the offsets, the residuals $r_k = y_k - \delta_j - \eta_{ij}$ of textile $i$ share one $c_i$. Integrating a Gaussian $c_i$ out of a product of Gaussians gives, with weights $w_k = 1/\sigma_k^2$, $W_i = \sum_k w_k$ and weighted mean $\bar r_i$,

$$\int \prod_k \mathcal N(r_k \mid c, \sigma_k^2)\, \mathcal N(c \mid \mu, s^2)\, dc = \Bigl[\prod_k \mathcal N(r_k \mid \bar r_i, \sigma_k^2)\Bigr] \sqrt{2\pi / W_i}\; \mathcal N\bigl(\bar r_i \mid \mu,\ 1/W_i + s^2\bigr).$$

Only the last factor depends on the calendar date, so the date is summed out with one logsumexp over the 2,901-year grid per textile. We check the identity by brute-force numerical integration before trusting it.

M_ONE = np.zeros((4, len(d)))
M_ONE[smp_i, np.arange(len(d))] = 1          # textile membership
W_I = M_ONE @ w_obs
WINDOWS = {"shroud": (-1000, 1900), "Nubia": (1001, 1200), "mummy": (-300, 200), "cope": (1290, 1310)}
LOG_PRIOR = np.full((4, len(YEAR)), -np.inf)
for i, t in enumerate(TEXTILES):
    inside = (YEAR >= WINDOWS[t][0]) & (YEAR <= WINDOWS[t][1])
    LOG_PRIOR[i, inside] = -np.log(inside.sum())


def textile_loglik_numpy(r, i, year_idx):
    """Closed form of log p(r_k of textile i | calendar year), curve value integrated out."""
    k = smp_i == i
    rbar = (w_obs[k] * r[k]).sum() / W_I[i]
    v = 1 / W_I[i] + S_Y[year_idx] ** 2
    within = (-0.5 * w_obs[k] * (r[k] - rbar) ** 2 - np.log(sd_obs[k]) - 0.5 * np.log(2 * np.pi)).sum()
    return within + 0.5 * np.log(2 * np.pi / W_I[i]) - 0.5 * (rbar - MU_Y[year_idx]) ** 2 / v - 0.5 * np.log(2 * np.pi * v)


# brute force: integrate over the curve value c on a fine grid at a few calendar years
r_test = y_obs - rng.normal(0, 20, len(y_obs))
c_grid = np.linspace(0, 3000, 300_001)
for i, yr in [(0, 1300), (3, 1295), (2, 50)]:
    j = int(np.flatnonzero(YEAR == yr)[0])
    k = smp_i == i
    log_int = (-0.5 * ((r_test[k][:, None] - c_grid) / sd_obs[k][:, None]) ** 2
               - np.log(sd_obs[k][:, None] * np.sqrt(2 * np.pi))).sum(0)
    log_int += -0.5 * ((c_grid - MU_Y[j]) / S_Y[j]) ** 2 - np.log(S_Y[j] * np.sqrt(2 * np.pi))
    brute = np.log(np.trapezoid(np.exp(log_int - log_int.max()), c_grid)) + log_int.max()
    print(f"{TEXTILES[i]:>6} at {ad_label(yr)}: closed form {textile_loglik_numpy(r_test, i, j):.6f}, "
          f"numerical {brute:.6f}")

shroud at AD 1300: closed form -67.950807, numerical -67.950807

  cope at AD 1295: closed form -81.156223, numerical -81.156223
 mummy at AD 50: closed form -65.171430, numerical -65.171430

The closed form matches numerical integration to six decimals. Before fitting, a prior predictive check in the units that matter: with HalfNormal(50)-year priors on $\omega$ and both $\tau$s, how far apart do three labs' means of one textile typically end up, and what chi-square would the 1989 test compute?

n_sim = 4000
omega_s, tau_s = np.abs(rng.normal(0, 50, n_sim)), np.abs(rng.normal(0, 50, n_sim))
lab_err = E_obs[0]                                   # shroud lab-mean errors
sim_means = (rng.normal(0, 1, (n_sim, 3)) * omega_s[:, None] + rng.normal(0, 1, (n_sim, 3)) * tau_s[:, None]
             + rng.normal(0, 1, (n_sim, 3)) * lab_err)
spread = sim_means.max(1) - sim_means.min(1)
T_prior = ward_wilson_T(sim_means, lab_err)
print(f"prior: spread of three lab means, median {np.median(spread):.0f} yr (90%: {np.quantile(spread, 0.05):.0f}-"
      f"{np.quantile(spread, 0.95):.0f}); chi-square > 6: {np.mean(T_prior > 6):.2f}; "
      f"observed shroud spread {M_obs[0].max() - M_obs[0].min():.0f} yr")

prior: spread of three lab means, median 97 yr (90%: 22-277); chi-square > 6: 0.61; observed shroud spread 103 yr

The prior allows anything from perfect agreement to labs 250 years apart; the observed spread on the shroud (about 100 radiocarbon years between Arizona and Oxford) is well inside it, and so is agreement. Now the model.

MU_G, S2_G = pt.as_tensor(MU_Y), pt.as_tensor(S_Y**2)

with pm.Model(coords={"lab": LABS, "textile": TEXTILES}) as lab_model:
    omega = pm.HalfNormal("omega", 50)
    delta = pm.Deterministic("delta", omega * pm.Normal("delta_z", 0, 1, dims="lab"), dims="lab")
    tau_shroud = pm.HalfNormal("tau_shroud", 50)
    tau_ctrl = pm.HalfNormal("tau_ctrl", 50)
    tau = pt.stack([tau_shroud, tau_ctrl, tau_ctrl, tau_ctrl])[:, None]
    eta = pm.Deterministic("eta", tau * pm.Normal("eta_z", 0, 1, dims=("textile", "lab")), dims=("textile", "lab"))
    r = y_obs - delta[lab_i] - eta[smp_i, lab_i]
    rbar = pm.Deterministic("rbar", pt.dot(M_ONE, w_obs * r) / W_I, dims="textile")
    within = pt.sum(-0.5 * w_obs * (r - rbar[smp_i]) ** 2)       # constants dropped
    v = 1 / W_I[:, None] + S2_G[None, :]
    log_grid = LOG_PRIOR - 0.5 * (rbar[:, None] - MU_G[None, :]) ** 2 / v - 0.5 * pt.log(v)
    pm.Potential("likelihood", within + pt.sum(pt.logsumexp(log_grid, axis=1)))
    lab_idata = pm.sample(random_seed=RANDOM_SEED, progressbar=False, target_accept=0.9)

print(f"divergences: {int(lab_idata.sample_stats['diverging'].sum())}, "
      f"tuning steps: {lab_idata.posterior.attrs.get('tuning_steps')}")
az.summary(lab_idata, var_names=["omega", "delta", "tau_shroud", "tau_ctrl"], round_to=2)

divergences: 0, tuning steps: 400

az.summary(lab_idata, var_names=["eta"], round_to=1).iloc[:, :4]

Fast (seconds), no divergences, r_hat at most 1.01. What the posterior says:

• The lab offsets are all slightly positive (the labs read a little old, anchored by the controls, most clearly the cope, whose measured age is older than the curve allows for 1290-1310), with Oxford highest, but every interval comfortably includes zero and the between-lab sd $\omega$ is poorly determined - three labs cannot pin down a population sd.
• Heterogeneity: $\tau_{\text{shroud}}$ is about twice $\tau_{\text{ctrl}}$ in the mean, and its interval is wide. The shroud-specific effects $\eta$ put Arizona low and Oxford high, the pattern in the raw data.

How strongly do the data say the shroud is more heterogeneous than the controls?

post_lab = az.extract(lab_idata)
p_tau = float((post_lab["tau_shroud"] > post_lab["tau_ctrl"]).mean())
print(f"P(tau_shroud > tau_ctrl) = {p_tau:.2f};  prior value 0.50")
print("tau_shroud quantiles (5/50/95%):", np.quantile(post_lab["tau_shroud"], [0.05, 0.5, 0.95]).round(0))
print("tau_ctrl   quantiles (5/50/95%):", np.quantile(post_lab["tau_ctrl"], [0.05, 0.5, 0.95]).round(0))

P(tau_shroud > tau_ctrl) = 0.75;  prior value 0.50
tau_shroud quantiles (5/50/95%): [ 6. 36. 87.]
tau_ctrl   quantiles (5/50/95%): [ 1. 17. 46.]

Moderately: the shroud looks more heterogeneous than the controls with probability around 3 in 4, not overwhelmingly. That is the honest reading of the 1989 chi-square debate: three lab means (and a dozen measurements) are too few to settle whether the shroud sample was inhomogeneous - Casabianca et al. (2019, Archaeometry) argued it was from the raw data - but the model absorbs the extra scatter into $\tau_{\text{shroud}}$, so the calendar date below already pays for it.

Posterior predictive check: the 1989 chi-square test. For each posterior draw, simulate a fresh set of 49 measurements (new $c_i$, the same labs' offsets and heterogeneity, the quoted measurement errors), recompute each textile's lab means and Ward & Wilson's chi-square, and see where the observed values fall.

def grid_weights(rbar_draws, i):
    """Posterior probability of each calendar year for textile i, one row per draw."""
    v = 1 / W_I[i] + S_Y**2
    lg = LOG_PRIOR[i] - 0.5 * (rbar_draws[:, None] - MU_Y) ** 2 / v - 0.5 * np.log(v)
    lg -= lg.max(1, keepdims=True)
    p = np.exp(lg)
    return p / p.sum(1, keepdims=True)


thin = az.extract(lab_idata, num_samples=1000, random_seed=RANDOM_SEED)
delta_d = thin["delta"].to_numpy().T                           # draws x lab
eta_d = thin["eta"].transpose("sample", "textile", "lab").to_numpy()
rbar_d = thin["rbar"].to_numpy().T                             # draws x textile
n_d = len(rbar_d)
theta_d = np.zeros((n_d, 4), int)
year_post = np.zeros((4, len(YEAR)))
for i in range(4):
    P = grid_weights(rbar_d[:, i], i)
    year_post[i] = P.mean(0)
    theta_d[:, i] = YEAR[(P.cumsum(1) < rng.uniform(size=(n_d, 1))).sum(1)]
year_idx_d = theta_d - YEAR[0]
c_rep = rng.normal(MU_Y[year_idx_d], S_Y[year_idx_d])
y_rep = (c_rep[:, smp_i] + delta_d[:, lab_i] + eta_d[:, smp_i, lab_i] + rng.normal(0, 1, (n_d, len(d))) * sd_obs)
T_rep = np.array([ward_wilson_T(*lab_means(yr, sd_obs, smp_i, lab_i)) for yr in y_rep])
ppp = (T_rep >= T_obs).mean(0)

print("posterior predictive P(replicated chi-square >= observed):", dict(zip(TEXTILES, ppp.round(2))))
print("classical p-values if the quoted errors were the whole story (chi-square, 2 df):",
      dict(zip(TEXTILES, np.exp(-T_obs / 2).round(3))))
fig, axes = plt.subplots(1, 4, figsize=(14, 3.4), sharey=True)
for i, ax in enumerate(axes):
    ax.hist(np.clip(T_rep[:, i], 0, 30), bins=np.arange(0, 30.5, 1), color=GREY, alpha=0.7, density=True)
    ax.axvline(T_obs[i], color=CAL, lw=2)
    ax.set_title(f"{TEXTILES[i]}: observed {T_obs[i]:.1f}, p = {ppp[i]:.2f}", fontsize=11)
    ax.set(xlabel="chi-square of lab means", yticks=[])
fig.suptitle("Posterior predictive check of the inter-laboratory chi-square (grey: replicates, 30 = 30+; "
             "red: observed)", fontsize=11);

posterior predictive P(replicated chi-square >= observed): {'shroud': np.float64(0.34), 'Nubia': np.float64(0.94), 'mummy': np.float64(0.5), 'cope': np.float64(0.41)}
classical p-values if the quoted errors were the whole story (chi-square, 2 df): {'shroud': np.float64(0.014), 'Nubia': np.float64(0.901), 'mummy': np.float64(0.28), 'cope': np.float64(0.139)}

With the heterogeneity terms in the model, the observed chi-squares are typical of what the model generates for every textile, the shroud included: the model has explained the disagreement rather than ignored it. (A model that trusted the quoted errors, $\tau = 0$, gives the shroud's chi-square a p-value of about 1% with these errors, 4% with the paper's - which is where the 1989 debate started - the second line above; the chi-square of 2 d.f. has the closed-form p-value $e^{-T/2}$.)

In plain words: The three labs disagree a bit more about the shroud than about the control cloths. Our model lets each lab have its own small bias and lets the shroud pieces differ from each other, and then asks what date is consistent with all of it.

#6 · The calendar date of the shroud

For every posterior draw the grid gives the exact conditional distribution of each textile's calendar year; averaging those over draws gives the marginal posterior (and one sampled year per draw gives draws of $\theta$ for anything else). First the controls, as a check.

def summarise_year(p, year=YEAR):
    c = np.cumsum(p)
    q = [int(year[np.searchsorted(c, v)]) for v in (0.025, 0.5, 0.975)]
    return {"median": ad_label(q[1]), "95% equal-tailed": f"{ad_label(q[0])} to {ad_label(q[2])}",
            "95% HPD": fmt_ranges(hpd_ranges(p, year))}


pd.DataFrame({t: summarise_year(year_post[i]) | {"prior window": f"{ad_label(WINDOWS[t][0])} to {ad_label(WINDOWS[t][1])}"}
              for i, t in enumerate(TEXTILES)}).T

The controls land where history puts them, and the mummy linen, which had only a loose prior, comes out in the 1st century BC to 2nd century AD, like the British Museum's own date. The cope's posterior piles up at the early end of its window: the labs measured it a little too old for 1290-1310, which is precisely the information that nudges the lab offsets up.

Now the shroud, compared with calibrating Damon et al.'s pooled 691 ± 31.

p_model = year_post[0]
prob = {
    "before AD 1260": (YEAR < 1260),
    "AD 1260-1390 (Damon et al.'s range)": (YEAR >= 1260) & (YEAR <= 1390),
    "after AD 1390": (YEAR > 1390),
    "before AD 1000": (YEAR < 1000),
    "before AD 400": (YEAR < 400),
}
pd.DataFrame({"hierarchical model": {k: p_model[m].sum() for k, m in prob.items()},
              "pooled 691 ± 31": {k: p_shroud_pooled[m].sum() for k, m in prob.items()}}).map(
    lambda v: f"{v:.4f}" if v >= 1e-4 else ("< 1e-4" if v > 0 else "0 (underflow)"))

fig, ax = plt.subplots(figsize=(11, 3.8))
m = (YEAR >= 1150) & (YEAR <= 1480)
ax.fill_between(YEAR[m], 0, p_model[m], color=CAL, alpha=0.45, lw=0, label="hierarchical model (lab offsets + heterogeneity)")
ax.plot(YEAR[m], p_shroud_pooled[m], color=INK, lw=1.5, ls="--", label="calibrating the pooled 691 ± 31")
for (a, b), yv in zip(hpd_ranges(p_model), [-0.0006] * 5):
    ax.plot([a, b], [yv, yv], color=CAL, lw=4, solid_capstyle="butt")
ax.set(xlabel="calendar year AD", ylabel="probability per year", title="When was the shroud's flax harvested?")
ax.legend(fontsize=9)
print("hierarchical model:", summarise_year(p_model))
print("pooled date       :", summarise_year(p_shroud_pooled));

hierarchical model: {'median': 'AD 1314', '95% equal-tailed': 'AD 1278 to AD 1389', '95% HPD': 'AD 1273-1326, AD 1353-1394'}
pooled date       : {'median': 'AD 1297', '95% equal-tailed': 'AD 1275 to AD 1385', '95% HPD': 'AD 1273-1317, AD 1360-1388'}

The hierarchical model moves the date slightly later than the pooled calibration (the labs' shared tendency to read old, learned from the controls, is subtracted) and gives the two bumps more equal weight; the chance of a date before 1260 stays well under 1%. The first-century question has an unambiguous answer within the model: the probability is not small but numerically zero, because a first-century linen should read about 1,950-2,000 radiocarbon years and the labs measured under 700, some forty standard errors away.

So the real question is outside the model: could something have made a first-century cloth look medieval? The usual candidate is contamination with younger carbon. For each posterior draw of the shroud's true radiocarbon content we can ask how much of the carbon would have to be younger contamination for a cloth from AD 30 to read as it did. Radiocarbon content is $F = e^{-\text{age}/8033}$ (the conventional Libby mean life), and a mixture with a fraction $f$ of contaminant gives $F_{\text{meas}} = (1-f)F_{\text{true}} + f F_{\text{cont}}$.

def fraction_needed(age_meas, age_true, age_cont):
    F = lambda a: np.exp(-np.asarray(a, float) / 8033)
    return (F(age_meas) - F(age_true)) / (F(age_cont) - F(age_true))


age_ad30 = float(np.interp(1950 - 30, curve.cal_bp, curve.c14_bp))
c_shroud = rng.normal(rbar_d[:, 0], 1 / np.sqrt(W_I[0]))        # the shroud's true 14C age, per draw
age_1532 = float(np.interp(1950 - 1532, curve.cal_bp, curve.c14_bp))
contam = {"carbon from 1950 (the 'youngest' possible before bomb tests)": fraction_needed(c_shroud, age_ad30, 0.0),
          "carbon from the 1532 Chambéry fire": fraction_needed(c_shroud, age_ad30, age_1532)}
print(f"radiocarbon age of AD 30 on IntCal20: {age_ad30:.0f} BP; of AD 1532: {age_1532:.0f} BP")
print(f"shroud's true radiocarbon age (lab effects removed): median {np.median(c_shroud):.0f} BP, "
      f"90% {np.quantile(c_shroud, 0.05):.0f}-{np.quantile(c_shroud, 0.95):.0f}")
for k, f in contam.items():
    print(f"fraction of the carbon that must be contamination, {k}: median {np.median(f):.0%} "
          f"(90%: {np.quantile(f, 0.05):.0%}-{np.quantile(f, 0.95):.0%})")

radiocarbon age of AD 30 on IntCal20: 1971 BP; of AD 1532: 299 BP
shroud's true radiocarbon age (lab effects removed): median 661 BP, 90% 607-710
fraction of the carbon that must be contamination, carbon from 1950 (the 'youngest' possible before bomb tests): median 64% (90%: 61%-67%)
fraction of the carbon that must be contamination, carbon from the 1532 Chambéry fire: median 77% (90%: 73%-80%)

For a first-century cloth to give these readings, about two thirds of its carbon would have to be twentieth-century contamination that survived three labs' different cleaning methods, and about three quarters if the contamination came from the 1532 fire. That is the scale of the claim a first-century date requires. What the radiocarbon data cannot rule out is that the strip was not representative of the cloth (the "medieval repair" hypothesis): that is a question about the sample, not the measurement, and no statistical model of these 49 numbers can answer it.

In plain words: The labs' readings, adjusted for each lab's small quirks, point to the late 1200s or the 1300s, with well under a 1-in-100 chance of anything before 1260. For the cloth to be from Jesus's time, well over half of the carbon in the tested strip would have to be modern dirt that three different cleaning methods all missed.

#7 · A Bayesian chronology: Sluggan Bog

A single date is often vague; a sequence of dates can be sharp, because the order of the layers is information. Our example is the peat core from Sluggan Bog, Co. Antrim, Northern Ireland, studied for its pollen record by Smith & Goddard (1991, New Phytologist 118) and used as a test case for Bayesian age-depth modelling in the R package Bchron (Haslett & Parnell 2008). 31 radiocarbon dates on bulk peat slices, from 44.5 cm to 518 cm depth. Several dates come from the same slice (the same depth, 2-5 cm thick). The file is an R data file; as in E23/E24 we read R's serialisation format directly.

def read_rdata(path):
    """Minimal reader for R's XDR serialisation (gzip, bzip2 or xz): named lists, vectors, strings."""
    raw = Path(path).read_bytes()
    buf = (bz2.decompress(raw) if raw[:2] == b"BZ" else
           lzma.decompress(raw) if raw[:6] == b"\xfd7zXZ\x00" else gzip.decompress(raw))
    assert buf[:7] in (b"RDX2\nX\n", b"RDX3\nX\n"), "not an XDR R data file"
    pos, refs = 7, []

    def i32():
        nonlocal pos
        pos += 4
        return struct.unpack(">i", buf[pos - 4:pos])[0]

    def item():
        nonlocal pos
        flags = i32()
        kind, has_attr, has_tag = flags & 0xFF, flags & (1 << 9), flags & (1 << 10)
        if kind == 254:                                     # NULL
            return None
        if kind == 255:                                     # reference to an earlier symbol
            return refs[(flags >> 8) - 1]
        if kind == 1:                                       # symbol
            refs.append(item())
            return refs[-1]
        if kind == 9:                                       # string
            n = i32()
            pos += max(n, 0)
            return None if n == -1 else buf[pos - n:pos].decode()
        if kind == 2:                                       # pairlist -> dict
            out = {}
            while kind == 2:
                _ = item() if has_attr else None
                tag = item() if has_tag else None
                out[tag] = item()
                flags = i32()
                kind, has_attr, has_tag = flags & 0xFF, flags & (1 << 9), flags & (1 << 10)
            return out
        n = i32()
        if kind in (10, 13):                                # logical / integer
            value, pos = np.frombuffer(buf, ">i4", n, pos).astype(int), pos + 4 * n
        elif kind == 14:                                    # double
            value, pos = np.frombuffer(buf, ">f8", n, pos).astype(float), pos + 8 * n
        elif kind in (16, 19):                              # character vector / list
            value = [item() for _ in range(n)]
        else:
            raise NotImplementedError(f"R type {kind}")
        attrs = item() if has_attr else None
        if kind == 19 and attrs and "names" in attrs:
            value = dict(zip(attrs["names"], value))
        return value

    version = i32()
    pos += 8
    if version == 3:
        pos += 4 + struct.unpack(">i", buf[pos:pos + 4])[0]
    return item()


data.describe("sluggan_moss")
bog = pd.DataFrame(read_rdata(data.path("sluggan_moss"))["Sluggan"])
bog = bog.drop(columns="calCurves").sort_values(["position", "ages"]).reset_index(drop=True)
DEPTHS = np.sort(bog.position.unique())
K = len(DEPTHS)
lev = np.searchsorted(DEPTHS, bog.position.to_numpy())
print(f"{len(bog)} dates at {K} depths; slices with several dates: "
      f"{ {float(x): int(n) for x, n in bog.position.value_counts().items() if n > 1} }")
bog.T

sluggan_moss: R data file (bzip2-compressed XDR): a data frame of 31 radiocarbon dates down a peat core: id, ages (14C yr BP), ageSds, position (depth in cm), thickness (cm), calCurves. Several dates share a depth.
  source: Smith & Goddard (1991), A 12,500 year record of vegetational history at Sluggan Bog, Co. Antrim, N. Ireland, New Phytologist 118:167-187; as distributed in the R package Bchron (Haslett & Parnell 2008), from the European Pollen Database.
  url:    https://raw.githubusercontent.com/cran/Bchron/master/data/Sluggan.rda
31 dates at 21 depths; slices with several dates: {232.5: 3, 272.5: 3, 367.5: 3, 447.5: 3, 484.5: 2, 518.0: 2}

Calibrating each date on its own (the same grid calculation, now in cal BP over 0-20,000 years at 5-year steps), and plotting against depth:

TOP, STEP = 20000.0, 5.0
BP = np.arange(0, TOP + 1, STEP)
MU_B = np.interp(BP, curve.cal_bp, curve.c14_bp)
S_B = np.interp(BP, curve.cal_bp, curve.c14_sd)
y_bog, e_bog = bog.ages.to_numpy(float), bog.ageSds.to_numpy(float)
V_B = e_bog[:, None] ** 2 + S_B**2
LIK_B = np.exp(-0.5 * (y_bog[:, None] - MU_B) ** 2 / V_B) / np.sqrt(2 * np.pi * V_B)  # density in y, per cal year
# cumulative integral over calendar years, scaled so each date's total is 1
CUM = np.concatenate([np.zeros((len(bog), 1)), np.cumsum((LIK_B[:, 1:] + LIK_B[:, :-1]) / 2 * STEP, axis=1)], axis=1)
CUM /= CUM[:, -1:]
indep_q = np.array([[np.interp(q, c, BP) for q in (0.025, 0.5, 0.975)] for c in CUM])
bog["indep_lo"], bog["indep_med"], bog["indep_hi"] = indep_q.T.round()


def depth_plot(ax, lo, med, hi, color, label, offset=0.0, lw=2):
    ax.plot(med, bog.position + offset, "o", color=color, ms=4, label=label)
    ax.hlines(bog.position + offset, lo, hi, color=color, lw=lw, alpha=0.8)


fig, ax = plt.subplots(figsize=(11, 5.5))
depth_plot(ax, *indep_q.T, GREY, "each date calibrated on its own (95%)")
ax.invert_yaxis()
ax.set(xlabel="calendar age (cal BP, years before 1950)", ylabel="depth (cm)",
       title="Sluggan Bog: 31 dates, each calibrated independently")
ax.legend(fontsize=9);
bog[["position", "ages", "ageSds", "indep_lo", "indep_med", "indep_hi"]].iloc[18:25]

Mostly the dates get older with depth, as they must - but not everywhere. The slice at 427 cm (9360 ± 150 BP) calibrates to about 10,600 cal BP, older than one of the three dates in the slice 20 cm below it (8895 ± 125, about 9,970 cal BP). And the three dates in the 367.5 cm slice span 8,830-9,370 cal BP between them. Bulk peat dates are notorious for this: roots grow down into older peat, and old carbon washes in.

#8 · A sequence of phases, with the dates integrated out

The OxCal-style model (Bronk Ramsey 2009) for a stratified sequence treats each slice as a phase bounded by two boundaries: every date in slice $k$ lies somewhere between boundary $B_{k-1}$ (the top of the slice) and $B_k$ (its bottom), uniformly, and the boundaries are in depth order, $B_0 < B_1 < \dots < B_K$ in cal BP. Dates within one slice are not ordered among themselves. $B_0$ is the end of the record (the top of the dated peat) and $B_K$ its start.

Priors. end $\sim U(0, 20000)$ and span $= B_K - B_0 \sim U(0, 20000)$ cal BP, and the 21 slice durations are the span times $\text{Dirichlet}(1, \dots, 1)$ proportions. Dirichlet(1) spacings are exactly the gaps of uniformly scattered boundaries given the span, so this is the uniform-phase prior - but because the span has its own prior, adding more boundaries does not push the total span up (the known bias of putting independent uniform priors on many ordered events; Nicholls & Jones 2001, Applied Statistics). This is the "ordered transform or equivalent": pm.Dirichlet + cumsum gives ordered boundaries by construction.

Integrating the dates out. Given the boundaries, date $i$ in slice $k$ contributes

$$p(y_i \mid B) = \frac{1}{B_k - B_{k-1}} \int_{B_{k-1}}^{B_k} \mathcal N\bigl(y_i \mid \mu(\theta),\ \sigma_i^2 + s(\theta)^2\bigr)\, d\theta ,$$

a difference of the date's cumulative calibrated likelihood, precomputed once on the grid (CUM above) and linearly interpolated at the boundaries. The integral is smooth in the boundaries - it averages the wiggles away - so NUTS sees 23 well-behaved parameters and no calendar dates at all. The dates are recovered afterwards by inverse-CDF sampling between their boundaries.

Outliers. OxCal's answer to dates like the 427 cm one is an outlier model (Bronk Ramsey 2009, Radiocarbon 51:1023): each date is, with prior probability 5%, unrelated to its slice. We make the outlier's calendar age uniform over the whole 0-20,000 range, so its likelihood is $\int \mathcal N(y_i \mid \mu(\theta), \dots)\, d\theta / 20000$, and mix: $p(y_i \mid B) = 0.95\, p_{\text{in}} + 0.05\, p_{\text{out}}$ (the indicator is summed out). First, the model without outliers.

CUM_T = pt.as_tensor(CUM)
LOG_OUT = np.log(1 / TOP)      # every date's total (scaled) mass is 1, spread over 20,000 years
IDX = np.arange(len(bog))


def cum_at(b):
    """Each date's cumulative calibrated likelihood at the boundaries b (dates x boundaries)."""
    x = pt.clip(b / STEP, 0, len(BP) - 1.000001)
    i0 = pt.cast(pt.floor(x), "int64")
    f = x - i0
    return CUM_T[:, i0] * (1 - f) + CUM_T[:, i0 + 1] * f


def sequence_model(p_outlier):
    with pm.Model(coords={"slice": DEPTHS, "boundary_id": np.arange(K + 1), "date": IDX}) as model:
        end = pm.Uniform("end", 0.0, TOP)
        span = pm.Uniform("span", 0.0, TOP)
        gaps = pm.Dirichlet("gaps", np.ones(K), dims="slice")
        b = pm.Deterministic("boundary", end + span * pt.concatenate([pt.zeros(1), pt.cumsum(gaps)]),
                             dims="boundary_id")
        cb = cum_at(b)
        mass = cb[IDX, lev + 1] - cb[IDX, lev]
        ll_in = pt.log(pt.maximum(mass, 1e-300)) - pt.log(b[lev + 1] - b[lev])
        ll = ll_in if p_outlier == 0 else pt.logaddexp(np.log1p(-p_outlier) + ll_in, np.log(p_outlier) + LOG_OUT)
        pm.Deterministic("loglik_date", ll, dims="date")
        pm.Deterministic("loglik_in", ll_in, dims="date")
        pm.Potential("likelihood", ll.sum())
    return model


# start: boundaries between the independent medians of neighbouring slices, forced into order
slice_med = np.maximum.accumulate([np.median(indep_q[lev == k, 1]) for k in range(K)]) + np.arange(K)
b_init = np.r_[max(slice_med[0] - 100, 1.0), (slice_med[1:] + slice_med[:-1]) / 2, slice_med[-1] + 300]
INITVALS = {"end": b_init[0], "span": b_init[-1] - b_init[0], "gaps": np.diff(b_init) / (b_init[-1] - b_init[0])}


def per_chain(idata):
    return pd.DataFrame({"mean log density": idata.sample_stats["logp"].mean("draw").to_numpy().round(1),
                         "span (median)": idata.posterior["span"].median("draw").to_numpy().round(0),
                         "start (median)": idata.posterior["boundary"].isel(boundary_id=-1).median("draw").to_numpy().round(0)},
                        index=[f"chain {c}" for c in range(idata.posterior.sizes["chain"])])


def diagnose(idata, name):
    ss = idata.sample_stats
    print(f"{name}: divergences {int(ss['diverging'].sum())}, mean tree depth {float(ss['depth'].mean()):.1f} "
          f"(max 10), boundaries: max r_hat {float(az.rhat(idata)['boundary'].max()):.2f}, "
          f"min bulk ESS {float(az.ess(idata)['boundary'].min()):.0f}")
    return per_chain(idata)


strict_idata = pm.sample(model=sequence_model(0.0), random_seed=RANDOM_SEED, progressbar=False,
                         initvals=INITVALS, target_accept=0.9, tune=1500,
                         var_names=["end", "span", "gaps", "boundary"])
diagnose(strict_idata, "no outliers")

no outliers: divergences 8, mean tree depth 10.0 (max 10), boundaries: max r_hat 1.02, min bulk ESS 255

All four chains agree (the per-chain log densities and start boundaries match, r_hat about 1.02), but with a few divergences and a mean tree depth at the maximum of 10. The strict model also has a trap that this run happened to avoid: while developing the notebook, with other seeds and a shorter warm-up, one chain in two runs out of three got stranded hundreds or thousands of log-density units below the others. nutpie's random start (jitter on the unconstrained scale) had put some date outside its slice, where its calibrated likelihood underflows to zero; the log density is then flat (we floor it at $10^{-300}$) and no gradient leads back. And the strict model must squeeze the conflicting 427 cm and 447.5 cm dates into order however badly they fit, which is exactly what an outlier model is for.

Now with a 5% outlier prior, same sampler settings:

seq_model = sequence_model(0.05)
diag_idata = pm.sample(model=seq_model, random_seed=RANDOM_SEED, progressbar=False, initvals=INITVALS,
                       target_accept=0.9, tune=1500, var_names=["end", "span", "gaps", "boundary"])
diagnose(diag_idata, "outliers, diagonal mass matrix")

outliers, diagonal mass matrix: divergences 0, mean tree depth 10.0 (max 10), boundaries: max r_hat 1.25, min bulk ESS 12

Worse, not better: one chain disagrees with the others about the span and the start, r_hat climbs to about 1.25 and the smallest ESS is around a dozen, without a single divergence.

The tell-tale sign is the tree depth: nearly every iteration hits the maximum of 10 (1,023 leapfrog steps). The posterior is long and thin in directions that mix many parameters - move one interior boundary and its neighbours, the span and the Dirichlet proportions all have to move together - and a diagonal mass matrix cannot rescale a diagonal ridge. nutpie's low-rank mass-matrix adaptation (adaptation="low_rank") learns the few dominant correlated directions during warm-up:

seq_idata = pm.sample(model=seq_model, random_seed=RANDOM_SEED, progressbar=False, initvals=INITVALS,
                      target_accept=0.9, tune=1500, nuts={"adaptation": "low_rank"})
del diag_idata
display(diagnose(seq_idata, "outliers, low-rank mass matrix"))
az.summary(seq_idata, var_names=["end", "span"], round_to=0)

outliers, low-rank mass matrix: divergences 1, mean tree depth 9.2 (max 10), boundaries: max r_hat 1.01, min bulk ESS 468

Every chain now agrees, r_hat is at most about 1.01 and the smallest bulk ESS over the 22 boundaries is in the hundreds, at a similar run time. The outlier component does two jobs: statistically it lets a conflicting date be discounted, and computationally it puts a floor under every date's likelihood ($0.05 / 20000$ per year), so there are no flat zero-likelihood plateaus to get stranded on.

Which dates are outliers? The posterior probability that date $i$ is an outlier, averaged over draws, is $0.05\,p_{\text{out}} / (0.05\,p_{\text{out}} + 0.95\,p_{\text{in}})$.

post_seq = az.extract(seq_idata)
ll_in_d = post_seq["loglik_in"].to_numpy().T                    # draws x dates
log_in = np.log(0.95) + ll_in_d
log_out = np.log(0.05) + LOG_OUT
p_out_d = 1 / (1 + np.exp(log_in - log_out))
bog["p_outlier"] = p_out_d.mean(0).round(3)
print("prior outlier probability 0.05; the six highest posterior values:")
bog.sort_values("p_outlier", ascending=False).head(6)[["position", "ages", "ageSds", "indep_med", "p_outlier"]]

prior outlier probability 0.05; the six highest posterior values: