2024-08-07 update : adding matrix normalization func. and using AR1 i…

…n spatial func. (#53)

R/generate_simulated_data.R

@@ -62,6 +62,11 @@
 #' @param scaling_factor Scaling factor for aux site
 #' @param aux_site_bool Boolean to use the aux site framework with
 #' scaling factor.
+#' @param init_stat Boolean. Should the initial value of the AR(1) be drawn
+#' from the process's stationary distribution (`TRUE`) or from the process's
+#' conditional error distribution (`FALSE`)? Note that the process only has
+#' a defined stationary distribution if `phi_rt` < 1.
+#' Default `TRUE`.
 #'
 #' @return a list containing three dataframes. hosp_data is a dataframe
 #' containing the number of daily hospital admissions by day for a theoretical
@@ -101,7 +106,8 @@
 #'   phi_rt = 0.6,
 #'   sigma_eps = sqrt(0.02),
 #'   scaling_factor = 0.01,
-#'   aux_site_bool = TRUE
+#'   aux_site_bool = TRUE,
+#'   init_stat = TRUE
 #' )
 #' hosp_data <- sim_data$hosp_data
 #' ww_data <- sim_data$ww_data
@@ -161,7 +167,8 @@ generate_simulated_data <- function(r_in_weeks = # nolint
                                     phi_rt = 0.6,
                                     sigma_eps = sqrt(0.02),
                                     scaling_factor = 0.01,
-                                    aux_site_bool = TRUE) {
+                                    aux_site_bool = TRUE,
+                                    init_stat = TRUE) {
   # Some quick checks to make sure the inputs work as expected
   stopifnot(
     "weekly R(t) passed in isn't long enough" =
@@ -366,19 +373,17 @@ generate_simulated_data <- function(r_in_weeks = # nolint
   )
   # Auxiliary Site
   if (aux_site_bool) {
-    state_deviation_noise_vec <- state_deviation_noise_vec_aux_rng(
-      scaling_factor,
-      sigma_eps,
-      n_weeks
-    )
-    state_deviation_init <- rnorm(
-      n = 1,
+    state_deviation_noise_vec <- rnorm(
+      n = n_weeks,
       mean = 0,
-      sd = scaling_factor * sigma_eps
+      sd = 1
     )
     log_r_site_aux <- construct_aux_rt(log_r_state_week,
       state_deviation_ar_coeff = phi_rt,
-      state_deviation_noise_vec
+      scaling_factor,
+      sigma_eps,
+      state_deviation_noise_vec,
+      init_stat
     )
     log_r_site <- rbind(
       log_r_site,

inst/stan/functions/aux_site_process_rng.stan

@@ -6,23 +6,33 @@
   * of spatial epsilon.
   * @param state_deviation_ar_coeff Coefficient for AR(1) temporal correlation on
   * state deviations for aux site.
+  * @param init_stat Boolean. Should the initial value of the AR(1) be drawn
+  * from the process's stationary stationary distribution (1) or from the
+  * process's conditional error distribution (0)? Note that the process only has
+  * a defined stationary distribution if `state_deviation_ar_coeff` < 1.
   * @return A vector for unadjusted Rt for aux site where columns are time points.
   */
 vector aux_site_process_rng(vector log_state_rt,
                               real scaling_factor,
                               real sigma_eps,
-                              real state_deviation_ar_coeff) {
+                              real state_deviation_ar_coeff,
+                              int init_stat) {
 
 
   //Presets
   int n_time = size(log_state_rt);
 
-  vector[n_time] state_deviation_noise_vec = state_deviation_noise_vec_aux_rng(scaling_factor,
-                                                                               sigma_eps,
-                                                                               n_time);
+  // creates vector of standard normal for contruct aux rt.
+  vector[n_time] z;
+  for (i in 1:n_time){
+    z[i] = std_normal_rng();
+  }
   vector[n_time] log_aux_site_rt = construct_aux_rt(log_state_rt,
                                                     state_deviation_ar_coeff,
-                                                    state_deviation_noise_vec);
+                                                    scaling_factor,
+                                                    sigma_eps,
+                                                    z,
+                                                    init_stat);
 
   return log_aux_site_rt;
 }

inst/stan/functions/construct_aux_rt.stan

@@ -3,31 +3,30 @@
   * @param log_state_rt "State" level unadjusted Rt, on log scale.
   * @param state_deviation_ar_coeff Coefficient for AR(1) temporal correlation on
   * subpopulation deviations.
-  * @param state_deviation_init Initial value of the spatial deviation
-  * @param state_deviation_noise_vec vector with random vectors of noise
-  * for state devaiation process for aux site.
+  * @param scale_factor Scaling factor for aux site process.
+  * @param init_stat Boolean. Should the initial value of the AR(1) be drawn
+  * from the process's stationary stationary distribution (1) or from the
+  * process's conditional error distribution (0)? Note that the process only has
+  * a defined stationary distribution if `state_deviation_ar_coeff` < 1.
   * @return A vector for unadjusted Rt for aux site where columns are time points.
   */
 vector construct_aux_rt(vector log_state_rt,
                         real state_deviation_ar_coeff,
-                        vector state_deviation_noise_vec) {
+                        real scaling_factor,
+                        real sigma_eps,
+                        vector z,
+                        int init_stat) {
 
 
   // presets
-  int n_time = dims(state_deviation_noise_vec)[1];
-  vector[n_time] log_aux_site_rt;
-  real state_deviation_t_i = 0;
-  // To explain state_deviation_t_i = 0 initialization,
-  // we need a state_deviation_t_i value for t_1 = 1, so we will
-  // initialize the process with index 1 of spatial
-  // deviation noise matrix w/out spatial_deviation ar coeff, so we set
-  // state_deviation_t_i to 0 initially.
-
-  for (t_i in 1:n_time) {
-    state_deviation_t_i = state_deviation_ar_coeff * state_deviation_t_i +
-        state_deviation_noise_vec[t_i];
-    log_aux_site_rt[t_i] = log_state_rt[t_i] + state_deviation_t_i;
-  }
+  int n_time = dims(z)[1];
+  vector[n_time] log_aux_site_rt = ar1(
+    log_state_rt,
+    state_deviation_ar_coeff,
+    scaling_factor * sigma_eps,
+    z,
+    init_stat
+  );
 
   return log_aux_site_rt;
 }

inst/stan/functions/matrix_normalization.stan

@@ -0,0 +1,11 @@
+/**
+  * Normalizes a matrix using its determinant to the power of 1 / n.
+  * @param matrx A matrix, usually the correlation matrix need to normalize
+  * @return A matrix that has been normalized.
+  *
+  */
+matrix matrix_normalization(matrix matrx) {
+  int n = cols(matrx);
+  matrix[n,n] norm_matrx = matrx / pow(determinant(matrx), 1.0 / n);
+  return norm_matrx;
+}

inst/stan/functions/spatial_functions.stan

@@ -2,6 +2,7 @@
   * File to call all desired spatial functions for use with expose_functions()
   */
 functions{
+  #include ar1.stan
   #include exponential_decay_corr_func.stan
   #include independence_corr_func.stan
   #include construct_spatial_rt.stan
@@ -10,4 +11,5 @@ functions{
   #include construct_aux_rt.stan
   #include state_deviation_noise_vec_aux_rng.stan
   #include aux_site_process_rng.stan
+  #include matrix_normalization.stan
 }