Fixed-Parameter Likelihood Validation

Particle filter log-likelihood distribution: R pomp vs. pypomp, all models

import os

import numpy as np
import pandas as pd
import report_utils as ru  # puts tests/ on sys.path
import checks
from IPython.display import HTML, display
from plotnine import aes, facet_wrap, geom_density, geom_vline, ggplot, labs

ru.nav_bar("loglik")

Likelihood Evaluation  |  Timing & Throughput  |  Scaling  |  All reports

Introduction

This report checks the particle filter’s likelihood estimate by running it many times at a single fixed parameter vector for each model. Only the random seed varies across replicates, so any systematic difference between pypomp and R (pomp, or panelPomp for the panel) would be a difference between the implementations.

• SPX: the estimates from Sun (2024).
• Dhaka cholera: the published maximum likelihood estimate (MLE) that pomp uses as the default for dacca().
• Measles: the He et al. (2010) MLEs for London, the largest unit, and Halesworth, one of the smallest. pypomp runs once in single and once in double precision. Since pomp uses double precision throughout, the 64-bit run is the direct comparison, while the 32-bit run shows the effect of pypomp’s default precision.
• Panel measles: the He et al. (2010) MLEs for London, Halesworth, Hastings and Cardiff. The parameters R0, sigma, gamma, sigmaSE, cohort and amplitude are shared across units and set to the average of the four towns’ estimates, while iota, rho, psi, S_0, E_0, I_0 and R_0 are unit-specific. Since the panel log-likelihood is the sum of the unit log-likelihoods, we compare each unit separately.

Why the panel covariates are imported from pomp

The Python and R smoothing splines interpolate the birth rate and population covariates differently, most noticeably for Cardiff, and this alone shifts the log-likelihood. Since we want to compare the filters rather than two slightly different models, we export the interpolated covariates from pomp to tests/panel_measles/R_covariates.csv and load them into the pypomp objects (model.align_covariates). As such, both packages filter the same model.

MODELS = ru.models("loglik")
runs = {m: ru.load_runs(m, "loglik") for m in MODELS}


def load_logliks(model):
    spec = ru.MODELS[model]["loglik"]
    frames = []
    for arm, sub in spec["arms"].items():
        path = os.path.join(ru.results_dir(spec["dir"], sub), "pfilter_logliks.csv")
        if not os.path.exists(path):
            continue
        df = pd.read_csv(path).dropna(subset=["logLik"])
        if "unit" not in df.columns:
            df["unit"] = ""
        if arm == spec["reference"]:
            role = "R"
        elif arm in spec.get("info", []):
            role = arm
        else:
            role = "pypomp"
        frames.append(df[["unit", "logLik"]].assign(model=model, source=arm, role=role))
    return pd.concat(frames, ignore_index=True) if frames else pd.DataFrame(
        columns=["unit", "logLik", "model", "source", "role"]
    )


combined = pd.concat([load_logliks(m) for m in MODELS], ignore_index=True)
combined["label"] = combined["model"].map(lambda m: ru.MODELS[m]["label"])
combined["panel"] = np.where(
    combined["unit"] == "", combined["label"], combined["label"] + "\n" + combined["unit"]
)
combined["panel"] = pd.Categorical(combined["panel"], categories=combined["panel"].unique())

ru.show_missing(runs)

---

Automated checks

records = {}
print(f"*{checks.CAVEAT}*\n")
for model in MODELS:
    spec = ru.MODELS[model]["loglik"]
    ref, info = spec["reference"], spec.get("info", [])
    df = combined[combined["model"] == model]
    units = sorted(df["unit"].unique()) or [""]
    arms = [a for a in spec["arms"] if a != ref]
    m = 2 * len(units) * sum(a not in info for a in arms)
    rows = []
    for arm in arms:
        for unit in units:
            sel = df["unit"] == unit
            compared = f"{arm} vs {ref}" if unit == "" else f"{unit}: {arm} vs {ref}"
            rows += checks.loglik_checks(
                df[sel & (df["source"] == arm)]["logLik"],
                df[sel & (df["source"] == ref)]["logLik"],
                compared,
                m,
                verdict=arm not in info,
            )
    level = checks.run_level(runs[model])
    rows = checks.gate_smoke(rows, level)
    version = runs[model][arms[0]]["meta"].get("pypomp_version")
    records[model] = checks.record(rows, level, pypomp_version=version)

    print(ru.heading(model))
    print(f"Verdict: {checks.badge(records[model]['verdict'])}\n")
    if info:
        print(f"Rows for {', '.join(info)} are for reference and receive no verdict.\n")
    print(ru.raw_html(checks.rows_html(rows)))

