fixes for periodic kernel dimension differences

epiforecasts · Aug 15, 2024 · e0c3849 · e0c3849
1 parent ecf8965
commit e0c3849
Show file tree

Hide file tree

Showing 4 changed files with 13 additions and 5 deletions.
diff --git a/R/create.R b/R/create.R
@@ -610,7 +610,8 @@ create_initial_conditions <- function(data) {
     out <- create_delay_inits(data)
 
     if (data$fixed == 0) {
-      out$eta <- array(rnorm(data$M, mean = 0, sd = 0.1))
+      out$eta <- array(rnorm(
+        ifelse(data$gp_type == 1, data$M * 2, data$M), mean = 0, sd = 0.1))
       out$rescaled_rho <- array(rlnorm(1,
         meanlog = data$ls_meanlog,
         sdlog = ifelse(data$ls_sdlog > 0, data$ls_sdlog, 0.01)
@@ -621,7 +622,7 @@ create_initial_conditions <- function(data) {
           out$rescaled_rho < data$ls_min, data$ls_min + 0.001,
           default = out$rescaled_rho
           ))
-      }else{
+      }else {
         out$rescaled_rho <- array(numeric(0))
       }
 

diff --git a/inst/stan/estimate_infections.stan b/inst/stan/estimate_infections.stan
@@ -26,7 +26,7 @@ transformed data {
   int noise_terms = setup_noise(
     ot_h, t, horizon, estimate_r, stationary, future_fixed, fixed_from
   );
-  matrix[noise_terms, M] PHI = setup_gp(M, L, noise_terms, gp_type == 1, w0);  // basis function
+  matrix[noise_terms, gp_type == 1 ? 2*M : M] PHI = setup_gp(M, L, noise_terms, gp_type == 1, w0);  // basis function
   // Rt
   real r_logmean = log(r_mean^2 / sqrt(r_sd^2 + r_mean^2));
   real r_logsd = sqrt(log(1 + (r_sd^2 / r_mean^2)));
@@ -44,7 +44,7 @@ parameters {
   // gaussian process
   array[fixed || gp_type == 3 ? 0 : 1] real<lower = ls_min, upper = ls_max> rescaled_rho;  // length scale of noise GP
   array[fixed ? 0 : 1] real<lower = 0> alpha;    // scale of noise GP
-  vector[fixed ? 0 : M] eta;               // unconstrained noise
+  vector[fixed ? 0 : gp_type == 1 ? 2*M : M] eta;               // unconstrained noise
   // Rt
   vector[estimate_r] log_R;                // baseline reproduction number estimate (log)
   array[estimate_r] real initial_infections;    // seed infections

diff --git a/inst/stan/functions/gaussian_process.stan b/inst/stan/functions/gaussian_process.stan
@@ -148,7 +148,7 @@ matrix setup_gp(int M, real L, int dimension, int is_periodic, real w0) {
   */
 vector update_gp(matrix PHI, int M, real L, real alpha,
                   array[] real rho, vector eta, int type, real nu) {
-  vector[M] diagSPD;    // spectral density
+  vector[type == 1 ? 2 * M : M] diagSPD;    // spectral density
 
   // GP in noise - spectral densities
   if (type == 0) {

diff --git a/tests/testthat/test-estimate_infections.R b/tests/testthat/test-estimate_infections.R
@@ -52,6 +52,13 @@ test_that("estimate_infections successfully returns estimates using a Linear ker
   )
 })
 
+test_that("estimate_infections successfully returns estimates using a periodic kernel", {
+  skip_on_cran()
+  test_estimate_infections(
+    reported_cases, gp = gp_opts(kernel = "periodic")
+  )
+})
+
 test_that("estimate_infections successfully returns estimates when passed NA values", {
   skip_on_cran()
   reported_cases_na <- data.table::copy(reported_cases)