SPX Comparison

Performance, Convergence, and Likelihood Validation: R pomp vs. pypomp

Published

June 24, 2026

Introduction

This report compares the performance, parameter convergence, and likelihood evaluation of the S&P 500 (SPX) stochastic volatility model used in Sun (2024) using the pypomp package (on CPU and GPU) and the R pomp package (on CPU).

The SPX model is a stochastic volatility model with the following parameter set:

$\mu$: Expected rate of return of the asset.
$\kappa$: Mean rate of reversion in variance.
$\theta$: Long term average variance.
$\xi$: Volatility coefficient of volatility.
$\rho$: Leverage effect parameter (correlation between return and volatility shocks).
$V_0$: Initial variance.

Because the SPX model uses simple random number generation (sampling from a normal distribution) and features a single rproc step per observation, it is highly sensitive to framework overhead in the IF2 optimization and particle filter (pfilter) algorithms. This benchmark compares runtimes, convergence trajectories, and likelihood distributions to verify the mathematical correctness and efficiency of pypomp.

Setup

System & Session Information (R)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.2     ✔ tibble    3.3.0
✔ lubridate 1.9.4     ✔ tidyr     1.3.1
✔ purrr     1.1.0     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

Attaching package: 'pomp'


The following object is masked from 'package:purrr':

    map

R version 4.4.0 (2024-04-24)
Platform: x86_64-pc-linux-gnu
Running under: Red Hat Enterprise Linux 8.10 (Ootpa)

Matrix products: default
BLAS:   /sw/pkgs/arc/stacks/gcc/13.2.0/R/4.4.0/lib64/R/lib/libRblas.so 
LAPACK: /sw/pkgs/arc/stacks/gcc/13.2.0/R/4.4.0/lib64/R/lib/libRlapack.so;  LAPACK version 3.12.0