print(f"\n<small>{checks.rule_text()}</small>\n")
checks.save_models(records, "loglik_checks.json")

These checks are computed from the results each time the report is rendered. They only compare summary statistics, so the plots can show problems that the checks miss. Visual review lists what the plots should show.

SPX

Verdict: PASS

Check	Compared	Result	Verdict
Mean log-likelihood	pypomp (GPU) vs R pomp	difference +0.050 nats (SE 0.049), z = +1.04	PASS
Distribution (KS)	pypomp (GPU) vs R pomp	D = 0.0275, p = 0.131	PASS

Dhaka cholera

Verdict: PASS

Check	Compared	Result	Verdict
Mean log-likelihood	pypomp (GPU) vs R pomp	difference +0.000 nats (SE 0.019), z = +0.02	PASS
Distribution (KS)	pypomp (GPU) vs R pomp	D = 0.0125, p = 0.941	PASS

Measles

Verdict: PASS

Rows for pypomp (32-bit) are for reference and receive no verdict.

Check	Compared	Result	Verdict
Mean log-likelihood	Halesworth: pypomp (64-bit) vs R pomp	difference -0.004 nats (SE 0.024), z = -0.15	PASS
Distribution (KS)	Halesworth: pypomp (64-bit) vs R pomp	D = 0.0161, p = 0.739	PASS
Mean log-likelihood	London: pypomp (64-bit) vs R pomp	difference -0.022 nats (SE 0.030), z = -0.74	PASS
Distribution (KS)	London: pypomp (64-bit) vs R pomp	D = 0.0219, p = 0.351	PASS
Mean log-likelihood	Halesworth: pypomp (32-bit) vs R pomp	difference -0.009 nats (SE 0.024), z = -0.39	info
Distribution (KS)	Halesworth: pypomp (32-bit) vs R pomp	D = 0.0181, p = 0.6	info
Mean log-likelihood	London: pypomp (32-bit) vs R pomp	difference -0.143 nats (SE 0.030), z = -4.68	info
Distribution (KS)	London: pypomp (32-bit) vs R pomp	D = 0.0567, p = 1.9e-05	info

Panel measles

Verdict: PASS

Check	Compared	Result	Verdict
Mean log-likelihood	Cardiff: pypomp (GPU) vs R panelPomp	difference +0.084 nats (SE 0.053), z = +1.58	PASS
Distribution (KS)	Cardiff: pypomp (GPU) vs R panelPomp	D = 0.0253, p = 0.2	PASS
Mean log-likelihood	Halesworth: pypomp (GPU) vs R panelPomp	difference +0.026 nats (SE 0.032), z = +0.83	PASS
Distribution (KS)	Halesworth: pypomp (GPU) vs R panelPomp	D = 0.0192, p = 0.523	PASS
Mean log-likelihood	Hastings: pypomp (GPU) vs R panelPomp	difference +0.025 nats (SE 0.043), z = +0.57	PASS
Distribution (KS)	Hastings: pypomp (GPU) vs R panelPomp	D = 0.0203, p = 0.45	PASS
Mean log-likelihood	London: pypomp (GPU) vs R panelPomp	difference -0.097 nats (SE 0.069), z = -1.40	PASS
Distribution (KS)	London: pypomp (GPU) vs R panelPomp	D = 0.0186, p = 0.561	PASS

PASS when |z| < 3, CHECK when 3 ≤ |z| < 5, and FAIL when |z| ≥ 5. KS p-values use the equivalent thresholds. Both are Bonferroni-adjusted within each model over the tests that receive a verdict.

Visual review

• Density comparison: in every panel, the pypomp and R densities should overlap heavily. For measles, the 32-bit and 64-bit runs differ only in precision, so any shift of the 32-bit density away from the other two should be due to precision rather than the filter.

---

Density comparison

Each panel is centered on the mean of the R replicates, so the horizontal axis is in nats from that mean.

