From bdc96ece3e41b34e4d6a4203b77fa2225b8c0bc2 Mon Sep 17 00:00:00 2001 From: Sebastian Funk Date: Thu, 20 Jul 2023 09:58:13 +0100 Subject: [PATCH 1/6] don't defer to Poisson --- inst/stan/functions/observation_model.stan | 24 ++++++++-------------- 1 file changed, 9 insertions(+), 15 deletions(-) diff --git a/inst/stan/functions/observation_model.stan b/inst/stan/functions/observation_model.stan index 3c6f450f5..917bdd74d 100644 --- a/inst/stan/functions/observation_model.stan +++ b/inst/stan/functions/observation_model.stan @@ -54,28 +54,22 @@ void truncation_lp(real[] truncation_mean, real[] truncation_sd, void report_lp(int[] cases, vector reports, real[] rep_phi, real phi_mean, real phi_sd, int model_type, real weight) { - real sqrt_phi = 1e5; if (model_type) { - // the reciprocal overdispersion parameter (phi) + real sqrt_phi; // the reciprocal overdispersion parameter (phi) rep_phi[model_type] ~ normal(phi_mean, phi_sd) T[0,]; sqrt_phi = 1 / sqrt(rep_phi[model_type]); - // defer to poisson if phi is large, to avoid overflow or - // if poisson specified - } - if (sqrt_phi > 1e4) { if (weight == 1) { - cases ~ poisson(reports); - }else{ - target += poisson_lpmf(cases | reports) * weight; + cases ~ neg_binomial_2(reports, sqrt_phi); + } else { + target += neg_binomial_2_lpmf(cases | reports, sqrt_phi) * weight; } } else { if (weight == 1) { - cases ~ neg_binomial_2(reports, sqrt_phi); - }else{ - target += neg_binomial_2_lpmf(cases | reports, sqrt_phi); + cases ~ poisson(reports); + } else { + target += poisson_lpmf(cases | reports) * weight; } } - } // update log likelihood (as above but not vectorised and returning log likelihood) vector report_log_lik(int[] cases, vector reports, @@ -85,11 +79,11 @@ vector report_log_lik(int[] cases, vector reports, real sqrt_phi = 1e5; if (model_type) { // the reciprocal overdispersion parameter (phi) - sqrt_phi = 1 / sqrt(rep_phi[model_type]); + sqrt_phi = 1 / sqrt(rep_phi[model_type]); } // defer to poisson if phi is large, to avoid overflow - if (sqrt_phi > 1e4) { + if (model_type == 0) { for (i in 1:t) { log_lik[i] = poisson_lpmf(cases[i] | reports[i]) * weight; } From e75ec497f4c51f57cba9847b126df1b8f17f40c9 Mon Sep 17 00:00:00 2001 From: Sebastian Funk Date: Thu, 20 Jul 2023 09:59:08 +0100 Subject: [PATCH 2/6] reinstate tests --- tests/testthat/test-estimate_secondary.R | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/tests/testthat/test-estimate_secondary.R b/tests/testthat/test-estimate_secondary.R index f896737f8..6da058f46 100644 --- a/tests/testthat/test-estimate_secondary.R +++ b/tests/testthat/test-estimate_secondary.R @@ -88,12 +88,12 @@ test_that("estimate_secondary can recover simulated parameters", { ) # These tests currently fail indicating the model is not recovering the # simulated parameters. - # expect_equal( - # prev_posterior[, mean], c(1.6, 0.8, 0.3), tolerance = 0.1 - # ) - # expect_equal( - # prev_posterior[, median], c(1.6, 0.8, 0.3), tolerance = 0.1 - # ) + expect_equal( + prev_posterior[, mean], c(1.6, 0.8, 0.3), tolerance = 0.1 + ) + expect_equal( + prev_posterior[, median], c(1.6, 0.8, 0.3), tolerance = 0.1 + ) }) test_that("forecast_secondary can return values from simulated data and plot From c1b076f47c0ffb226f6d3e6ab768904e68a649dd Mon Sep 17 00:00:00 2001 From: Sebastian Funk Date: Thu, 20 Jul 2023 09:59:35 +0100 Subject: [PATCH 3/6] reduce adapt_delta --- tests/testthat/test-estimate_secondary.R | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/testthat/test-estimate_secondary.R b/tests/testthat/test-estimate_secondary.R index 6da058f46..9e11fb2a5 100644 --- a/tests/testthat/test-estimate_secondary.R +++ b/tests/testthat/test-estimate_secondary.R @@ -23,7 +23,6 @@ cases[ # with a secondary case inc <- estimate_secondary(cases[1:60], obs = obs_opts(scale = list(mean = 0.2, sd = 0.2), week_effect = FALSE), - control = list(adapt_delta = 0.98), verbose = FALSE ) From a966ec6b0c46d2be3ab49fb28d8b362341b737a5 Mon Sep 17 00:00:00 2001 From: Sebastian Funk Date: Fri, 21 Jul 2023 07:31:38 +0100 Subject: [PATCH 4/6] update log-likelihood generation --- inst/stan/functions/observation_model.stan | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/inst/stan/functions/observation_model.stan b/inst/stan/functions/observation_model.stan index 917bdd74d..419528374 100644 --- a/inst/stan/functions/observation_model.stan +++ b/inst/stan/functions/observation_model.stan @@ -76,11 +76,6 @@ vector report_log_lik(int[] cases, vector reports, real[] rep_phi, int model_type, real weight) { int t = num_elements(reports); vector[t] log_lik; - real sqrt_phi = 1e5; - if (model_type) { - // the reciprocal overdispersion parameter (phi) - sqrt_phi = 1 / sqrt(rep_phi[model_type]); - } // defer to poisson if phi is large, to avoid overflow if (model_type == 0) { @@ -88,10 +83,11 @@ vector report_log_lik(int[] cases, vector reports, log_lik[i] = poisson_lpmf(cases[i] | reports[i]) * weight; } } else { + real sqrt_phi = 1 / sqrt(rep_phi[model_type]); for (i in 1:t) { log_lik[i] = neg_binomial_2_lpmf(cases[i] | reports[i], sqrt_phi) * weight; } - } + } return(log_lik); } // sample reported cases from the observation model From 9c57d60538128abca78a4c1ba343c09e550d598d Mon Sep 17 00:00:00 2001 From: Sebastian Funk Date: Fri, 21 Jul 2023 07:32:12 +0100 Subject: [PATCH 5/6] remove outdated comments --- tests/testthat/test-estimate_secondary.R | 2 -- 1 file changed, 2 deletions(-) diff --git a/tests/testthat/test-estimate_secondary.R b/tests/testthat/test-estimate_secondary.R index 9e11fb2a5..95fc48702 100644 --- a/tests/testthat/test-estimate_secondary.R +++ b/tests/testthat/test-estimate_secondary.R @@ -85,8 +85,6 @@ test_that("estimate_secondary can recover simulated parameters", { inc_posterior[, median], c(1.8, 0.5, 0.4), tolerance = 0.1 ) - # These tests currently fail indicating the model is not recovering the - # simulated parameters. expect_equal( prev_posterior[, mean], c(1.6, 0.8, 0.3), tolerance = 0.1 ) From 776ceb8294294fc0ce3c4e6578739745e3c67dc7 Mon Sep 17 00:00:00 2001 From: Sebastian Funk Date: Fri, 21 Jul 2023 07:32:33 +0100 Subject: [PATCH 6/6] increase tolerance for parameter recovery --- tests/testthat/test-estimate_secondary.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/testthat/test-estimate_secondary.R b/tests/testthat/test-estimate_secondary.R index 95fc48702..223ecca18 100644 --- a/tests/testthat/test-estimate_secondary.R +++ b/tests/testthat/test-estimate_secondary.R @@ -86,10 +86,10 @@ test_that("estimate_secondary can recover simulated parameters", { tolerance = 0.1 ) expect_equal( - prev_posterior[, mean], c(1.6, 0.8, 0.3), tolerance = 0.1 + prev_posterior[, mean], c(1.6, 0.8, 0.3), tolerance = 0.2 ) expect_equal( - prev_posterior[, median], c(1.6, 0.8, 0.3), tolerance = 0.1 + prev_posterior[, median], c(1.6, 0.8, 0.3), tolerance = 0.2 ) })