locale:
 [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8       
 [4] LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8   
 [7] LC_PAPER=C.UTF-8       LC_NAME=C              LC_ADDRESS=C          
[10] LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   

time zone: America/Detroit
tzcode source: system (glibc)

attached base packages:
[1] stats     graphics  grDevices datasets  utils     methods   base     

other attached packages:
 [1] pomp_6.3          lubridate_1.9.4   forcats_1.0.0     stringr_1.5.1    
 [5] dplyr_1.1.4       purrr_1.1.0       readr_2.1.5       tidyr_1.3.1      
 [9] tibble_3.3.0      ggplot2_3.5.2     tidyverse_2.0.0   reticulate_1.42.0

loaded via a namespace (and not attached):
 [1] Matrix_1.7-0       gtable_0.3.6       jsonlite_2.0.0     compiler_4.4.0    
 [5] renv_1.0.11        tidyselect_1.2.1   Rcpp_1.1.0         scales_1.4.0      
 [9] png_0.1-8          yaml_2.3.10        fastmap_1.2.0      here_1.0.1        
[13] coda_0.19-4.1      lattice_0.22-6     R6_2.6.1           generics_0.1.4    
[17] knitr_1.50         rprojroot_2.0.4    tzdb_0.5.0         pillar_1.11.0     
[21] RColorBrewer_1.1-3 rlang_1.1.6        deSolve_1.40       stringi_1.8.7     
[25] xfun_0.52          timechange_0.3.0   cli_3.6.5          withr_3.0.2       
[29] magrittr_2.0.3     digest_0.6.37      grid_4.4.0         mvtnorm_1.3-3     
[33] hms_1.1.3          lifecycle_1.0.4    vctrs_0.6.5        data.table_1.17.8 
[37] evaluate_1.0.4     glue_1.8.0         farver_2.1.2       rmarkdown_2.29    
[41] tools_4.4.0        pkgconfig_2.0.3    htmltools_0.5.8.1

[1] "2026-06-24 12:10:11 EDT"

Metadata & Execution Runtimes

System Configurations

GPU Pomp Object Metadata:

Model Initialization Metadata:
------------------------------
pypomp version: 0.4.6.0
JAX version:    0.9.0.1
Python version: 3.12.13
Platform info:  Linux-4.18.0-553.123.1.el8_10.x86_64-x86_64-with-glibc2.28
Default device: gpu
Created at:     2026-06-16T19:32:29.339498+00:00


CPU Pomp Object Metadata:

Model Initialization Metadata:
------------------------------
pypomp version: 0.4.6.0
JAX version:    0.9.0.1
Python version: 3.12.13
Platform info:  Linux-4.18.0-553.123.1.el8_10.x86_64-x86_64-with-glibc2.28
Default device: cpu
Created at:     2026-06-16T19:32:30.337740+00:00

Performance & Throughput Comparison

The table below presents the execution runtimes, speedup factors, and throughput of each configuration relative to R pomp. The GPU used is a NVIDIA V100, and the CPU used is a 3.0 GHz Intel Xeon Gold 6154.

Configuration	IF2 Time	IF2 Speedup	Pfilter Time	Pfilter Speedup	Total Time	Total Speedup	Throughput (vs 1 CPU Core)
R pomp (CPU, 36 Cores)	6893.52s (114.89m)	1.00x	289.63s (4.83m)	1.00x	7183.41s (119.72m)	1.00x	36.00x
pypomp (CPU, 36 Cores)	3226.21s (53.77m)	2.14x	141.54s (2.36m)	2.05x	3367.75s (56.13m)	2.13x	76.79x
pypomp (GPU, 1 V100 GPU)	372.22s (6.20m)	18.52x	40.59s (0.68m)	7.14x	412.82s (6.88m)	17.40x	626.43x

Performance Expectations

GPU Acceleration: We expect pypomp on the V100 GPU to achieve at least a 16x speedup over the 36-core CPU reference R run.
CPU JAX Overhead: pypomp on the CPU (using JAX) achieve at least a 2x speedup over the R run.
Throughput: Throughput measures the equivalent serial CPU core count required to match the execution speed. A GPU run should demonstrate throughput equivalent to hundreds of CPU cores.

Parameter Convergence & Traces (IF2)

Comparison of Python and R Parameter Traces

Here, we overlay the trajectories from R, Python (CPU), and Python (GPU) to compare the search paths directly.

Expectation

Since pypomp is designed to be a direct Python/JAX implementation of R’s pomp package, the paths taken by the individual replicates in Python (both CPU and GPU) should be similar to those in R, overlapping substantially.

Comparison of Python and R Parameter Trace Quantiles

To abstract away the noise of individual runs, this plot compares the 10th and 90th percentile bounds of the parameters across all runs.

Expectation

While comparing individual noisy trajectories can be difficult, the 10th and 90th percentile ribbons provide a more robust view of convergence. If the package is working properly, the ribbons for Python (GPU), Python (CPU), and R should overlap closely throughout the iterations.

Parameter Estimates Distribution

Comparison of Final Parameter Estimates

This plot compares the distribution of final estimated parameters across replicates.

Expectation

The density of the final parameter values (at the last iteration) across replicates should overlap heavily for each configuration.

Comparison of Starting Parameter Values (Sanity Check)

This plot shows the distribution of initial parameter values to verify that the configurations started from identical distributions.

Expectation

This density plot serves as a sanity check. Because the CPU, GPU, and R runs are initialized using the same random starting box, the starting densities should be be very similar across the configurations.

Log-Likelihood Comparisons

Summary of Final Log-Likelihood Estimates

The table below summarizes the log-likelihood (LL) statistics across replicates at the final iteration.

Configuration	Min	Median	Mean	Max	SD
pypomp (GPU)	1.082571e+04	11841.80	1.178364e+04	11850.08	2.335500e+02
pypomp (CPU)	-9.814143e+30	11841.58	-2.726151e+28	11852.80	5.172507e+29
R pomp	-7.886362e+28	11841.38	-2.190656e+26	11849.36	4.156478e+27

View Top 10 Best Replicates for Each Configuration

Python (GPU) Top 10:

         LL        sd
1  11850.08 0.5501608
2  11849.95 0.8087098
3  11849.69 0.8691682
4  11849.62 0.9652821
5  11849.35 0.6062785
6  11849.21 1.3459180
7  11849.09 1.3932870
8  11848.95 0.7935572
9  11848.84 1.3342843
10 11848.76 0.4291984

Python (CPU) Top 10:

         LL        sd
1  11852.80 5.4839632
2  11852.40 2.9185789
3  11850.47 4.8009773
4  11850.13 3.8902133
5  11849.47 0.6087979
6  11849.36 0.5509179
7  11849.35 0.4042187
8  11849.17 0.6684040
9  11849.04 1.3936672
10 11848.88 1.2137035

R Top 10:

                est        se
result.155 11849.36 0.3611692
result.292 11849.29 2.3713004
result.184 11849.03 0.6089239
result.212 11848.95 0.7779408
result.284 11848.90 1.1797977
result.127 11848.84 0.4285825
result.43  11848.72 0.7901084
result.220 11848.65 0.4942647
result.334 11848.53 2.3266285
result.210 11848.37 1.4432874

Density Plot Comparison (LL > 11,800)

This plot shows the density distribution of final log-likelihood estimates.

Expectation

We would like the final log-likelihood estimates from each configuration to be identical, so we want to see all 3 density estimates overlap heavily.

Likelihood Evaluation at Fixed Parameter (rep3600 Check)

This section evaluates the particle filter’s likelihood calculation itself by running it multiple times on a single, fixed parameter vector (Sun (2024) estimates).

Expectation

Again, we want the two density curves to be about the same, overlapping heavily.

Comparison with Sun (2024) Estimates

This table compares the parameter estimates reported in Sun (2024) to our best Python (GPU) parameter estimates.

quantity	param_value	python_value	Weizhe_value	difference
mu	0.0004	-7.7555	-7.9074	0.1519
kappa	0.0315	-3.4580	-3.4609	0.0030
theta	0.0001	-9.1970	-9.0970	-0.1000
xi	0.0023	-6.0827	-6.0880	0.0053
rho	-0.7656	-2.0192	-1.8921	-0.1270
V_0	0.0000	-10.8398	-9.7435	-1.0963

References

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

--- title: "SPX Comparison" subtitle: "Performance, Convergence, and Likelihood Validation: R pomp vs. pypomp" date: now format: html: theme: light: cosmo dark: darkly page-layout: full toc: true toc-depth: 3 code-fold: true code-summary: "Show Code" code-tools: true embed-resources: true bibliography: references.bib engine: knitr --- ## Introduction This report compares the performance, parameter convergence, and likelihood evaluation of the S&P 500 (SPX) stochastic volatility model used in @sun2024 using the `pypomp` package (on CPU and GPU) and the R `pomp` package (on CPU). The SPX model is a stochastic volatility model with the following parameter set: - $\mu$: Expected rate of return of the asset. - $\kappa$: Mean rate of reversion in variance. - $\theta$: Long term average variance. - $\xi$: Volatility coefficient of volatility. - $\rho$: Leverage effect parameter (correlation between return and volatility shocks). - $V_0$: Initial variance. Because the SPX model uses simple random number generation (sampling from a normal distribution) and features a single `rproc` step per observation, it is highly sensitive to framework overhead in the IF2 optimization and particle filter (`pfilter`) algorithms. This benchmark compares runtimes, convergence trajectories, and likelihood distributions to verify the mathematical correctness and efficiency of `pypomp`. --- ## Setup ::: {.callout-note collapse="true"} ### System & Session Information (R) ```{r Load R libraries, echo = FALSE} library(reticulate) library(tidyverse) library(pomp) setwd(file.path(here::here(), "tests", "spx", "performance")) knitr::opts_chunk$set(echo = FALSE, warning = FALSE, message = FALSE) sessionInfo() print(Sys.time()) ``` ::: ```{python Import python modules} import os os.environ["JAX_PLATFORMS"] = "cpu" import logging # Hide the noisy JAX CUDA initialization error on machines without a GPU logging.getLogger("jax._src.xla_bridge").setLevel(logging.CRITICAL) import pickle import numpy as np import pandas as pd import pypomp as pp import jax ``` ```{python Load data in python} GPU_PKL_PATH = "performance/gpu_results/spx_results_rl4.pkl" CPU_PKL_PATH = "performance/cpu_results/spx_results_rl4.pkl" gpu_pkl = pickle.load(open(GPU_PKL_PATH, "rb")) cpu_pkl = pickle.load(open(CPU_PKL_PATH, "rb")) ``` ```{python Make basic frame} def process_pkl(pkl): res = pkl.results() ll_frame = ( pd.DataFrame( { "LL": res["logLik"], "sd": res["se"], } ) .sort_values(by="LL", ascending=False) .reset_index(drop=True) ) traces = pkl.traces() return ll_frame, traces LL_frame_gpu, traces_gpu = process_pkl(gpu_pkl) LL_frame_cpu, traces_cpu = process_pkl(cpu_pkl) ``` ```{r Make traces frame in R} process_traces <- function(tr) { if ("se" %in% colnames(tr)) { tr <- tr |> select(-se) } tr |> rename(rep = theta_idx, iter = iteration) |> pivot_longer( cols = c(-rep, -iter, -logLik, -method), names_to = "quantity", values_to = "param_value" ) } traces_gpu <- process_traces(py$traces_gpu) traces_cpu <- process_traces(py$traces_cpu) ``` ```{r Set custom ggplot theme} # Define premium plot styles theme_premium <- function() { theme_minimal(base_size = 11) + theme( plot.title = element_text( face = "bold", size = 13, color = "#2c3e50", margin = margin(b = 10) ), plot.subtitle = element_text( size = 10, color = "#7f8c8d", margin = margin(b = 15) ), axis.title = element_text( face = "bold", size = 10, color = "#34495e" ), axis.text = element_text(size = 9, color = "#2c3e50"), legend.title = element_text( face = "bold", size = 9, color = "#34495e" ), legend.text = element_text(size = 9, color = "#2c3e50"), legend.position = "bottom", strip.background = element_rect(fill = "#f8f9fa", color = NA), strip.text = element_text( face = "bold", size = 9, color = "#2c3e50" ), panel.grid.major = element_line(color = "#eaeded"), panel.grid.minor = element_line(color = "#f4f6f6") ) } scale_color_premium <- function(...) { scale_color_manual( values = c( "python (GPU)" = "#1abc9c", "python (CPU)" = "#3498db", "R" = "#e74c3c", "python" = "#1abc9c" ), ... ) } scale_fill_premium <- function(...) { scale_fill_manual( values = c( "python (GPU)" = "#1abc9c", "python (CPU)" = "#3498db", "R" = "#e74c3c", "python" = "#1abc9c" ), ... ) } ``` --- ## Metadata & Execution Runtimes ::: {.callout-note collapse="true"} ### System Configurations ```{python} print("GPU Pomp Object Metadata:") gpu_pkl.print_metadata() print("\nCPU Pomp Object Metadata:") cpu_pkl.print_metadata() ``` ::: ### Performance & Throughput Comparison The table below presents the execution runtimes, speedup factors, and throughput of each configuration relative to R pomp. The GPU used is a NVIDIA V100, and the CPU used is a 3.0 GHz Intel Xeon Gold 6154. ```{r R Execution time} load(here::here( "tests", "spx", "performance", "R_results", "search360_hidden", "1d_global_search360.rda" )) # Extract times r_mif <- t.if.box["elapsed"] r_pf <- t.L.box["elapsed"] r_total <- t.box["elapsed"] gpu_time <- py$gpu_pkl$time() gpu_mif <- gpu_time$time[gpu_time$method == "mif"] gpu_pf <- gpu_time$time[gpu_time$method == "pfilter"] gpu_total <- sum(gpu_time$time) cpu_time <- py$cpu_pkl$time() cpu_mif <- cpu_time$time[cpu_time$method == "mif"] cpu_pf <- cpu_time$time[cpu_time$method == "pfilter"] cpu_total <- sum(cpu_time$time) # Speedups (relative to R pomp) gpu_mif_speedup <- r_mif / gpu_mif gpu_pf_speedup <- r_pf / gpu_pf gpu_total_speedup <- r_total / gpu_total cpu_mif_speedup <- r_mif / cpu_mif cpu_pf_speedup <- r_pf / cpu_pf cpu_total_speedup <- r_total / cpu_total # Throughput vs 1 CPU core # R uses 36 cores r_throughput <- 36.0 gpu_throughput <- gpu_total_speedup * 36.0 cpu_throughput <- cpu_total_speedup * 36.0 timing_df <- data.frame( Configuration = c( "R pomp (CPU, 36 Cores)", "pypomp (CPU, 36 Cores)", "pypomp (GPU, 1 V100 GPU)" ), IF2_Time = c( sprintf("%.2fs (%.2fm)", r_mif, r_mif / 60), sprintf("%.2fs (%.2fm)", cpu_mif, cpu_mif / 60), sprintf("%.2fs (%.2fm)", gpu_mif, gpu_mif / 60) ), IF2_Speedup = c( "1.00x", sprintf("%.2fx", cpu_mif_speedup), sprintf("%.2fx", gpu_mif_speedup) ), Pfilter_Time = c( sprintf("%.2fs (%.2fm)", r_pf, r_pf / 60), sprintf("%.2fs (%.2fm)", cpu_pf, cpu_pf / 60), sprintf("%.2fs (%.2fm)", gpu_pf, gpu_pf / 60) ), Pfilter_Speedup = c( "1.00x", sprintf("%.2fx", cpu_pf_speedup), sprintf("%.2fx", gpu_pf_speedup) ), Total_Time = c( sprintf("%.2fs (%.2fm)", r_total, r_total / 60), sprintf("%.2fs (%.2fm)", cpu_total, cpu_total / 60), sprintf("%.2fs (%.2fm)", gpu_total, gpu_total / 60) ), Total_Speedup = c( "1.00x", sprintf("%.2fx", cpu_total_speedup), sprintf("%.2fx", gpu_total_speedup) ), Throughput_vs_1_Core = c( "36.00x", sprintf("%.2fx", cpu_throughput), sprintf("%.2fx", gpu_throughput) ) ) colnames(timing_df) <- c( "Configuration", "IF2 Time", "IF2 Speedup", "Pfilter Time", "Pfilter Speedup", "Total Time", "Total Speedup", "Throughput (vs 1 CPU Core)" ) knitr::kable( timing_df, format = "html", table.attr = "class='table table-striped table-hover'" ) ``` ::: {.callout-note} ### Performance Expectations - **GPU Acceleration**: We expect `pypomp` on the V100 GPU to achieve at least a 16x speedup over the 36-core CPU reference R run. - **CPU JAX Overhead**: `pypomp` on the CPU (using JAX) achieve at least a 2x speedup over the R run. - **Throughput**: Throughput measures the equivalent serial CPU core count required to match the execution speed. A GPU run should demonstrate throughput equivalent to hundreds of CPU cores. ::: --- ## Parameter Convergence & Traces (IF2) ### Comparison of Python and R Parameter Traces ```{r Set best_rep} # Set best_rep to the rep with the largest logLik on the last iteration for GPU last_iter <- max(traces_gpu$iter) best_rep <- traces_gpu |> filter(iter == last_iter, method == "pfilter") |> group_by(rep) |> summarise(final_logLik = logLik[1]) |> arrange(desc(final_logLik)) |> slice(1) |> pull(rep) ``` Here, we overlay the trajectories from R, Python (CPU), and Python (GPU) to compare the search paths directly. ```{r Load R estimates} load(here::here( "tests", "spx", "performance", "R_results", "search360_hidden", "1d_global_search360.rda" )) LL_df_R <- as.data.frame(L.box) traces_long <- lapply(seq_along(if.box), function(x) { tr <- traces(if.box[[x]]) as.data.frame(tr) |> mutate(iter = 1:nrow(tr), rep = x) }) |> bind_rows() -> traces_long transform_traces <- function(df) { df |> mutate( value_T = case_when( quantity == "logLik" ~ param_value, quantity %in% c("mu", "kappa", "theta", "xi", "V_0") ~ log(if_else( param_value > 0, param_value, NA_real_ )), quantity == "rho" ~ log(if_else( abs(param_value) < 1, (1 + param_value) / (1 - param_value), NA_real_ )), TRUE ~ param_value ) ) } traces_long <- traces_long |> pivot_longer( names_to = "quantity", values_to = "param_value", cols = c(-iter, -rep) ) |> mutate(language = "R") |> transform_traces() # Memory management rm(if.box) invisible(gc()) traces_pyr <- bind_rows( traces_gpu |> mutate(language = "python (GPU)") |> transform_traces(), traces_cpu |> mutate(language = "python (CPU)") |> transform_traces(), traces_long ) ``` ```{r Compare python and R traces, fig.width = 11, fig.height = 8} traces_pyr |> filter(quantity != "logLik") |> filter(value_T >= -12, value_T <= 12) |> ggplot( mapping = aes( x = iter, y = value_T, group = interaction(rep, language), color = language ) ) + geom_line(alpha = 0.15) + facet_wrap(vars(quantity), scales = "free") + labs( title = "Comparison of Python and R Parameter Traces", subtitle = "Overlay of all replicate trajectories (transformed values)", x = "Iteration", y = "Transformed Value", color = "Configuration" ) + scale_color_premium() + theme_premium() + guides(color = guide_legend(override.aes = list(alpha = 1))) ``` ::: {.callout-note} ### Expectation Since `pypomp` is designed to be a direct Python/JAX implementation of R's `pomp` package, the paths taken by the individual replicates in Python (both CPU and GPU) should be similar to those in R, overlapping substantially. ::: --- ### Comparison of Python and R Parameter Trace Quantiles To abstract away the noise of individual runs, this plot compares the 10th and 90th percentile bounds of the parameters across all runs. ```{r Compare python and R traces quantiles, fig.width = 11, fig.height = 8} traces_pyr |> filter(quantity != "logLik") |> filter(value_T >= -12, value_T <= 12) |> group_by(iter, quantity, language) |> summarise( ymin = quantile(value_T, 0.1, na.rm = TRUE), ymax = quantile(value_T, 0.9, na.rm = TRUE), .groups = "drop" ) |> ggplot( mapping = aes( x = iter, ymin = ymin, ymax = ymax, fill = language, color = language ) ) + geom_ribbon(alpha = 0.3, color = NA) + geom_line(aes(y = (ymin + ymax) / 2), linetype = "dashed", alpha = 0.7) + facet_wrap(vars(quantity), scales = "free") + labs( title = "Comparison of Python and R Parameter Trace Quantiles", subtitle = "Shaded regions represent the 10th to 90th percentile interval; dashed line is the midpoint", x = "Iteration", y = "Transformed Value Range", fill = "Configuration", color = "Configuration" ) + scale_fill_premium() + scale_color_premium() + theme_premium() ``` ::: {.callout-note} ### Expectation While comparing individual noisy trajectories can be difficult, the 10th and 90th percentile ribbons provide a more robust view of convergence. If the package is working properly, the ribbons for Python (GPU), Python (CPU), and R should overlap closely throughout the iterations. ::: --- ## Parameter Estimates Distribution ### Comparison of Final Parameter Estimates This plot compares the distribution of final estimated parameters across replicates. ```{r Compare python and R estimates, fig.width = 10, fig.height = 7} traces_pyr |> group_by(language, quantity) |> filter(iter == max(iter)) |> summarise( p10 = quantile(value_T, 0.1, na.rm = TRUE), p90 = quantile(value_T, 0.9, na.rm = TRUE), .groups = "drop" ) |> right_join( traces_pyr |> group_by(language) |> filter(iter == max(iter)), by = c("language", "quantity") ) |> filter(value_T >= -10, value_T <= 10) |> ggplot(mapping = aes(x = value_T, fill = language, color = language)) + geom_density(alpha = 0.3) + facet_wrap(vars(quantity), scales = "free") + labs( title = "Comparison of Final Parameter Estimates", subtitle = "Density distribution of parameters at the final IF2 iteration (transformed, zoomed)", x = "Transformed Parameter Value", y = "Density", fill = "Configuration", color = "Configuration" ) + scale_fill_premium() + scale_color_premium() + theme_premium() ``` ::: {.callout-note} ### Expectation The density of the final parameter values (at the last iteration) across replicates should overlap heavily for each configuration. ::: --- ### Comparison of Starting Parameter Values (Sanity Check) This plot shows the distribution of initial parameter values to verify that the configurations started from identical distributions. ```{r Compare python and R estimates (boxplot), fig.width = 10, fig.height = 7} traces_pyr |> filter(quantity != "logLik", quantity != "loglik") |> group_by(language, quantity) |> filter(iter == min(iter)) |> ggplot(mapping = aes(x = value_T, fill = language, color = language)) + geom_density(alpha = 0.3) + facet_wrap(vars(quantity), scales = "free") + labs( title = "Comparison of Initial Parameter Distributions", subtitle = "Density of starting parameter configurations across restarts (transformed, zoomed)", x = "Transformed Parameter Value", y = "Density", fill = "Configuration", color = "Configuration" ) + scale_fill_premium() + scale_color_premium() + theme_premium() ``` ::: {.callout-note} ### Expectation This density plot serves as a sanity check. Because the CPU, GPU, and R runs are initialized using the same random starting box, the starting densities should be be very similar across the configurations. ::: --- ## Log-Likelihood Comparisons ### Summary of Final Log-Likelihood Estimates The table below summarizes the log-likelihood (LL) statistics across replicates at the final iteration. ```{r LL summary statistics} ll_summary_table <- data.frame( Configuration = c("pypomp (GPU)", "pypomp (CPU)", "R pomp"), Min = c(min(py$LL_frame_gpu$LL), min(py$LL_frame_cpu$LL), min(LL_df_R$est)), Median = c( median(py$LL_frame_gpu$LL), median(py$LL_frame_cpu$LL), median(LL_df_R$est) ), Mean = c( mean(py$LL_frame_gpu$LL), mean(py$LL_frame_cpu$LL), mean(LL_df_R$est) ), Max = c(max(py$LL_frame_gpu$LL), max(py$LL_frame_cpu$LL), max(LL_df_R$est)), SD = c(sd(py$LL_frame_gpu$LL), sd(py$LL_frame_cpu$LL), sd(LL_df_R$est)) ) knitr::kable( ll_summary_table, digits = 2, format = "html", table.attr = "class='table table-striped table-hover'" ) ``` ::: {.callout-tip collapse="true"} ### View Top 10 Best Replicates for Each Configuration **Python (GPU) Top 10:** ```{r} head(arrange(py$LL_frame_gpu, desc(LL)), 10) ``` **Python (CPU) Top 10:** ```{r} head(arrange(py$LL_frame_cpu, desc(LL)), 10) ``` **R Top 10:** ```{r} head(arrange(LL_df_R, desc(est)), 10) ``` ::: --- ### Density Plot Comparison (LL > 11,800) This plot shows the density distribution of final log-likelihood estimates. ```{r LL density comparison, fig.width = 8, fig.height = 5} data.frame( GPU = py$LL_frame_gpu$LL, CPU = py$LL_frame_cpu$LL ) |> pivot_longer(everything(), names_to = "language", values_to = "LL") |> mutate(language = paste0("python (", language, ")")) |> bind_rows(data.frame(LL = LL_df_R$est, language = "R")) |> filter(LL > 11800) |> ggplot(aes(x = LL, fill = language, color = language)) + geom_density(alpha = 0.3) + labs( title = "Comparison of Final Log-Likelihood Estimates", subtitle = "Density of estimated log-likelihoods at final parameter sets (LL > 11,800)", x = "Log-Likelihood (LL)", y = "Density", fill = "Configuration", color = "Configuration" ) + scale_fill_premium() + scale_color_premium() + theme_premium() ``` ::: {.callout-note} ### Expectation We would like the final log-likelihood estimates from each configuration to be identical, so we want to see all 3 density estimates overlap heavily. ::: --- ## Likelihood Evaluation at Fixed Parameter (rep3600 Check) This section evaluates the particle filter's likelihood calculation itself by running it multiple times on a single, fixed parameter vector (@sun2024 estimates). ```{python Load rep3600 eval results} EVAL_PKL_PATH = os.path.join("pfilter_check", "py_results", "spx_results_eval.pkl") pkl_eval = pickle.load(open(EVAL_PKL_PATH, "rb")) LL_rep3600 = pkl_eval.results_history[0].logLiks[0].values.tolist() ``` ```{r Load rep3600 R eval and plot combined density, fig.width = 8, fig.height = 5} load(file.path("pfilter_check", "R_results", "spx_results_eval.rda")) # Extract numeric log-likelihoods from L.box regardless of shape LL_rep3600_R <- as.numeric(L.box) data.frame( LL = c(unlist(py$LL_rep3600), LL_rep3600_R), source = c( rep("python", length(py$LL_rep3600)), rep("R", length(LL_rep3600_R)) ) ) |> ggplot(aes(x = LL, fill = source, color = source)) + geom_density(alpha = 0.3) + labs( title = "Likelihood Evaluation at Fixed Parameter (rep3600)", subtitle = "Log-Likelihood density distribution evaluated at a single benchmark parameter set", x = "Log-Likelihood (LL)", y = "Density", fill = "Source", color = "Source" ) + scale_fill_premium() + scale_color_premium() + theme_premium() ``` ::: {.callout-note} ### Expectation Again, we want the two density curves to be about the same, overlapping heavily. ::: --- ::: {.callout-tip collapse="true"} ### Comparison with @sun2024 Estimates This table compares the parameter estimates reported in @sun2024 to our best Python (GPU) parameter estimates. ```{r Sun estimates} rho_transform <- function(x) { out <- log(if_else(abs(x) < 1, (1 + x) / (1 - x), NA_real_)) return(out) } data <- list( Parameter = c("mu", "kappa", "theta", "xi", "rho", "V_0"), Value = c( log(3.68e-4), log(3.14e-2), log(1.12e-4), log(2.27e-3), rho_transform(-7.38e-1), log(7.66e-3^2) ) ) df_weizhe <- as.data.frame(data) traces_gpu |> filter(iter == last_iter, rep == best_rep, method == "pfilter") |> select(quantity, param_value) |> mutate( param_value_T = case_when( quantity == "logLik" ~ param_value, quantity %in% c("mu", "kappa", "theta", "xi", "V_0") ~ log(if_else( param_value > 0, param_value, NA_real_ )), quantity == "rho" ~ log(if_else( abs(param_value) < 1, (1 + param_value) / (1 - param_value), NA_real_ )), TRUE ~ param_value ) ) |> left_join(df_weizhe, by = c("quantity" = "Parameter")) |> mutate(difference = param_value_T - Value) |> rename(Weizhe_value = Value, python_value = param_value_T) |> knitr::kable( digits = 4, format = "html", table.attr = "class='table table-striped table-hover'" ) ``` ::: --- # References