r_mean = (
    combined[combined["role"] == "R"]
    .groupby("panel", observed=True)["logLik"]
    .mean()
    .rename("r_mean")
)
plot_df = combined.join(r_mean, on="panel")
plot_df["centered"] = plot_df["logLik"] - plot_df["r_mean"]

(
    ggplot(plot_df, aes(x="centered", fill="role", color="role"))
    + geom_density(alpha=0.3)
    + geom_vline(xintercept=0, linetype="dashed", color="#7f8c8d")
    + facet_wrap("~panel", scales="free", ncol=4)
    + labs(
        title="Particle filter log-likelihood at fixed parameters",
        subtitle="One density per implementation, per model and unit",
        x="Log-likelihood minus the R mean (nats)",
        y="Density",
        fill="Source",
        color="Source",
    )
    + ru.scale_fill_premium()
    + ru.scale_color_premium()
    + ru.theme_premium
)

---

Summary statistics

We include logmeanexp, the log of the mean likelihood, since it is the likelihood estimate that pomp reports.

summary = (
    combined.groupby(["model", "label", "unit", "source"], sort=False)["logLik"]
    .agg(n="count", mean="mean", sd="std", min="min", max="max", logmeanexp=ru.logmeanexp)
    .reset_index()
    .drop(columns="model")
    .rename(columns={"label": "model"})
)
display(HTML(summary.round(3).to_html(classes="table table-striped table-hover", index=False)))

model	unit	source	n	mean	sd	min	max	logmeanexp
SPX		R pomp	3600	11848.080	2.047	11841.868	11856.854	11850.702
SPX		pypomp (GPU)	3600	11848.131	2.091	11841.511	11855.983	11850.594
Dhaka cholera		R pomp	3600	-3748.430	0.807	-3751.499	-3745.534	-3748.108
Dhaka cholera		pypomp (GPU)	3600	-3748.430	0.799	-3751.178	-3745.343	-3748.110
Measles	London	R pomp	3600	-3831.169	1.284	-3835.509	-3827.471	-3830.400
Measles	Halesworth	R pomp	3600	-319.387	1.023	-324.472	-316.214	-318.913
Measles	London	pypomp (64-bit)	3600	-3831.191	1.283	-3836.151	-3826.560	-3830.398
Measles	Halesworth	pypomp (64-bit)	3600	-319.391	1.032	-324.317	-316.086	-318.899
Measles	London	pypomp (32-bit)	3600	-3831.311	1.302	-3836.896	-3826.956	-3830.500
Measles	Halesworth	pypomp (32-bit)	3600	-319.397	1.013	-323.847	-316.276	-318.921
Panel measles	London	R panelPomp	3600	-3953.893	2.950	-3966.563	-3944.690	-3950.535
Panel measles	Halesworth	R panelPomp	3600	-323.290	1.342	-328.727	-319.271	-322.526
Panel measles	Hastings	R panelPomp	3600	-1592.498	1.855	-1599.050	-1586.051	-1590.693
Panel measles	Cardiff	R panelPomp	3600	-2381.299	2.219	-2391.755	-2375.465	-2379.522
Panel measles	London	pypomp (GPU)	3600	-3953.990	2.906	-3969.566	-3945.263	-3950.864
Panel measles	Halesworth	pypomp (GPU)	3600	-323.264	1.345	-332.314	-319.286	-322.487
Panel measles	Hastings	pypomp (GPU)	3600	-1592.473	1.786	-1598.612	-1585.347	-1590.684
Panel measles	Cardiff	pypomp (GPU)	3600	-2381.216	2.272	-2391.078	-2375.206	-2379.365

Panel log-likelihood

The sum of the unit logmeanexp values is the panel log-likelihood estimate that each package would report.

panel = (
    summary[summary["model"] == ru.MODELS["panel_measles"]["label"]]
    .groupby("source")["logmeanexp"]
    .sum()
    .reset_index()
    .rename(columns={"logmeanexp": "panel logLik"})
)
display(HTML(panel.round(3).to_html(classes="table table-striped table-hover", index=False)))

source	panel logLik
R panelPomp	-8243.275
pypomp (GPU)	-8243.400

---

Appendix: settings and environment

print(ru.settings_sections(runs))

SPX

