Does pypomp’s ABC-MCMC sample the right ABC posterior at each tolerance?

Published

October 16, 2026

Show Code 
import os
import sys

sys.path.append("..")
sys.path.append("../..")
import checks
import report_utils as ru
from IPython.display import display, HTML, Markdown

display(HTML(f"<div style='margin-bottom: 20px;'>{ru.nav_bar('abc')}</div>"))
PMCMC  |  ABC

Tests Pomp._abc at each tolerance in the ladder ε = 10⁶, 4, 2.5, 1.5 against the exact ABC posterior from rejection sampling (reference/run.py) and against R pomp::abc. Data: 208 weeks simulated from pomp::sir at known parameters. Free parameters beta1 and rho have a flat prior on model.PRIOR_BOX, and all others are fixed at truth. Probes: mean reports in each quarter of the year.

Verdicts use z = difference / Monte Carlo SE, with SEs from bulk ESS: PASS |z| < 3, CHECK 3 ≤ |z| < 5, FAIL |z| ≥ 5. Draws are pooled over chains after discarding the first half of each chain.

Show Code 
import glob
import re

import numpy as np
import pandas as pd
from plotnine import (
    aes,
    geom_density,
    geom_line,
    geom_point,
    geom_vline,
    ggplot,
    labs,
    scale_x_log10,
)

os.environ["JAX_PLATFORMS"] = "cpu"
import model

GPU = os.path.join("results", "gpu")
R_DIR = os.path.join("results", "R")
REF_DIR = os.path.join("..", "reference", "results", "gpu")
PARAMS = list(model.FREE)
OTHER = {"beta1": "rho", "rho": "beta1"}
BURN_FRAC = 0.5
PRIOR_EPS = max(model.ABC_EPS_LADDER)

truth = model.true_theta()

sweeps = {}
for path in sorted(glob.glob(os.path.join(GPU, "eps*"))):
    m = re.search(r"eps([0-9.e+-]+)$", path)
    tr = ru.load_mcmc_traces(os.path.join(path, "traces.csv.gz"), "pypomp")
    if m and tr is not None:
        sweeps[float(m.group(1))] = {
            "traces": tr,
            "post": ru.drop_burnin(tr, BURN_FRAC),
            "acceptance": ru.read_if_exists(os.path.join(path, "acceptance.csv")),
            "meta": ru.load_json(os.path.join(path, "latest.json")),
        }

r_traces = ru.load_mcmc_traces(os.path.join(R_DIR, "traces.csv.gz"), "R pomp")
r_meta = ru.load_json(os.path.join(R_DIR, "latest.json"))
r_acc = ru.read_if_exists(os.path.join(R_DIR, "acceptance.csv"))
r_sweeps = {}
if r_traces is not None:
    for eps, tr in r_traces.groupby("epsilon"):
        r_sweeps[float(eps)] = {"traces": tr, "post": ru.drop_burnin(tr, BURN_FRAC)}

rejection = ru.read_if_exists(os.path.join(REF_DIR, "abc_rejection.csv.gz"))
ref_meta = ru.load_json(os.path.join(REF_DIR, "latest.json"))
ref_cfg = ref_meta.get("run_config", {}) if ref_meta else {}
grid = ru.read_if_exists(os.path.join(REF_DIR, "grid_posterior.csv"))

have_py = len(sweeps) > 0
have_r = len(r_sweeps) > 0
have_ref = rejection is not None
Show Code 
scorecard = []


def ref_moments(eps):
    if eps > ref_cfg.get("EPS_KEEP", np.inf):
        return {p: ru.uniform_moments(*model.PRIOR_BOX[p]) for p in PARAMS}
    acc = rejection[rejection["distance"] < eps**2]
    return {p: ru.iid_moments(acc[p]) for p in PARAMS}


def comparison(a, b, la, lb):
    return pd.DataFrame([r for p in PARAMS for r in ru.compare(a[p], b[p], p, la, lb)])


def score(frames, check, question):
    df = pd.concat(frames, ignore_index=True)
    worst = ru.worst_row(df, "z")
    if worst is None:
        detail = "no standard errors (run level too low)"
    else:
        eps_txt = f" at ε={worst['epsilon']:g}" if "epsilon" in df.columns else ""
        detail = f"largest gap: {worst['statistic']} {worst['parameter']}{eps_txt}, {worst['z']:+.1f} SE"
    scorecard.append({
        "check": check, "question": question,
        "verdict": ru.worst_verdict(df["verdict"]), "detail": detail,
    })


