Merge 776ceb8 into 1205d1a

epiforecasts · Jul 21, 2023 · 2833027 · 2833027
2 parents 1205d1a + 776ceb8
commit 2833027
Show file tree

Hide file tree

Showing 2 changed files with 16 additions and 29 deletions.
diff --git a/inst/stan/functions/observation_model.stan b/inst/stan/functions/observation_model.stan
@@ -54,50 +54,40 @@ 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,
                       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 (sqrt_phi > 1e4) {
+  if (model_type == 0) {
     for (i in 1:t) {
       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

diff --git 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
 )
 
@@ -86,14 +85,12 @@ 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
-  # )
-  # 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.2
+  )
+  expect_equal(
+    prev_posterior[, median], c(1.6, 0.8, 0.3), tolerance = 0.2
+  )
 })
 
 test_that("forecast_secondary can return values from simulated data and plot