Skip to content

Commit 7a0326e

Browse files
sbfnk-botsbfnk
andcommitted
Revert default alpha to Normal(0, 0.2)
Normal(0, 0.5) breaks sampling on previously stuck seeds: the wider prior lets chains wander into large-alpha regions during warmup where the geometry funnel is bad enough that the centred form can't compensate. Verified empirically — seed=8 went from clean (R-hat=1.022, td=0) at sd=0.2 to catastrophic (R-hat=1020, td_hits=500) at sd=0.5. Document the prior-sensitivity but keep the default conservative. Co-authored-by: sbfnk <sebastian.funk@lshtm.ac.uk>
1 parent 71112fd commit 7a0326e

4 files changed

Lines changed: 16 additions & 12 deletions

File tree

NEWS.md

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -20,7 +20,7 @@
2020
## Model changes
2121

2222
- Improved identifiability of the non-stationary Gaussian process used to model Rt over time. Three changes in `inst/stan`: (1) GP increments are rescaled by `1/sqrt(gp_n)` so the `alpha` hyperparameter controls the standard deviation of the cumulated log-Rt trajectory rather than of the increment, matching what `gp_opts()` already documents; (2) the cumulated GP is centred so that `R0` is the mean Rt over the trajectory rather than the initial value, eliminating the `(R0, drift)` ridge in the joint posterior; (3) the GP coefficients are now sampled in centred form (`eta ~ N(0, diagSPD)` rather than unit-normal eta scaled by `diagSPD`), which avoids the `(alpha, eta)` funnel and gives substantially smoother HMC geometry in the small-`alpha` regime. Together the changes eliminate stuck chains and divergent transitions seen on previously pathological seeds, while sampling roughly 4× faster.
23-
- The default prior on `alpha` in `gp_opts()` has been widened from `Normal(0, 0.01)` to `Normal(0, 0.5)` to compensate for the rescaling and centring above. Under the old parameterisation, `alpha` represented the increment SD and the implied trajectory SD scaled with `sqrt(gp_n)`; under the new parameterisation, `alpha` represents the trajectory SD directly. The data is only weakly informative about `alpha` (a 13× shift in prior median typically moves the posterior median by ~50%), so the wider default lets users with informed beliefs about Rt variability supply tighter priors as needed. Users with custom `alpha` priors should rescale by approximately `sqrt(gp_n) * sqrt(3)` to match the previous prior's implied trajectory SD.
23+
- The default prior on `alpha` in `gp_opts()` has been widened from `Normal(0, 0.01)` to `Normal(0, 0.2)` to compensate for the rescaling and centring above. Under the old parameterisation, `alpha` represented the increment SD and the implied trajectory SD scaled with `sqrt(gp_n)`; under the new parameterisation, `alpha` represents the trajectory SD directly. The data is only weakly informative about `alpha` (a 13× shift in prior median typically moves the posterior median by ~50%), and very wide priors (e.g. sd > 0.3) can allow chains to wander into pathological regions during warmup. Users with custom `alpha` priors should rescale by approximately `sqrt(gp_n) * sqrt(3)` to match the previous prior's implied trajectory SD.
2424

2525
## Bug fixes
2626

R/opts.R

Lines changed: 7 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -457,12 +457,14 @@ backcalc_opts <- function(prior = c("reports", "none", "infections"),
457457
#' parameter of the Gaussian process kernel. After rescaling of the
458458
#' non-stationary GP, this is approximately the standard deviation of the
459459
#' cumulated log-Rt trajectory over the observation window. Defaults to a
460-
#' half-normal distribution with mean 0 and sd 0.5:
461-
#' `Normal(mean = 0, sd = 0.5)` (a lower limit of 0 will be enforced
460+
#' half-normal distribution with mean 0 and sd 0.2:
461+
#' `Normal(mean = 0, sd = 0.2)` (a lower limit of 0 will be enforced
462462
#' automatically to ensure positivity). Note: the data is only weakly
463463
#' informative about `alpha` — the posterior is sensitive to the prior
464-
#' choice. Users with strong domain knowledge of the expected scale of
465-
#' log-Rt variation should consider supplying an informed prior.
464+
#' choice — and very wide priors (e.g. sd > 0.3) may allow chains to
465+
#' wander into pathological regions during warmup. Users with strong
466+
#' domain knowledge of the expected scale of log-Rt variation should
467+
#' consider supplying an informed prior.
466468
#'
467469
#' @param kernel Character string, the type of kernel required. Currently
468470
#' supporting the Matern kernel ("matern"), squared exponential kernel ("se"),
@@ -504,7 +506,7 @@ backcalc_opts <- function(prior = c("reports", "none", "infections"),
504506
gp_opts <- function(basis_prop = 0.2,
505507
boundary_scale = 1.5,
506508
ls = LogNormal(mean = 21, sd = 7, max = 60),
507-
alpha = Normal(mean = 0, sd = 0.5),
509+
alpha = Normal(mean = 0, sd = 0.2),
508510
kernel = c("matern", "se", "ou", "periodic"),
509511
matern_order = 3 / 2,
510512
w0 = 1.0) {

man/gp_opts.Rd

Lines changed: 7 additions & 5 deletions
Some generated files are not rendered by default. Learn more about customizing how changed files appear on GitHub.

tests/testthat/test-gp_opts.R

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,7 @@ test_that("gp_opts returns correct default values", {
22
gp <- gp_opts()
33
expect_equal(gp$basis_prop, 0.2)
44
expect_equal(gp$boundary_scale, 1.5)
5-
expect_equal(gp$alpha, Normal(0, 0.5))
5+
expect_equal(gp$alpha, Normal(0, 0.2))
66
expect_equal(gp$ls, LogNormal(mean = 21, sd = 7, max = 60))
77
expect_equal(gp$kernel, "matern")
88
expect_equal(gp$matern_order, 3 / 2)

0 commit comments

Comments
 (0)