py_probes = ru.read_if_exists(os.path.join(GPU, "probe_values.csv"))
r_probes = ru.read_if_exists(os.path.join(R_DIR, "probe_values.csv"))
probe_cmp = None
if py_probes is not None and r_probes is not None:
    probe_cmp = py_probes.merge(r_probes, on="probe", suffixes=(" pypomp", " R"))
    probe_cmp["relative difference"] = (
        (probe_cmp["value pypomp"] - probe_cmp["value R"]).abs() / probe_cmp["value R"].abs()
    )
    worst = probe_cmp["relative difference"].max()
    scorecard.append({
        "check": "1. Same probes",
        "question": "Do pypomp and R compute the same summary statistics on the data?",
        "verdict": "PASS" if worst < 1e-6 else "FAIL",
        "detail": f"largest relative difference {worst:.1e}",
    })

py_m = {e: {p: ru.mcmc_moments(s["post"], p) for p in PARAMS} for e, s in sweeps.items()}
r_m = {e: {p: ru.mcmc_moments(s["post"], p) for p in PARAMS} for e, s in r_sweeps.items()}

vs_ref = {}
if have_py and have_ref:
    for eps in sorted(sweeps, reverse=True):
        vs_ref[eps] = comparison(py_m[eps], ref_moments(eps), "pypomp", "exact").assign(epsilon=eps)
    tight = [vs_ref[e] for e in vs_ref if e != PRIOR_EPS]
    if tight:
        score(tight, "2. Right ABC posterior",
              "At each finite ε, does pypomp match rejection sampling?")

vs_R = {}
for eps in sorted(set(sweeps) & set(r_sweeps), reverse=True):
    vs_R[eps] = comparison(py_m[eps], r_m[eps], "pypomp", "R pomp").assign(epsilon=eps)
if vs_R:
    score(list(vs_R.values()), "3. Agrees with R", "At each ε, does pypomp match R pomp's abc?")

prior_cmp = None
if PRIOR_EPS in sweeps:
    prior_cmp = comparison(
        py_m[PRIOR_EPS], {p: ru.uniform_moments(*model.PRIOR_BOX[p]) for p in PARAMS},
        "pypomp", "prior"
    ).assign(epsilon=PRIOR_EPS)
    score([prior_cmp], "4. Huge ε gives the prior",
          f"At ε={PRIOR_EPS:g}, where every move is accepted, is the answer the flat prior?")

stuck_rows = []
for label, acc in [("pypomp", s["acceptance"].assign(epsilon=e)) for e, s in sweeps.items()
                   if s["acceptance"] is not None] + ([("R pomp", r_acc)] if r_acc is not None else []):
    for eps, g in acc.groupby("epsilon"):
        stuck_rows.append({"source": label, "epsilon": float(eps), "chains": len(g),
                           "mean acceptance": g["acceptance_rate"].mean(),
                           "chains that never moved": int((g["acceptance_rate"] == 0).sum())})
stuck = pd.DataFrame(stuck_rows)
if len(stuck):
    stuck = stuck.sort_values(["source", "epsilon"], ascending=[True, False])
    frac = stuck["chains that never moved"].sum() / stuck["chains"].sum()
    scorecard.append({
        "check": "5. No stuck chains",
        "question": "Did every chain accept at least one move at every ε?",
        "verdict": "PASS" if frac == 0 else ("CHECK" if frac < 0.01 else "FAIL"),
        "detail": f"{int(stuck['chains that never moved'].sum())} of {int(stuck['chains'].sum())} chain-runs never moved",
    })

conv_rows = []
for label, arms in (("pypomp", sweeps), ("R pomp", r_sweeps)):
    for eps, s in sorted(arms.items(), reverse=True):
        for p in PARAMS:
            arr = ru.chains_array(s["post"], p)
            conv_rows.append({"source": label, "epsilon": eps, "parameter": p,
                              "chains": arr.shape[0], "draws per chain": arr.shape[1],
                              "split R-hat": ru.split_rhat(arr), "bulk ESS": ru.bulk_ess(arr)})
