Measles Fixed-Parameter Likelihood Validation

Particle filter log-likelihood distribution: 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>"))

Likelihood Evaluation  |  Parameter Estimation  |  Timing & Throughput  |  IF2 vs IFAD

Introduction

This report evaluates the particle filter’s likelihood calculation itself, by running it many times at the He et al. (2010) maximum likelihood estimates for each town. Nothing varies across replicates except the random seed, so any systematic difference between pypomp and pomp is a difference between the implementations.

pypomp is run twice, in single and in double precision. We have seen some indication that single precision slightly shifts the estimated log-likelihood for the measles model so we check that here. pomp computes in double precision throughout, so the 64-bit arm is the like-for-like comparison and the 32-bit arm measures what the default costs.

Two towns are used: London, the largest unit, and Halesworth, among the smallest.

---

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 (32-bit)": os.path.join("results", "f32"),
    "pypomp (64-bit)": os.path.join("results", "f64"),
}

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

Setting / Parameter	R pomp	pypomp (32-bit)	pypomp (64-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)	32-bit (float32)	64-bit (float64)
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.0rc1	pypomp 1.0.0rc1
Backend / Engine	R 4.4.0	JAX 0.11.1	JAX 0.11.1
Quant Git Commit	7f9ea97	630cc23	630cc23
Run Timestamp	2026-08-10 21:51:33	2026-09-02 19:52:20	2026-09-02 21:04:40
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	59687793	59687795

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, facet_wrap, labs

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

frames = []
for label, d in PLATFORMS.items():
    df = ru.read_if_exists(os.path.join(d, "pfilter_logliks.csv"))
    if df is not None:
        df = df[["unit", "replicate", "logLik"]].copy()
        df["source"] = label
        frames.append(df)

combined = pd.concat(frames, ignore_index=True) if frames else pd.DataFrame()

---

Summary Statistics

logmeanexp is the log of the mean likelihood and sits above the mean log-likelihood by roughly sd^2/2. It is the quantity pomp reports as the likelihood estimate, so it is listed alongside the raw mean.

summary = (
    combined.groupby(["unit", "source"])["logLik"]
    .agg(
        n="count",
        mean="mean",
        sd="std",
        min="min",
        max="max",
        logmeanexp=ru.logmeanexp,
    )
    .reset_index()
    .sort_values(["unit", "source"])
)

display(HTML(summary.round(3).to_html(classes="table table-striped table-hover", index=False)))

unit	source	n	mean	sd	min	max	logmeanexp
Halesworth	R pomp	3600	-319.387	1.023	-324.472	-316.214	-318.913
Halesworth	pypomp (32-bit)	3600	-319.422	1.002	-323.313	-316.020	-318.960
Halesworth	pypomp (64-bit)	3600	-319.379	1.026	-323.924	-316.504	-318.903
London	R pomp	3600	-3831.169	1.284	-3835.509	-3827.471	-3830.400
London	pypomp (32-bit)	3600	-3831.327	1.273	-3836.108	-3827.097	-3830.558
London	pypomp (64-bit)	3600	-3831.119	1.257	-3835.992	-3826.338	-3830.346

Kolmogorov-Smirnov tests against the R baseline

rows = []
if not combined.empty and "R" in set(combined["source"]):
    for unit in sorted(combined["unit"].unique()):
        r_ll = combined[(combined["unit"] == unit) & (combined["source"] == "R")]["logLik"]
        for label in ["pypomp (32-bit)", "pypomp (64-bit)"]:
            py_ll = combined[(combined["unit"] == unit) & (combined["source"] == label)]["logLik"]
            if len(py_ll) == 0 or len(r_ll) == 0:
                continue
            d, p = ks_2samp(py_ll, r_ll)
            rows.append({
                "Unit": unit,
                "Comparison": f"{label} vs R",
                "Mean difference (nats)": py_ll.mean() - r_ll.mean(),
                "KS D": d,
                "KS p": p,
            })

ks_table = pd.DataFrame(rows)
display(HTML(ks_table.round(4).to_html(classes="table table-striped table-hover", index=False)))

---

Density Comparison

(
    ggplot(combined, aes(x="logLik", fill="source", color="source"))
    + geom_density(alpha=0.35)
    + facet_wrap("~unit", scales="free")
    + labs(
        title="Particle filter log-likelihood at the He et al. (2010) estimates",
        subtitle="One density per implementation, per town",
        x="Log-likelihood",
        y="Density",
        fill="Source",
        color="Source",
    )
    + ru.scale_fill_premium()
    + ru.scale_color_premium()
    + ru.theme_premium
)

Expectation

The three densities should overlap heavily within each town. A shift in the 32-bit curve relative to the other two localises the difference to floating point precision rather than to the filter.