From 923051c05e8f2f194efcba35665a1107d1e164c4 Mon Sep 17 00:00:00 2001 From: Sam Abbott Date: Wed, 25 Oct 2023 14:22:54 +0100 Subject: [PATCH] Update phi parameterisation (#487) * Update observation_model.stan * Update NEWS.md * Update inst/stan/functions/observation_model.stan Co-authored-by: Sebastian Funk * Update inst/stan/functions/observation_model.stan Co-authored-by: Sebastian Funk * Update inst/stan/functions/observation_model.stan Co-authored-by: Sebastian Funk * Update inst/stan/functions/observation_model.stan Co-authored-by: Sebastian Funk * Update inst/stan/functions/observation_model.stan Co-authored-by: Sebastian Funk * Update inst/stan/functions/observation_model.stan Co-authored-by: Sebastian Funk * Update inst/stan/functions/observation_model.stan Co-authored-by: Sebastian Funk * sqrt_phi -> dispersion * Catch another sqrt_phi * Update NEWS.md Co-authored-by: Sebastian Funk --------- Co-authored-by: Sebastian Funk --- NEWS.md | 4 ++++ inst/stan/functions/observation_model.stan | 19 +++++++++---------- 2 files changed, 13 insertions(+), 10 deletions(-) diff --git a/NEWS.md b/NEWS.md index eee0a5c1f..fabb049d2 100644 --- a/NEWS.md +++ b/NEWS.md @@ -16,6 +16,10 @@ * Changed all instances of arguments that refer to the maximum of a distribution to reflect the maximum. Previously this did, in some instance, refer to the length of the PMF. By @sbfnk in #468. * Fixed a bug in the bounds of delays when setting initial conditions. By @sbfnk in #474. +## Model changes + +* Updated the parameterisation of the dispersion term `phi` to be `phi = 1 / sqrt_phi ^ 2` rather than the previous parameterisation `phi = 1 / sqrt(sqrt_phi)` based on the suggested prior [here](https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations#story-when-the-generic-prior-fails-the-case-of-the-negative-binomial) and the performance benefits seen in the `epinowcast` package (see [here](https://github.com/epinowcast/epinowcast/blob/8eff560d1fd8305f5fb26c21324b2bfca1f002b4/inst/stan/epinowcast.stan#L314)). By @seabbs in # and reviewed by @sbfnk. + # EpiNow2 1.4.0 This release contains some bug fixes, minor new features, and the initial stages of some broader improvement to future handling of delay distributions. diff --git a/inst/stan/functions/observation_model.stan b/inst/stan/functions/observation_model.stan index 3968d3aa8..ed3364ae6 100644 --- a/inst/stan/functions/observation_model.stan +++ b/inst/stan/functions/observation_model.stan @@ -55,13 +55,12 @@ void report_lp(array[] int cases, vector reports, array[] real rep_phi, real phi_mean, real phi_sd, int model_type, real weight) { if (model_type) { - real sqrt_phi; // the reciprocal overdispersion parameter (phi) + real dispersion = 1 / pow(rep_phi[model_type], 2); rep_phi[model_type] ~ normal(phi_mean, phi_sd) T[0,]; - sqrt_phi = 1 / sqrt(rep_phi[model_type]); if (weight == 1) { - cases ~ neg_binomial_2(reports, sqrt_phi); + cases ~ neg_binomial_2(reports, dispersion); } else { - target += neg_binomial_2_lpmf(cases | reports, sqrt_phi) * weight; + target += neg_binomial_2_lpmf(cases | reports, dispersion) * weight; } } else { if (weight == 1) { @@ -83,9 +82,9 @@ vector report_log_lik(array[] int cases, vector reports, log_lik[i] = poisson_lpmf(cases[i] | reports[i]) * weight; } } else { - real sqrt_phi = 1 / sqrt(rep_phi[model_type]); + real dispersion = 1 / pow(rep_phi[model_type], 2); for (i in 1:t) { - log_lik[i] = neg_binomial_2_lpmf(cases[i] | reports[i], sqrt_phi) * weight; + log_lik[i] = neg_binomial_2_lpmf(cases[i] | reports[i], dispersion) * weight; } } return(log_lik); @@ -94,9 +93,9 @@ vector report_log_lik(array[] int cases, vector reports, array[] int report_rng(vector reports, array[] real rep_phi, int model_type) { int t = num_elements(reports); array[t] int sampled_reports; - real sqrt_phi = 1e5; + real dispersion = 1e5; if (model_type) { - sqrt_phi = 1 / sqrt(rep_phi[model_type]); + dispersion = 1 / pow(rep_phi[model_type], 2); } for (s in 1:t) { @@ -104,10 +103,10 @@ array[] int report_rng(vector reports, array[] real rep_phi, int model_type) { sampled_reports[s] = 0; } else { // defer to poisson if phi is large, to avoid overflow - if (sqrt_phi > 1e4) { + if (dispersion > 1e4) { sampled_reports[s] = poisson_rng(reports[s] > 1e8 ? 1e8 : reports[s]); } else { - sampled_reports[s] = neg_binomial_2_rng(reports[s] > 1e8 ? 1e8 : reports[s], sqrt_phi); + sampled_reports[s] = neg_binomial_2_rng(reports[s] > 1e8 ? 1e8 : reports[s], dispersion); } } }