conv = pd.DataFrame(conv_rows)
row = ru.ess_check(conv, "6. Enough draws")
if row:
    scorecard.append(row)

Summary

Show Code 
for label, ok, path in (
    ("pypomp ABC runs", have_py, os.path.join(GPU, "eps*")),
    ("R pomp traces", have_r, R_DIR),
    ("Rejection reference", have_ref, os.path.join(REF_DIR, "abc_rejection.csv.gz")),
):
    if not ok:
        display(HTML(ru.missing_note(label, path)))
if scorecard:
    display(HTML(ru.scorecard_html(scorecard)))
py_meta = sweeps[min(sweeps)]["meta"] if have_py else {}
checks.save(
    [{"check": r["check"], "compared": r["question"], "result": r["detail"], "verdict": r["verdict"]}
     for r in scorecard],
    (py_meta.get("run_config") or {}).get("RUN_LEVEL"),
    pypomp_version=py_meta.get("pypomp_version"),
)
check question verdict detail
1. Same probes Do pypomp and R compute the same summary statistics on the data? PASS largest relative difference 3.6e-08
2. Right ABC posterior At each finite ε, does pypomp match rejection sampling? PASS largest gap: mean beta1 at ε=4, +1.2 SE
3. Agrees with R At each ε, does pypomp match R pomp's abc? PASS largest gap: mean rho at ε=2.5, -2.2 SE
4. Huge ε gives the prior At ε=1e+06, where every move is accepted, is the answer the flat prior? PASS largest gap: sd beta1 at ε=1e+06, +0.6 SE
5. No stuck chains Did every chain accept at least one move at every ε? PASS 0 of 4240 chain-runs never moved
6. Enough draws Does every run have at least 400 effective independent draws? PASS smallest bulk ESS 4,870 (R pomp, ε=1.5, beta1); split R-hat up to 1.136, see diagnostics

Visual review

Tolerance ladder

The ladder runs loosest first. ε = 10⁶ starts from the prior. Each tighter ε starts from points drawn from the previous run whose single simulation already falls within the new ε (ABC-SMC selection), so chains do not start where they can never move.

1. Same probes?

Probes from model.py and model.R on the data must agree to 1e-6 relative error.

Show Code 
if probe_cmp is not None:
    display(HTML(ru.table_html(probe_cmp, fmt="%.6g")))
else:
    display(Markdown("_Needs both run.py and run.R results._"))
probe value pypomp value R relative difference
qtr1 352.5 352.5 0
qtr2 700.981 700.981 6.69778e-09
qtr3 542.519 542.519 8.65409e-09
qtr4 194.846 194.846 3.6144e-08

2. Right ABC posterior?

pypomp vs rejection sampling at each finite ε.

Show Code 
if have_py and have_ref:
    eps_finite = [e for e in sorted(sweeps, reverse=True) if e != PRIOR_EPS]
    for param in PARAMS:
        frames = []
        for eps in eps_finite:
            frames.append(ru.thin_for_plot(sweeps[eps]["post"][[param]]).assign(
                epsilon=f"{eps:g}", source="pypomp ABC-MCMC"))
            acc = rejection[rejection["distance"] < eps**2]
            frames.append(ru.thin_for_plot(acc[[param]]).assign(
                epsilon=f"{eps:g}", source="exact (rejection)"))
        stacked = pd.concat(frames)
        display(
            ggplot(stacked, aes(x=param, color="epsilon", linetype="source"))
            + geom_density(size=0.9)
            + geom_vline(xintercept=truth[param], color=ru.color_palette["truth"], linetype="dotted")
            + labs(title=f"ABC-MCMC against the exact ABC posterior: {param}",
                   subtitle="Solid and dashed curves of the same color should coincide")
            + ru.theme_premium
        )
    display(HTML(ru.table_html(pd.concat([vs_ref[e] for e in eps_finite], ignore_index=True))))

