UK Measles Model Comparison

Log-likelihood and parameter estimates in pypomp vs. R pomp

Introduction

This report compares the distributions of the log-likelihood and parameter estimates obtained for the measles model using the particle filter and iterated filtering 2 (IF2) in pypomp and pomp.

---

Setup

import os
import logging
import pickle
import numpy as np
import pandas as pd
import pyreadr
from scipy.special import logit, logsumexp
from plotnine import *
from IPython.display import HTML

# Hide the noisy JAX CUDA initialization error on machines without a GPU
logging.getLogger("jax._src.xla_bridge").setLevel(logging.CRITICAL)

# ---------------------------------------------------------
# Load Parameter Comparison Data
# ---------------------------------------------------------
# Load Python parameter results
py_param_path = "parameter_comparison/results/mif_coefs.pkl"
with open(py_param_path, "rb") as f:
    py_results = pickle.load(f)

py_results["source"] = "pypomp"
py_results = py_results.drop(
    columns=["logLik", "se", "theta_idx"], errors="ignore"
).melt(id_vars=["source", "unit"], var_name="param", value_name="value")

# Load R parameter results
r_param_path = "parameter_comparison/results/mif_coefs.rds"
r_data = pyreadr.read_r(r_param_path)
r_results = r_data[None]
r_results["source"] = "R"
r_results = r_results[["source", "unit", "coef", "names"]].rename(
    columns={"coef": "value", "names": "param"}
)

# ---------------------------------------------------------
# Load Particle Filter Likelihood Data
# ---------------------------------------------------------
py_pf_32_path = "logLik_comparison/results/pfilter_logliks_f32.pkl"
py_pf_64_path = "logLik_comparison/results/pfilter_logliks_f64.pkl"
r_pf_path = "logLik_comparison/results/pfilter_logliks_f64.rds"

with open(py_pf_32_path, "rb") as f:
    pf_py_results = pickle.load(f)
pf_py_results["source"] = "pypomp (32-bit)"

with open(py_pf_64_path, "rb") as f:
    pf_py_exact_results = pickle.load(f)
pf_py_exact_results["source"] = "pypomp (64-bit)"

pf_r_results = pyreadr.read_r(r_pf_path)[None]
pf_r_results["source"] = "R"

---

Parameter Convergence & Distributions

This section compares the final parameter values obtained from the IF2 algorithm across replicates. Ideally, these should overlap as much as possible.

Density Plots

We overlay the densities of 360 final parameter estimates for each unit to assess the agreement between pypomp and R.

---

Particle Filter Likelihood Comparison

In this section, we compare 3600 likelihood evaluations from the particle filters across different configurations. The configurations are as follows:

1. pypomp run using 32-bit floating point precision.
2. pypomp run using 64-bit floating point precision.
3. pomp run using 64-bit floating point precision.

We include pypomp with 64-bit floating point precision due to some previous results suggesting precision makes a difference in this model. For any given unit, we use the same parameter set for each evaluation.

Summary Statistics

The table below presents the summary statistics of the log-likelihood (LL) calculations.

unit	source	mean_logLik	sd_logLik	min_logLik	max_logLik	logmeanexp_logLik
Halesworth	R	-319.349	1.013	-323.392	-315.126	-318.865
Halesworth	pypomp (32-bit)	-319.435	1.024	-323.886	-316.213	-318.962
Halesworth	pypomp (64-bit)	-319.380	1.029	-323.703	-316.189	-318.893
London	R	-3831.162	1.276	-3836.289	-3827.186	-3830.392
London	pypomp (32-bit)	-3831.301	1.278	-3836.055	-3827.163	-3830.541
London	pypomp (64-bit)	-3831.162	1.281	-3836.358	-3827.332	-3830.390

Log-Likelihood Distributions

Below are the density plots of the calculated log-likelihoods for each unit (London and Halesworth). Ideally, these overlap as much as possible.