Does pypomp’s particle MCMC sample the right posterior?

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('pmcmc')}</div>"))
PMCMC  |  ABC

Tests Pomp._pmcmc against a grid posterior (reference/run.py) and R pomp::pmcmc. 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.

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_vline,
    ggplot,
    labs,
)

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

GPU = os.path.join("results", "gpu")
R_DIR = os.path.join("results", "R")
REF_DIR = os.path.join("..", "reference", "results", "gpu")
REFERENCE = os.path.join(REF_DIR, "grid_posterior.csv")
PARAMS = ["beta1", "rho"]
OTHER = {"beta1": "rho", "rho": "beta1"}
BURN_FRAC = 0.5

truth = pd.read_csv("../data/true_theta.csv")
truth = dict(zip(truth["parameter"], truth["value"]))

sweeps = {}
for path in sorted(glob.glob(os.path.join(GPU, "J*"))):
    m = re.search(r"J(\d+)$", path)
    tr = ru.load_mcmc_traces(os.path.join(path, "traces.csv.gz"), "pypomp")
    if m and tr is not None:
        sweeps[int(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_post = ru.drop_burnin(r_traces, BURN_FRAC) if r_traces is not None else None
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"))
grid = ru.read_if_exists(REFERENCE)
ref_meta = ru.load_json(os.path.join(REF_DIR, "latest.json"))
py_pf = ru.read_if_exists(os.path.join(GPU, "pfilter_logliks.csv"))
r_pf = ru.read_if_exists(os.path.join(R_DIR, "pfilter_logliks.csv"))
pf_vs_J = ru.read_if_exists(os.path.join(GPU, "pfilter_vs_J.csv"))

have_py = len(sweeps) > 0
have_r = r_post is not None
have_grid = grid is not None
J_HI = max(sweeps) if have_py else None
J_LO = min(sweeps) if have_py else None
Show Code 
# Everything the scorecard summarizes is computed here; later sections display it.
scorecard = []

precond = None
if py_pf is not None and r_pf is not None:
    a, b = py_pf["logLik"], r_pf["logLik"]
    se = float(np.hypot(a.std(ddof=1) / np.sqrt(len(a)), b.std(ddof=1) / np.sqrt(len(b))))
    z = (a.mean() - b.mean()) / se if se > 0 else float("nan")
    precond = {"pypomp": a.mean(), "R pomp": b.mean(), "difference": a.mean() - b.mean(),
               "SE": se, "z": z, "verdict": ru.verdict(z)}
    scorecard.append({
        "check": "1. Same model",
        "question": "Do pypomp and R give the same likelihood at the true parameters?",
        "verdict": precond["verdict"],
        "detail": f"logLik differs by {precond['difference']:.2f} ({z:+.1f} SE)",
    })

py_m = {J: {p: ru.mcmc_moments(s["post"], p) for p in PARAMS} for J, s in sweeps.items()}
grid_m = {p: ru.grid_moments(grid, p, OTHER[p]) for p in PARAMS} if have_grid else None
r_m = {p: ru.mcmc_moments(r_post, p) for p in PARAMS} if have_r else None


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


def score(df, check, question):
    worst = ru.worst_row(df, "z")
    detail = ("no standard errors (run level too low)" if worst is None else
              f"largest gap: {worst['statistic']} {worst['parameter']}, {worst['z']:+.1f} SE")
    scorecard.append({
        "check": check, "question": question,
        "verdict": ru.worst_verdict(df["verdict"]), "detail": detail,
    })


vs_grid = vs_J = vs_R = None
if have_py and have_grid:
    vs_grid = {J: comparison(py_m[J], grid_m, f"pypomp J={J}", "grid") for J in sweeps}
    score(vs_grid[J_HI], "2. Right posterior",
          f"Does pypomp (J={J_HI}) match the grid reference?")
if have_py and len(sweeps) > 1:
    vs_J = comparison(py_m[J_LO], py_m[J_HI], f"J={J_LO}", f"J={J_HI}")
    score(vs_J, "3. Particle count irrelevant",
          f"Do J={J_LO} and J={J_HI} give the same posterior?")
if have_py and have_r:
    vs_R = comparison(py_m[J_HI], r_m, f"pypomp J={J_HI}", "R pomp")
    score(vs_R, "4. Agrees with R", "Does pypomp match R pomp's pmcmc?")

conv_rows = []
for label, post in [(f"pypomp J={J}", s["post"]) for J, s in sorted(sweeps.items())] + (
    [("R pomp", r_post)] if have_r else []
):
    for p in PARAMS:
        arr = ru.chains_array(post, p)
        conv_rows.append({"source": label, "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, "5. Enough draws")
if row:
    scorecard.append(row)

Summary

Show Code 
for label, ok, path in (
    ("pypomp PMCMC runs", have_py, os.path.join(GPU, "J*")),
    ("R pomp traces", have_r, R_DIR),
    ("Grid reference", have_grid, REFERENCE),
):
    if not ok:
        display(HTML(ru.missing_note(label, path)))
if scorecard:
    display(HTML(ru.scorecard_html(scorecard)))
py_meta = sweeps[J_HI]["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 model Do pypomp and R give the same likelihood at the true parameters? PASS logLik differs by -0.01 (-2.2 SE)
2. Right posterior Does pypomp (J=2000) match the grid reference? PASS largest gap: sd rho, -1.4 SE
3. Particle count irrelevant Do J=25 and J=2000 give the same posterior? PASS largest gap: mean rho, +1.7 SE
4. Agrees with R Does pypomp match R pomp's pmcmc? PASS largest gap: sd beta1, -0.9 SE
5. Enough draws Does every run have at least 400 effective independent draws? PASS smallest bulk ESS 2,239 (R pomp, rho); split R-hat up to 1.033, see diagnostics

The grid uses pypomp’s own simulator and filter, so checks 2 and 3 test the sampler only. Model translation is covered by checks 1 and 4.

Visual review

1. Same model?

Particle-filter log-likelihood at the true parameters, pypomp vs R, on the same data (3600 filters × 2000 particles each at run level 4).

Show Code 
if precond is not None:
    display(HTML(ru.table_html(pd.DataFrame([precond]), fmt="%.3f")))
    both = pd.concat([
        py_pf[["logLik"]].assign(source="pypomp"),
        r_pf[["logLik"]].assign(source="R pomp"),
    ])
    display(
        ggplot(both, aes(x="logLik", fill="source", color="source"))
        + geom_density(alpha=0.4)
        + labs(title="Particle filter log-likelihood at the true parameters",
               subtitle="Each curve is the spread over replicate filters", x="log likelihood")
        + ru.scale_fill_premium()
        + ru.scale_color_premium()
        + ru.theme_premium
    )
else:
    display(Markdown("_Needs both run.py and run.R results._"))
pypomp R pomp difference SE z verdict
-1250.135 -1250.128 -0.007 0.003 -2.174 PASS

2. Right posterior?

pypomp draws vs the grid posterior. The grid is the particle-filter likelihood on a lattice over the prior box, normalized (60 × 60 points, 3 filters × 5000 particles per point at run level 4). For the grid to be adequate, it needs at least ~2 cells per posterior sd, and filter noise must be small compared with the log-density step between neighbours.

Show Code 
if have_grid:
    rows = []
    for param in PARAMS:
        m = grid_m[param]
        spacing = float(np.mean(np.diff(np.sort(grid[param].unique()))))
        rows.append({"parameter": param, "posterior sd": m["sd"],
                     "grid spacing": spacing, "cells per sd": m["sd"] / spacing})
    display(HTML(ru.table_html(pd.DataFrame(rows))))
    if grid["logLik_sd"].notna().any():
        near = grid[grid["log_post"] > grid["log_post"].max() - 10]
        step = np.abs(np.diff(near.sort_values(["beta1", "rho"])["log_post"])).mean()
        display(HTML(ru.table_html(pd.DataFrame([
            {"quantity": "median sd of one logLik replicate near the mode (nats)",
             "value": near["logLik_sd"].median()},
            {"quantity": "mean log-density step between neighbours near the mode", "value": step},
        ]), fmt="%.3f")))
parameter posterior sd grid spacing cells per sd
beta1 9.8854 3.3333 2.9656
rho 0.0194 0.0058 3.3214
quantity value
median sd of one logLik replicate near the mode (nats) 0.086
mean log-density step between neighbours near the mode 0.548
Show Code 
if have_py and have_grid:
    for param in PARAMS:
        axis, dens = ru.grid_marginals(grid, param, OTHER[param])
        frames = [ru.thin_for_plot(s["post"][[param]]).assign(source=f"pypomp J={J}")
                  for J, s in sorted(sweeps.items())]
        display(
            ggplot(pd.concat(frames), aes(x=param, color="source"))
            + geom_density(size=1)
            + geom_line(pd.DataFrame({param: axis, "density": dens}),
                        aes(x=param, y="density"), inherit_aes=False,
                        color=ru.color_palette["reference"], linetype="dashed", size=1)
            + geom_vline(xintercept=truth[param], color=ru.color_palette["truth"], linetype="dotted")
            + labs(title=f"PMCMC against the grid reference: {param}",
                   subtitle="Dashed black: grid reference. Dotted orange: true value.")
            + ru.theme_premium
        )
    for J in sorted(sweeps):
        display(Markdown(f"**pypomp J={J} vs grid**"))
        display(HTML(ru.table_html(vs_grid[J])))

pypomp J=25 vs grid

parameter statistic pypomp J=25 grid difference SE z verdict
beta1 mean 398.5806 398.5619 0.0187 0.1348 0.1388 PASS
beta1 sd 9.7732 9.8854 -0.1122 0.0953 -1.1777 PASS
rho mean 0.6307 0.6300 0.0007 0.0003 2.1584 PASS
rho sd 0.0193 0.0194 -0.0001 0.0002 -0.2517 PASS

pypomp J=2000 vs grid

parameter statistic pypomp J=2000 grid difference SE z verdict
beta1 mean 398.4822 398.5619 -0.0797 0.1027 -0.7758 PASS
beta1 sd 9.8912 9.8854 0.0058 0.0726 0.0795 PASS
rho mean 0.6301 0.6300 0.0000 0.0002 0.1702 PASS
rho sd 0.0191 0.0194 -0.0002 0.0002 -1.3913 PASS

3. Does the particle count matter?

Low-J and high-J chains must target the same posterior. The low-J arm is chosen so that the log-likelihood estimate is noisy (second table).

Show Code 
if vs_J is not None:
    display(HTML(ru.table_html(vs_J)))
if pf_vs_J is not None:
    noise = pf_vs_J.groupby("J")["logLik"].agg(["mean", "std"]).reset_index()
    noise.columns = ["J", "mean logLik at truth", "sd of logLik estimate"]
    display(Markdown("Log-likelihood estimate at the truth, by particle count:"))
    display(HTML(ru.table_html(noise, fmt="%.3f")))
parameter statistic J=25 J=2000 difference SE z verdict
beta1 mean 398.5806 398.4822 0.0984 0.1695 0.5807 PASS
beta1 sd 9.7732 9.8912 -0.1180 0.1198 -0.9848 PASS
rho mean 0.6307 0.6301 0.0007 0.0004 1.6670 PASS
rho sd 0.0193 0.0191 0.0002 0.0003 0.5950 PASS

Log-likelihood estimate at the truth, by particle count:

J mean logLik at truth sd of logLik estimate
10 -1252.381 2.486
25 -1250.776 1.513
100 -1250.330 0.584
2000 -1250.154 0.135
Show Code 
if have_py:
    rows = []
    for J, s in sorted(sweeps.items()):
        secs = (s["meta"].get("run_config") or {}).get("execution_time", np.nan)
        acc = s["acceptance"]
        row = {"source": f"pypomp J={J}",
               "acceptance": acc["acceptance_rate"].mean() if acc is not None else np.nan,
               "seconds": secs}
        for p in PARAMS:
            ess = ru.bulk_ess(ru.chains_array(s["post"], p))
            row[f"ESS {p}"] = ess
            row[f"ESS/s {p}"] = ess / secs if secs else np.nan
        rows.append(row)
    if have_r:
        secs = (r_meta.get("run_config") or {}).get("execution_time", np.nan)
        row = {"source": "R pomp", "acceptance": r_acc["acceptance_rate"].mean()
               if r_acc is not None else np.nan, "seconds": secs}
        for p in PARAMS:
            ess = ru.bulk_ess(ru.chains_array(r_post, p))
            row[f"ESS {p}"] = ess
            row[f"ESS/s {p}"] = ess / secs if secs else np.nan
        rows.append(row)
    display(Markdown("**Cost and efficiency.** ESS counts draws kept after burn-in."))
    display(HTML(ru.table_html(pd.DataFrame(rows), fmt="%.3g")))

Cost and efficiency. ESS counts draws kept after burn-in.

source acceptance seconds ESS beta1 ESS/s beta1 ESS rho ESS/s rho
pypomp J=25 0.307 336 5.26e+03 15.6 3.28e+03 9.75
pypomp J=2000 0.652 809 9.27e+03 11.5 6.5e+03 8.03
R pomp 0.521 1.09e+03 3.42e+03 3.13 2.24e+03 2.05

4. Agrees with R?

pypomp at the highest J vs R pomp::pmcmc. pomp’s random walk is on the natural scale, and pypomp’s is on log (beta1) and logit (rho) scales, so acceptance rates are not comparable. The posteriors should match.

Show Code 
if vs_R is not None:
    for param in PARAMS:
        both = pd.concat([
            ru.thin_for_plot(sweeps[J_HI]["post"][[param]]).assign(source="pypomp"),
            ru.thin_for_plot(r_post[[param]]).assign(source="R pomp"),
        ])
        display(
            ggplot(both, aes(x=param, fill="source", color="source"))
            + geom_density(alpha=0.35)
            + geom_vline(xintercept=truth[param], color=ru.color_palette["truth"], linetype="dashed")
            + labs(title=f"pypomp against R pomp: {param}",
                   subtitle=f"pypomp at J={J_HI}; dashed line is the true value")
            + ru.scale_fill_premium()
            + ru.scale_color_premium()
            + ru.theme_premium
        )
    display(HTML(ru.table_html(vs_R)))
else:
    display(Markdown("_Needs both run.py and run.R results._"))

parameter statistic pypomp J=2000 R pomp difference SE z verdict
beta1 mean 398.4822 398.5222 -0.0400 0.1998 -0.2001 PASS
beta1 sd 9.8912 10.0211 -0.1299 0.1413 -0.9195 PASS
rho mean 0.6301 0.6300 0.0001 0.0005 0.2190 PASS
rho sd 0.0191 0.0192 -0.0001 0.0003 -0.3152 PASS

5. Are there enough draws?

Every run must have bulk ESS ≥ 400 for its standard errors to be trusted. Split R-hat is for information only.

Show Code 
if len(conv):
    display(HTML(ru.table_html(conv, fmt="%.3f")))
source parameter chains draws per chain split R-hat bulk ESS
pypomp J=25 beta1 128 100 1.017 5259.145
pypomp J=25 rho 128 100 1.033 3280.236
pypomp J=2000 beta1 128 100 1.005 9268.525
pypomp J=2000 rho 128 100 1.008 6501.527
R pomp beta1 36 1500 1.010 3419.139
R pomp rho 36 1500 1.013 2239.134
Show Code 
if have_py:
    for J, s in sorted(sweeps.items()):
        tr = s["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 J={J}: {param}",
                       subtitle="Left of the dotted line is discarded as burn-in",
                       x="iteration", color="chain")
                + ru.theme_premium
            )

6. Do the credible intervals cover the truth?

Sanity check, not scored: one dataset gives a single draw of coverage.

Show Code 
if have_py:
    rows = []
    sources = [(f"pypomp J={J_HI}", sweeps[J_HI]["post"])] + ([("R pomp", r_post)] if have_r else [])
    for label, post in sources:
        for param in PARAMS:
            lo, hi = np.quantile(post[param], [0.05, 0.95])
            rows.append({"source": label, "parameter": param, "truth": truth[param],
                         "5%": lo, "95%": hi, "covers truth": bool(lo <= truth[param] <= hi)})
    display(HTML(ru.table_html(pd.DataFrame(rows))))
source parameter truth 5% 95% covers truth
pypomp J=2000 beta1 400.0000 382.0563 414.6694 True
pypomp J=2000 rho 0.6000 0.5992 0.6624 True
R pomp beta1 400.0000 382.1794 415.1275 True
R pomp rho 0.6000 0.5991 0.6625 True

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("**Grid reference run**\n\n" + ru.format_metadata(ref_meta)))

pypomp run

  • Run at: 2026-10-08T00:43:20+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=pmcmc, model=sir, RUN_LEVEL=4, USE_CPU=False, MAIN_SEED=3141593, NCHAINS=128, M=2000, J=2000, free_params=[‘beta1’, ‘rho’], rw_sd_estimation_scale={‘beta1’: 0.02, ‘rho’: 0.05}, prior_box={‘beta1’: [300.0, 500.0], ‘rho’: [0.45, 0.8]}, execution_time=809.3712360858917, platform=gpu, trace_thin=10

R pomp run

  • Run at: 2026-10-06T12:06:18
  • pomp: 6.3 / R: 4.4.0
  • Commit: 7221fdf231
  • CPU: Intel(R) Xeon(R) Gold 6154 CPU @ 3.00GHz (36 cores)
  • Config: kind=pmcmc, model=sir, RUN_LEVEL=4, NCHAINS=36, Nmcmc=3000, Np=100, free_params=[‘beta1’, ‘rho’], rw_sd_natural_scale={‘beta1’: 8, ‘rho’: 0.012}, execution_time=1092.382

Grid 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