parameter statistic pypomp exact difference SE z verdict epsilon
beta1 mean 394.8691 394.6441 0.2250 0.1918 1.1728 PASS 4.0000
beta1 sd 35.7421 35.7578 -0.0157 0.1357 -0.1160 PASS 4.0000
rho mean 0.6242 0.6243 -0.0000 0.0003 -0.1484 PASS 4.0000
rho sd 0.0571 0.0569 0.0002 0.0002 0.8638 PASS 4.0000
beta1 mean 393.2944 393.1643 0.1302 0.2056 0.6331 PASS 2.5000
beta1 sd 27.3432 27.1981 0.1452 0.1454 0.9986 PASS 2.5000
rho mean 0.6276 0.6280 -0.0003 0.0003 -1.0392 PASS 2.5000
rho sd 0.0406 0.0404 0.0002 0.0002 0.8493 PASS 2.5000
beta1 mean 392.3623 392.1115 0.2509 0.3254 0.7710 PASS 1.5000
beta1 sd 23.1945 23.2883 -0.0938 0.2301 -0.4076 PASS 1.5000
rho mean 0.6305 0.6303 0.0002 0.0004 0.4367 PASS 1.5000
rho sd 0.0323 0.0322 0.0001 0.0003 0.2929 PASS 1.5000

3. Agrees with R?

pypomp vs R pomp::abc with the same probes, scales, tolerances, prior and data. pomp’s random walk is on the natural scale, and pypomp’s is on log and logit scales, so acceptance rates are not comparable. The posteriors should match.

Show Code 
if vs_R:
    eps_shared = [e for e in vs_R if e != PRIOR_EPS]
    for param in PARAMS:
        frames = []
        for eps in eps_shared:
            frames.append(ru.thin_for_plot(sweeps[eps]["post"][[param]]).assign(epsilon=f"{eps:g}", source="pypomp"))
            frames.append(ru.thin_for_plot(r_sweeps[eps]["post"][[param]]).assign(epsilon=f"{eps:g}", source="R pomp"))
        display(
            ggplot(pd.concat(frames), aes(x=param, color="epsilon", linetype="source"))
            + geom_density(size=0.9)
            + geom_vline(xintercept=truth[param], color=ru.color_palette["truth"], linetype="dotted")
            + labs(title=f"pypomp against R pomp: {param}",
                   subtitle="Solid and dashed curves of the same color should coincide")
            + ru.theme_premium
        )
    display(HTML(ru.table_html(pd.concat(list(vs_R.values()), ignore_index=True))))
else:
    display(Markdown("_Needs both run.py and run.R results._"))

parameter statistic pypomp R pomp difference SE z verdict epsilon
beta1 mean 399.9755 399.5035 0.4720 0.5819 0.8111 PASS 1000000.0000
beta1 sd 57.9278 57.6501 0.2776 0.4115 0.6747 PASS 1000000.0000
rho mean 0.6252 0.6253 -0.0001 0.0013 -0.0403 PASS 1000000.0000
rho sd 0.1007 0.1009 -0.0001 0.0009 -0.1492 PASS 1000000.0000
beta1 mean 394.8691 394.8728 -0.0037 0.3173 -0.0116 PASS 4.0000
beta1 sd 35.7421 35.8644 -0.1223 0.2243 -0.5452 PASS 4.0000
rho mean 0.6242 0.6244 -0.0002 0.0005 -0.3061 PASS 4.0000
rho sd 0.0571 0.0573 -0.0002 0.0004 -0.6023 PASS 4.0000
beta1 mean 393.2944 392.7693 0.5251 0.2993 1.7544 PASS 2.5000
beta1 sd 27.3432 26.8933 0.4500 0.2116 2.1260 PASS 2.5000
rho mean 0.6276 0.6286 -0.0010 0.0005 -2.1996 PASS 2.5000
rho sd 0.0406 0.0404 0.0002 0.0003 0.6589 PASS 2.5000
beta1 mean 392.3623 392.1509 0.2114 0.4580 0.4616 PASS 1.5000
beta1 sd 23.1945 22.7306 0.4639 0.3238 1.4325 PASS 1.5000
rho mean 0.6305 0.6306 -0.0002 0.0006 -0.3058 PASS 1.5000
rho sd 0.0323 0.0320 0.0003 0.0004 0.8024 PASS 1.5000

4. Does a huge tolerance give the prior?

At ε = 10⁶, every move is accepted, so draws must be uniform on the prior box.

Show Code 
if prior_cmp is not None:
    display(HTML(ru.table_html(prior_cmp)))
