Dhaka Cholera: Fixed-Parameter Likelihood Validation

Particle Filter Distribution Validation: R pomp vs. pypomp

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

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

IF2 vs IFAD  |  Timing & Throughput  |  Likelihood Evaluation

Introduction

To check that pypomp’s particle filter is doing the same arithmetic as R pomp’s, both are run many times at one fixed parameter vector — the published MLE that pomp ships as dacca()’s default — and their log-likelihood distributions are compared. Nothing varies across replicates except the seed, so any systematic gap between the two curves is an implementation difference rather than Monte Carlo noise.

---

Benchmark Settings & Environment

The table below summarizes algorithmic parameters, software environments, and compute hardware recorded in latest.json for each configuration.

PLATFORMS = {
    "R pomp": os.path.join("results", "R"),
    "pypomp (GPU)": os.path.join("results", "gpu"),
}

runs = ru.load_timing_data(PLATFORMS)
display(HTML(ru.build_settings_comparison_html(runs, is_panel=False)))

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

missing = [(label, r["dir"]) for label, r in runs.items() if not r["available"]]
for label, path in missing:
    display(HTML(f"<div class='alert alert-warning'><strong>Missing results for {label}:</strong> Expected at <code>{path}</code>.</div>"))

import numpy as np
import pandas as pd
from scipy.stats import ks_2samp
from plotnine import ggplot, aes, geom_density, labs

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

py_csv = os.path.join("results", "gpu", "pfilter_logliks.csv")
r_csv = os.path.join("results", "R", "pfilter_logliks.csv")

PY = pd.read_csv(py_csv)["logLik"].values if os.path.exists(py_csv) else np.array([])
R = pd.read_csv(r_csv)["logLik"].dropna().values if os.path.exists(r_csv) else np.array([])

---

Likelihood Evaluation at the MLE

Kolmogorov-Smirnov (KS) test

ks_stat, ks_p = ks_2samp(PY, R)

ks_table = pd.DataFrame({
    "Source": ["pypomp (GPU)", "R pomp"],
    "Replicates": [len(PY), len(R)],
    "Mean": [np.mean(PY), np.mean(R)],
    "SD": [np.std(PY, ddof=1), np.std(R, ddof=1)],
    "Min": [np.min(PY), np.min(R)],
    "Median": [np.median(PY), np.median(R)],
    "Max": [np.max(PY), np.max(R)],
})

display(HTML(ks_table.round(3).to_html(classes="table table-striped table-hover", index=False)))
print(
    f"Mean difference: {np.mean(PY) - np.mean(R):+.4f} nats\n"
    f"Two-sample KS:   D = {ks_stat:.4f}, p = {ks_p:.4g}"
)

Source	Replicates	Mean	SD	Min	Median	Max
pypomp (GPU)	3600	-3748.43	0.799	-3751.178	-3748.436	-3745.343
R pomp	3600	-3748.43	0.807	-3751.499	-3748.424	-3745.534

Mean difference: +0.0003 nats
Two-sample KS:   D = 0.0125, p = 0.9412

---

Density comparison

df = pd.concat([
    pd.DataFrame({"logLik": PY, "source": "python"}),
    pd.DataFrame({"logLik": R, "source": "R"}),
], ignore_index=True)

(
    ggplot(df, aes(x="logLik", fill="source", color="source"))
    + geom_density(alpha=0.3)
    + labs(
        title="Log-likelihood density at the MLE: pypomp vs R pomp",
        subtitle="Monte Carlo variation across identical independent evaluations",
        x="Log-Likelihood", y="Density", fill="Source", color="Source",
    )
    + ru.scale_fill_premium()
    + ru.scale_color_premium()
    + ru.theme_premium
)

Expectation

The two densities should overlap heavily. A KS test that fails to reject is the positive result here; a visible shift in location would point at a difference in the transition process or the measurement model rather than at sampling noise.