Setting / Parameter	R pomp	pypomp (GPU)
Algorithmic & Workload Settings
Run Level	4	4
Pfilter Particles ($N_{p,eval}$)	1,000	1,000
Evaluation Replicates ($N_{reps}$)	3,600	3,600
Random Seed	—	631409
Software & Environment
Pomp Framework	pomp 6.3	pypomp 1.0.0rc1
Backend / Engine	R 4.4.0	JAX 0.11.1
Quant Git Commit	8853fc5	630cc23
Run Timestamp	2026-08-05 16:32:31	2026-09-02 18:59:37
Hardware & Compute
Compute Device	Intel(R) Xeon(R) Gold 6254 CPU @ 3.10GHz (36 cores)	Tesla V100-PCIE-16GB (1 GPU)
Slurm Partition	standard	gpu
Slurm Job ID	56537488	59673767

Dhaka cholera

Setting / Parameter	R pomp	pypomp (GPU)
Algorithmic & Workload Settings
Run Level	4	4
Pfilter Particles ($N_{p,eval}$)	5,000	5,000
Evaluation Replicates ($N_{reps}$)	3,600	3,600
Random Seed	—	631409
Software & Environment
Pomp Framework	pomp 6.3	pypomp 1.0.0rc1
Backend / Engine	R 4.4.0	JAX 0.11.1
Quant Git Commit	3646a0c	630cc23
Run Timestamp	2026-08-06 12:38:32	2026-09-02 18:48:10
Hardware & Compute
Compute Device	Intel(R) Xeon(R) Gold 6154 CPU @ 3.00GHz (36 cores)	Tesla V100-PCIE-16GB (1 GPU)
Slurm Partition	standard	gpu
Slurm Job ID	56607557	59673748

Measles

Setting / Parameter	R pomp	pypomp (64-bit)	pypomp (32-bit)
Algorithmic & Workload Settings
Run Level	4	4	4
Pfilter Particles ($N_{p,eval}$)	5,000	5,000	5,000
Evaluation Replicates ($N_{reps}$)	3,600	3,600	3,600
Floating-Point Precision	64-bit (double)	64-bit (float64)	32-bit (float32)
Unit(s)	London, Halesworth (2 units)	London, Halesworth (2 units)	London, Halesworth (2 units)
Random Seed	—	594709947	594709947
Software & Environment
Pomp Framework	pomp 6.3	pypomp 1.0.6	pypomp 1.0.6
Backend / Engine	R 4.4.0	JAX 0.11.2	JAX 0.11.2
Quant Git Commit	7f9ea97	9c649ea	9c649ea
Run Timestamp	2026-08-10 21:51:33	2026-09-25 15:54:38	2026-09-24 20:35:50
Hardware & Compute
Compute Device	Intel(R) Xeon(R) Gold 6154 CPU @ 3.00GHz (36 cores)	NVIDIA RTX PRO 6000 Blackwell Server Edition (1 GPU)	NVIDIA RTX PRO 6000 Blackwell Server Edition (1 GPU)
Slurm Partition	standard	gpu-rtx6000	gpu-rtx6000
Slurm Job ID	56880073	61873140	61826878

Panel measles

Setting / Parameter	R panelPomp	pypomp (GPU)
Algorithmic & Workload Settings
Run Level	4	4
Pfilter Particles / Unit ($N_{p,eval}$)	5,000	5,000
Evaluation Replicates ($N_{reps}$)	3,600	3,600
Unit(s)	London, Halesworth, Hastings, Cardiff (4 units)	London, Halesworth, Hastings, Cardiff (4 units)
Shared Parameters	—	R0, sigma, gamma, sigmaSE, cohort, amplitude (6 params)
Random Seed	—	594709947
Software & Environment
Pomp Framework	panelPomp 1.7.0.0 (pomp 6.3)	pypomp 1.0.6
Backend / Engine	R 4.4.0	JAX 0.11.2
Quant Git Commit	b22249b	9c649ea
Run Timestamp	2026-08-12 13:00:06	2026-09-24 22:13:13
Hardware & Compute
Compute Device	Intel(R) Xeon(R) Gold 6154 CPU @ 3.00GHz (36 cores)	NVIDIA RTX PRO 6000 Blackwell Server Edition (1 GPU)
Slurm Partition	standard	gpu-rtx6000
Slurm Job ID	57205142	61839613

References

Sun, Weizhe. 2024. “Model Based Inference of Stochastic Volatility via Iterated Filtering.” Undergraduate Honors Thesis, University of Michigan.