parameter statistic pypomp prior difference SE z verdict epsilon
beta1 mean 399.9755 400.0000 -0.0245 0.4293 -0.0570 PASS 1000000.0000
beta1 sd 57.9278 57.7350 0.1927 0.3036 0.6349 PASS 1000000.0000
rho mean 0.6252 0.6250 0.0002 0.0010 0.1935 PASS 1000000.0000
rho sd 0.1007 0.1010 -0.0003 0.0007 -0.3907 PASS 1000000.0000

5. How the ABC posterior changes with ε

Acceptance and posterior sd by ε. The last row is the full posterior from the grid. Not scored.

Show Code 
if have_py:
    rows = []
    for eps, s in sorted(sweeps.items(), reverse=True):
        acc = s["acceptance"]
        row = {"epsilon": eps,
               "acceptance (pypomp)": acc["acceptance_rate"].mean() if acc is not None else np.nan}
        for p in PARAMS:
            row[f"sd {p} (pypomp)"] = py_m[eps][p]["sd"]
            if have_ref:
                row[f"sd {p} (exact)"] = ref_moments(eps)[p]["sd"]
        rows.append(row)
    if grid is not None:
        row = {"epsilon": 0.0}
        for p in PARAMS:
            row[f"sd {p} (exact)"] = ru.grid_moments(grid, p, OTHER[p])["sd"]
        rows.append(row)
    width = pd.DataFrame(rows)
    shown = width.copy()
    shown["epsilon"] = shown["epsilon"].map(lambda e: "full posterior (grid)" if e == 0 else f"{e:g}")
    display(HTML(ru.table_html(shown, fmt="%.4g")))

    finite = width[(width["epsilon"] > 0) & (width["epsilon"] < PRIOR_EPS)]
    if len(finite) > 1:
        display(
            ggplot(finite, aes(x="epsilon", y="acceptance (pypomp)"))
            + geom_line(color=ru.color_palette["pypomp"], size=1)
            + geom_point(color=ru.color_palette["pypomp"], size=2)
            + scale_x_log10()
            + labs(title="ABC-MCMC acceptance rate against tolerance",
                   subtitle="Tighter tolerance: sharper answer, more rejected moves",
                   x="epsilon (log scale)", y="acceptance rate")
            + ru.theme_premium
        )
epsilon acceptance (pypomp) sd beta1 (pypomp) sd beta1 (exact) sd rho (pypomp) sd rho (exact)
1e+06 0.8156 57.93 57.74 0.1007 0.101
4 0.4988 35.74 35.76 0.05711 0.05691
2.5 0.2422 27.34 27.2 0.04059 0.0404
1.5 0.06201 23.19 23.29 0.03231 0.03221
full posterior (grid) — — 9.885 — 0.01938

6. Did the chains mix?

The first table counts chains that never accepted a move. The second requires bulk ESS ≥ 400 per run. Split R-hat is for information only.

Show Code 
if len(stuck):
    display(HTML(ru.table_html(stuck, fmt="%.3g")))
source epsilon chains mean acceptance chains that never moved
R pomp 1e+06 36 0.83 0
R pomp 4 36 0.508 0
R pomp 2.5 36 0.245 0
R pomp 1.5 36 0.062 0
pypomp 1e+06 1024 0.816 0
pypomp 4 1024 0.499 0
pypomp 2.5 1024 0.242 0
pypomp 1.5 1024 0.062 0
Show Code 
if len(conv):
    display(HTML(ru.table_html(conv, fmt="%.3f")))
source epsilon parameter chains draws per chain split R-hat bulk ESS
pypomp 1000000.000 beta1 1024 50 1.047 18208.815
pypomp 1000000.000 rho 1024 50 1.079 9332.644
pypomp 4.000 beta1 1024 100 1.023 34991.358
pypomp 4.000 rho 1024 100 1.028 30288.307
pypomp 2.500 beta1 1024 100 1.042 17939.306
pypomp 2.500 rho 1024 100 1.041 17132.963
pypomp 1.500 beta1 1024 100 1.136 5189.937
pypomp 1.500 rho 1024 100 1.130 5361.259
R pomp 1000000.000 beta1 36 1000 1.001 21538.410
R pomp 1000000.000 rho 36 1000 1.001 16905.574
R pomp 4.000 beta1 36 1000 1.000 20054.242
R pomp 4.000 rho 36 1000 1.001 19333.220
R pomp 2.500 beta1 36 1000 1.001 15094.724
R pomp 2.500 rho 36 1000 1.001 14331.295
R pomp 1.500 beta1 36 1000 1.006 4869.868
R pomp 1.500 rho 36 1000 1.004 5844.363
Show Code 
if have_py:
    eps = min(sweeps)
    tr = sweeps[eps]["traces"]
    shown = tr[tr["chain"] < 16]
    for param in PARAMS:
        display(
            ggplot(shown, aes(x="iter", y=param, group="chain", color="factor(chain)"))
            + geom_line(alpha=0.7, size=0.4)
            + geom_vline(xintercept=tr["iter"].max() * BURN_FRAC, linetype="dotted", color="#7f8c8d")
            + labs(title=f"First 16 chains at ε={eps:g}: {param}",
                   subtitle="Left of the dotted line is discarded as burn-in",
                   x="iteration", color="chain")
            + ru.theme_premium
        )

Appendix: run provenance

Show Code 
if have_py:
    display(Markdown("**pypomp run**\n\n" + ru.format_metadata(py_meta)))
if have_r:
    display(Markdown("**R pomp run**\n\n" + ru.format_r_metadata(r_meta)))
if ref_meta:
    display(Markdown("**Rejection reference run**\n\n" + ru.format_metadata(ref_meta)))

pypomp run

  • Run at: 2026-10-08T00:41:13+00:00
  • pypomp: 1.0.8 / JAX: 0.11.2
  • Commit: c6ac80ca36
  • Devices: [‘cuda:0’]
  • GPU: NVIDIA RTX PRO 6000 Blackwell Server Edition
  • Config: kind=abc, model=sir, RUN_LEVEL=4, USE_CPU=False, MAIN_SEED=3141593, NCHAINS=1024, M=2000, epsilon=1.5, rung=3, probes=[‘qtr1’, ‘qtr2’, ‘qtr3’, ‘qtr4’], probe_scale={‘qtr1’: 25.4946033314898, ‘qtr2’: 42.5390499244671, ‘qtr3’: 42.4112568674804, ‘qtr4’: 24.9736557636163}, free_params=[‘beta1’, ‘rho’], rw_sd_estimation_scale={‘beta1’: 0.06, ‘rho’: 0.15}, prior_box={‘beta1’: [300.0, 500.0], ‘rho’: [0.45, 0.8]}, execution_time=357.2224667072296, platform=gpu, trace_thin=10

R pomp run

  • Run at: 2026-10-06T03:57:18
  • pomp: 6.3 / R: 4.4.0
  • Commit: f4bad06939
  • CPU: Intel(R) Xeon(R) Gold 6154 CPU @ 3.00GHz (36 cores)
  • Config: kind=abc, model=sir, RUN_LEVEL=4, NCHAINS=36, Nabc=40000, epsilon_ladder=[1000000, 4, 2.5, 1.5], trace_thin=20, probes=[‘qtr1’, ‘qtr2’, ‘qtr3’, ‘qtr4’], probe_scale={‘qtr1’: 25.4946, ‘qtr2’: 42.539, ‘qtr3’: 42.4113, ‘qtr4’: 24.9737}, free_params=[‘beta1’, ‘rho’], rw_sd_natural_scale={‘beta1’: 24, ‘rho’: 0.036}, execution_time=572.717

Rejection reference run

  • Run at: 2026-10-08T00:18:45+00:00
  • pypomp: 1.0.8 / JAX: 0.11.2
  • Commit: c6ac80ca36
  • Devices: [‘cuda:0’]
  • GPU: NVIDIA RTX PRO 6000 Blackwell Server Edition
  • Config: kind=reference, model=sir, RUN_LEVEL=4, USE_CPU=False, MAIN_SEED=3141593, GRID_N=60, NP_REF=5000, NREPS_REF=3, CHUNK=200, N_REJ=20000000, EPS_KEEP=4.0, free_params=[‘beta1’, ‘rho’], prior_box={‘beta1’: [300.0, 500.0], ‘rho’: [0.45, 0.8]}, grid_time=145.90228652954102, rejection_time=75.43277430534363, execution_time=221.33506083488464, platform=gpu