diff --git a/DESCRIPTION b/DESCRIPTION index 87c107bdc..980aee604 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -109,7 +109,7 @@ Imports: R.utils (>= 2.0.0), Rcpp (>= 0.12.0), rlang (>= 0.4.7), - rstan (>= 2.21.1), + rstan (>= 2.26.0), rstantools (>= 2.2.0), runner, scales, @@ -135,8 +135,8 @@ LinkingTo: Rcpp (>= 0.12.0), RcppEigen (>= 0.3.3.3.0), RcppParallel (>= 5.0.1), - rstan (>= 2.21.1), - StanHeaders (>= 2.21.0-5) + rstan (>= 2.26.0), + StanHeaders (>= 2.26.0) Biarch: true Config/testthat/edition: 3 Encoding: UTF-8 diff --git a/inst/stan/data/delays.stan b/inst/stan/data/delays.stan index 06a3d8317..5fa8207b8 100644 --- a/inst/stan/data/delays.stan +++ b/inst/stan/data/delays.stan @@ -1,18 +1,18 @@ int delay_n; // number of delay distributions int delay_n_p; // number of parametric delay distributions int delay_n_np; // number of nonparametric delay distributions - real delay_mean_mean[delay_n_p]; // prior mean of mean delay distribution - real delay_mean_sd[delay_n_p]; // prior sd of mean delay distribution - real delay_sd_mean[delay_n_p]; // prior sd of sd of delay distribution - real delay_sd_sd[delay_n_p]; // prior sd of sd of delay distribution - int delay_max[delay_n_p]; // maximum delay distribution - int delay_dist[delay_n_p]; // 0 = lognormal; 1 = gamma + array[delay_n_p] real delay_mean_mean; // prior mean of mean delay distribution + array[delay_n_p] real delay_mean_sd; // prior sd of mean delay distribution + array[delay_n_p] real delay_sd_mean; // prior sd of sd of delay distribution + array[delay_n_p] real delay_sd_sd; // prior sd of sd of delay distribution + array[delay_n_p] int delay_max; // maximum delay distribution + array[delay_n_p] int delay_dist; // 0 = lognormal; 1 = gamma int delay_np_pmf_max; // number of nonparametric pmf elements vector[delay_np_pmf_max] delay_np_pmf; // ragged array of fixed PMFs - int delay_np_pmf_groups[delay_n_np + 1]; // links to ragged array - int delay_weight[delay_n_p]; + array[delay_n_np + 1] int delay_np_pmf_groups; // links to ragged array + array[delay_n_p] int delay_weight; int delay_types; // number of delay types - int delay_types_p[delay_n]; // whether delay types are parametric - int delay_types_id[delay_n]; // whether delay types are parametric - int delay_types_groups[delay_types + 1]; // index of each delay (parametric or non) + array[delay_n] int delay_types_p; // whether delay types are parametric + array[delay_n] int delay_types_id; // whether delay types are parametric + array[delay_types + 1] int delay_types_groups; // index of each delay (parametric or non) diff --git a/inst/stan/data/generation_time.stan b/inst/stan/data/generation_time.stan index 6b372b078..73e5ef64b 100644 --- a/inst/stan/data/generation_time.stan +++ b/inst/stan/data/generation_time.stan @@ -1,8 +1,8 @@ - real gt_mean_sd[1]; // prior sd of mean generation time - real gt_mean_mean[1]; // prior mean of mean generation time - real gt_sd_mean[1]; // prior mean of sd of generation time - real gt_sd_sd[1]; // prior sd of sd of generation time - int gt_max[1]; // maximum generation time - int gt_fixed[1]; // 0 = variable gt; 1 = fixed gt - int gt_dist[1]; // distribution (0 = lognormal, 1 = gamma) + array[1] real gt_mean_sd; // prior sd of mean generation time + array[1] real gt_mean_mean; // prior mean of mean generation time + array[1] real gt_sd_mean; // prior mean of sd of generation time + array[1] real gt_sd_sd; // prior sd of sd of generation time + array[1] int gt_max; // maximum generation time + array[1] int gt_fixed; // 0 = variable gt; 1 = fixed gt + array[1] int gt_dist; // distribution (0 = lognormal, 1 = gamma) int gt_weight; diff --git a/inst/stan/data/observation_model.stan b/inst/stan/data/observation_model.stan index c5a7f5085..0ce9ef3bb 100644 --- a/inst/stan/data/observation_model.stan +++ b/inst/stan/data/observation_model.stan @@ -1,4 +1,4 @@ - int day_of_week[t - seeding_time]; // day of the week indicator (1 - 7) + array[t - seeding_time] int day_of_week; // day of the week indicator (1 - 7) int model_type; // type of model: 0 = poisson otherwise negative binomial real phi_mean; // Mean and sd of the normal prior for the real phi_sd; // reporting process diff --git a/inst/stan/data/observations.stan b/inst/stan/data/observations.stan index 845a6751a..654304f85 100644 --- a/inst/stan/data/observations.stan +++ b/inst/stan/data/observations.stan @@ -2,5 +2,5 @@ int seeding_time; // time period used for seeding and not observed int horizon; // forecast horizon int future_time; // time in future for Rt - int cases[t - horizon - seeding_time]; // observed cases + array[t - horizon - seeding_time] int cases; // observed cases vector[t] shifted_cases; // prior infections (for backcalculation) diff --git a/inst/stan/data/rt.stan b/inst/stan/data/rt.stan index 0b00c2af6..11b1989ae 100644 --- a/inst/stan/data/rt.stan +++ b/inst/stan/data/rt.stan @@ -4,7 +4,7 @@ real r_mean; // prior mean of reproduction number real r_sd; // prior standard deviation of reproduction number int bp_n; // no of breakpoints (0 = no breakpoints) - int breakpoints[t - seeding_time]; // when do breakpoints occur + array[t - seeding_time] int breakpoints; // when do breakpoints occur int future_fixed; // is underlying future Rt assumed to be fixed int fixed_from; // Reference date for when Rt estimation should be fixed int pop; // Initial susceptible population diff --git a/inst/stan/data/simulation_delays.stan b/inst/stan/data/simulation_delays.stan index 7e5251626..0f3fe7d9e 100644 --- a/inst/stan/data/simulation_delays.stan +++ b/inst/stan/data/simulation_delays.stan @@ -1,18 +1,18 @@ int delay_n; // number of delay distribution distributions int delay_n_p; // number of parametric delay distributions int delay_n_np; // number of nonparametric delay distributions - real delay_mean[n, delay_n_p]; // prior mean of mean delay distribution - real delay_sd[n, delay_n_p]; // prior sd of sd of delay distribution - int delay_max[delay_n_p]; // maximum delay distribution - int delay_dist[delay_n_p]; // 0 = lognormal; 1 = gamma + array[n, delay_n_p] real delay_mean; // prior mean of mean delay distribution + array[n, delay_n_p] real delay_sd; // prior sd of sd of delay distribution + array[delay_n_p] int delay_max; // maximum delay distribution + array[delay_n_p] int delay_dist; // 0 = lognormal; 1 = gamma int delay_np_pmf_max; // number of nonparametric pmf elements vector[delay_np_pmf_max] delay_np_pmf; // ragged array of fixed PMFs - int delay_np_pmf_groups[delay_n_np + 1]; // links to ragged array - int delay_weight[delay_n_p]; + array[delay_n_np + 1] int delay_np_pmf_groups; // links to ragged array + array[delay_n_p] int delay_weight; int delay_types; // number of delay types - int delay_types_p[delay_n]; // whether delay types are parametric - int delay_types_id[delay_n]; // whether delay types are parametric - int delay_types_groups[delay_types + 1]; // index of each delay (parametric or non) + array[delay_n] int delay_types_p; // whether delay types are parametric + array[delay_n] int delay_types_id; // whether delay types are parametric + array[delay_types + 1] int delay_types_groups; // index of each delay (parametric or non) int delay_id; // id of generation time diff --git a/inst/stan/data/simulation_observation_model.stan b/inst/stan/data/simulation_observation_model.stan index a42f6529f..c8cab6b35 100644 --- a/inst/stan/data/simulation_observation_model.stan +++ b/inst/stan/data/simulation_observation_model.stan @@ -1,8 +1,8 @@ - int day_of_week[t - seeding_time]; // day of the week indicator (1 - 7) + array[t - seeding_time] int day_of_week; // day of the week indicator (1 - 7) int week_effect; // should a day of the week effect be estimated - real day_of_week_simplex[n, week_effect]; + array[n, week_effect] real day_of_week_simplex; int obs_scale; - real frac_obs[n, obs_scale]; + array[n, obs_scale] real frac_obs; int model_type; - real rep_phi[n, model_type]; // overdispersion of the reporting process + array[n, model_type] real rep_phi; // overdispersion of the reporting process int trunc_id; // id of truncation diff --git a/inst/stan/data/simulation_rt.stan b/inst/stan/data/simulation_rt.stan index 93b1206a1..a97300451 100644 --- a/inst/stan/data/simulation_rt.stan +++ b/inst/stan/data/simulation_rt.stan @@ -1,5 +1,5 @@ - real initial_infections[seeding_time ? n : 0, 1]; // initial logged infections - real initial_growth[seeding_time > 1 ? n : 0, 1]; //initial growth + array[seeding_time ? n : 0, 1] real initial_infections; // initial logged infections + array[seeding_time > 1 ? n : 0, 1] real initial_growth; //initial growth matrix[n, t - seeding_time] R; // reproduction number int pop; // susceptible population diff --git a/inst/stan/dist_fit.stan b/inst/stan/dist_fit.stan index b247d14f1..3285ff59a 100644 --- a/inst/stan/dist_fit.stan +++ b/inst/stan/dist_fit.stan @@ -3,15 +3,15 @@ data { int N; vector[N] low; vector[N] up; - real lam_mean[dist == 0]; - real prior_mean[dist > 0]; - real prior_sd[dist > 0]; - real par_sigma[dist == 2]; + array[dist == 0] real lam_mean; + array[dist > 0] real prior_mean; + array[dist > 0] real prior_sd; + array[dist == 2] real par_sigma; } transformed data { - real prior_alpha[dist == 2]; - real prior_beta[dist == 2]; + array[dist == 2] real prior_alpha; + array[dist == 2] real prior_beta; if (dist == 2) { prior_alpha[1] = (prior_mean[1] / prior_sd[1])^2; @@ -20,16 +20,16 @@ transformed data { } parameters { - real lambda[dist == 0]; - real mu[dist == 1]; - real sigma[dist == 1]; - real alpha_raw[dist == 2]; - real beta_raw[dist == 2]; + array[dist == 0] real lambda; + array[dist == 1] real mu; + array[dist == 1] real sigma; + array[dist == 2] real alpha_raw; + array[dist == 2] real beta_raw; } transformed parameters{ - real alpha[dist == 2]; - real beta[dist == 2]; + array[dist == 2] real alpha; + array[dist == 2] real beta; if (dist == 2) { alpha[1] = prior_alpha[1] + par_sigma[1] * alpha_raw[1]; @@ -50,19 +50,19 @@ model { for(i in 1:N){ if (dist == 0) { - target += log( - exponential_cdf(up[i] , lambda) - - exponential_cdf(low[i], lambda) + target += log_diff_exp( + exponential_lcdf(up[i] | lambda), + exponential_lcdf(low[i] | lambda) ); } else if (dist == 1) { - target += log( - lognormal_cdf(up[i], mu, sigma) - - lognormal_cdf(low[i], mu, sigma) + target += log_diff_exp( + lognormal_lcdf(up[i] | mu, sigma), + lognormal_lcdf(low[i] | mu, sigma) ); } else if (dist == 2) { - target += log( - gamma_cdf(up[i], alpha, beta) - - gamma_cdf(low[i], alpha, beta) + target += log_diff_exp( + gamma_lcdf(up[i] | alpha, beta), + gamma_lcdf(low[i] | alpha, beta) ); } } diff --git a/inst/stan/estimate_infections.stan b/inst/stan/estimate_infections.stan index 3a4bac547..308256c2d 100644 --- a/inst/stan/estimate_infections.stan +++ b/inst/stan/estimate_infections.stan @@ -30,7 +30,7 @@ transformed data{ 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))); - int delay_type_max[delay_types] = get_delay_type_max( + array[delay_types] int delay_type_max = get_delay_type_max( delay_types, delay_types_p, delay_types_id, delay_types_groups, delay_max, delay_np_pmf_groups ); @@ -38,21 +38,21 @@ transformed data{ parameters{ // gaussian process - real rho[fixed ? 0 : 1]; // length scale of noise GP - real alpha[fixed ? 0 : 1]; // scale of of noise GP + array[fixed ? 0 : 1] real rho; // length scale of noise GP + array[fixed ? 0 : 1] real alpha; // scale of of noise GP vector[fixed ? 0 : M] eta; // unconstrained noise // Rt vector[estimate_r] log_R; // baseline reproduction number estimate (log) - real initial_infections[estimate_r] ; // seed infections - real initial_growth[estimate_r && seeding_time > 1 ? 1 : 0]; // seed growth rate - real bp_sd[bp_n > 0 ? 1 : 0]; // standard deviation of breakpoint effect - real bp_effects[bp_n]; // Rt breakpoint effects + array[estimate_r] real initial_infections ; // seed infections + array[estimate_r && seeding_time > 1 ? 1 : 0] real initial_growth; // seed growth rate + array[bp_n > 0 ? 1 : 0] real bp_sd; // standard deviation of breakpoint effect + array[bp_n] real bp_effects; // Rt breakpoint effects // observation model - real delay_mean[delay_n_p]; // mean of delays - real delay_sd[delay_n_p]; // sd of delays + array[delay_n_p] real delay_mean; // mean of delays + array[delay_n_p] real delay_sd; // sd of delays simplex[week_effect] day_of_week_simplex;// day of week reporting effect - real frac_obs[obs_scale]; // fraction of cases that are ultimately observed - real rep_phi[model_type]; // overdispersion of the reporting process + array[obs_scale] real frac_obs; // fraction of cases that are ultimately observed + array[model_type] real rep_phi; // overdispersion of the reporting process } transformed parameters { @@ -154,9 +154,9 @@ model { } generated quantities { - int imputed_reports[ot_h]; + array[ot_h] int imputed_reports; vector[estimate_r > 0 ? 0: ot_h] gen_R; - real r[ot_h]; + array[ot_h] real r; real gt_mean; real gt_var; vector[return_likelihood ? ot : 0] log_lik; @@ -167,9 +167,9 @@ generated quantities { r = R_to_growth(R, gt_mean, gt_var); } else { // sample generation time - real delay_mean_sample[delay_n_p] = + array[delay_n_p] real delay_mean_sample = normal_rng(delay_mean_mean, delay_mean_sd); - real delay_sd_sample[delay_n_p] = + array[delay_n_p] real delay_sd_sample = normal_rng(delay_sd_mean, delay_sd_sd); vector[delay_type_max[gt_id]] sampled_gt_rev_pmf = get_delay_rev_pmf( gt_id, delay_type_max[gt_id], delay_types_p, delay_types_id, diff --git a/inst/stan/estimate_secondary.stan b/inst/stan/estimate_secondary.stan index 0b4770a3e..51468f8d9 100644 --- a/inst/stan/estimate_secondary.stan +++ b/inst/stan/estimate_secondary.stan @@ -8,7 +8,7 @@ functions { data { int t; // time of observations - int obs[t]; // observed secondary data + array[t] int obs; // observed secondary data vector[t] primary; // observed primary data int burn_in; // time period to not use for fitting #include data/secondary.stan @@ -17,7 +17,7 @@ data { } transformed data{ - int delay_type_max[delay_types] = get_delay_type_max( + array[delay_types] int delay_type_max = get_delay_type_max( delay_types, delay_types_p, delay_types_id, delay_types_groups, delay_max, delay_np_pmf_groups ); @@ -25,11 +25,11 @@ transformed data{ parameters{ // observation model - real delay_mean[delay_n_p]; - real delay_sd[delay_n_p]; // sd of delays + array[delay_n_p] real delay_mean; + array[delay_n_p] real delay_sd; // sd of delays simplex[week_effect] day_of_week_simplex; // day of week reporting effect - real frac_obs[obs_scale]; // fraction of cases that are ultimately observed - real rep_phi[model_type]; // overdispersion of the reporting process + array[obs_scale] real frac_obs; // fraction of cases that are ultimately observed + array[model_type] real rep_phi; // overdispersion of the reporting process } transformed parameters { @@ -89,7 +89,7 @@ model { } generated quantities { - int sim_secondary[t - burn_in]; + array[t - burn_in] int sim_secondary; vector[return_likelihood > 1 ? t - burn_in : 0] log_lik; // simulate secondary reports sim_secondary = report_rng(secondary[(burn_in + 1):t], rep_phi, model_type); diff --git a/inst/stan/estimate_truncation.stan b/inst/stan/estimate_truncation.stan index 3eaf801a9..0df4006f4 100644 --- a/inst/stan/estimate_truncation.stan +++ b/inst/stan/estimate_truncation.stan @@ -7,16 +7,16 @@ functions { data { int t; int obs_sets; - int obs[t, obs_sets]; - int obs_dist[obs_sets]; + array[t, obs_sets] int obs; + array[obs_sets] int obs_dist; #include data/delays.stan } transformed data{ int trunc_id = 1; - int end_t[obs_sets]; - int start_t[obs_sets]; + array[obs_sets] int end_t; + array[obs_sets] int start_t; - int delay_type_max[delay_types]; + array[delay_types] int delay_type_max; delay_type_max = get_delay_type_max( delay_types, delay_types_p, delay_types_id, delay_types_groups, delay_max, delay_np_pmf_groups @@ -28,8 +28,8 @@ transformed data{ } } parameters { - real delay_mean[delay_n_p]; - real delay_sd[delay_n_p]; // sd of delays + array[delay_n_p] real delay_mean; + array[delay_n_p] real delay_sd; // sd of delays real phi; real sigma; } diff --git a/inst/stan/functions/delays.stan b/inst/stan/functions/delays.stan index 69a88ba17..01cbd98b5 100644 --- a/inst/stan/functions/delays.stan +++ b/inst/stan/functions/delays.stan @@ -1,8 +1,8 @@ -int[] get_delay_type_max( - int delay_types, int[] delay_types_p, int[] delay_types_id, - int[] delay_types_groups, int[] delay_max, int[] delay_np_pmf_groups +array[] int get_delay_type_max( + int delay_types, array[] int delay_types_p, array[] int delay_types_id, + array[] int delay_types_groups, array[] int delay_max, array[] int delay_np_pmf_groups ) { - int ret[delay_types]; + array[delay_types] int ret; for (i in 1:delay_types) { ret[i] = 1; for (j in delay_types_groups[i]:(delay_types_groups[i + 1] - 1)) { @@ -18,10 +18,10 @@ int[] get_delay_type_max( } vector get_delay_rev_pmf( - int delay_id, int len, int[] delay_types_p, int[] delay_types_id, - int[] delay_types_groups, int[] delay_max, - vector delay_np_pmf, int[] delay_np_pmf_groups, - real[] delay_mean, real[] delay_sigma, int[] delay_dist, + int delay_id, int len, array[] int delay_types_p, array[] int delay_types_id, + array[] int delay_types_groups, array[] int delay_max, + vector delay_np_pmf, array[] int delay_np_pmf_groups, + array[] real delay_mean, array[] real delay_sigma, array[] int delay_dist, int left_truncate, int reverse_pmf, int cumulative ) { // loop over delays @@ -75,9 +75,9 @@ vector get_delay_rev_pmf( } -void delays_lp(real[] delay_mean, real[] delay_mean_mean, real[] delay_mean_sd, - real[] delay_sd, real[] delay_sd_mean, real[] delay_sd_sd, - int[] delay_dist, int[] weight) { +void delays_lp(array[] real delay_mean, array[] real delay_mean_mean, array[] real delay_mean_sd, + array[] real delay_sd, array[] real delay_sd_mean, array[] real delay_sd_sd, + array[] int delay_dist, array[] int weight) { int mean_delays = num_elements(delay_mean); int sd_delays = num_elements(delay_sd); if (mean_delays) { diff --git a/inst/stan/functions/generated_quantities.stan b/inst/stan/functions/generated_quantities.stan index e2bc1d064..003fbee52 100644 --- a/inst/stan/functions/generated_quantities.stan +++ b/inst/stan/functions/generated_quantities.stan @@ -29,9 +29,9 @@ vector calculate_Rt(vector infections, int seeding_time, return(sR); } // Convert an estimate of Rt to growth -real[] R_to_growth(vector R, real gt_mean, real gt_var) { +array[] real R_to_growth(vector R, real gt_mean, real gt_var) { int t = num_elements(R); - real r[t]; + array[t] real r; if (gt_var > 0) { real k = gt_var / pow(gt_mean, 2); for (s in 1:t) { diff --git a/inst/stan/functions/infections.stan b/inst/stan/functions/infections.stan index c02d9e22e..db1de31a7 100644 --- a/inst/stan/functions/infections.stan +++ b/inst/stan/functions/infections.stan @@ -19,7 +19,7 @@ real update_infectiousness(vector infections, vector gt_rev_pmf, } // generate infections by using Rt = Rt-1 * sum(reversed generation time pmf * infections) vector generate_infections(vector oR, int uot, vector gt_rev_pmf, - real[] initial_infections, real[] initial_growth, + array[] real initial_infections, array[] real initial_growth, int pop, int ht) { // time indices and storage int ot = num_elements(oR); diff --git a/inst/stan/functions/observation_model.stan b/inst/stan/functions/observation_model.stan index 30d2b013c..3968d3aa8 100644 --- a/inst/stan/functions/observation_model.stan +++ b/inst/stan/functions/observation_model.stan @@ -1,5 +1,5 @@ // apply day of the week effect -vector day_of_week_effect(vector reports, int[] day_of_week, vector effect) { +vector day_of_week_effect(vector reports, array[] int day_of_week, vector effect) { int t = num_elements(reports); int wl = num_elements(effect); // scale day of week effect @@ -35,9 +35,9 @@ vector truncate(vector reports, vector trunc_rev_cmf, int reconstruct) { return(trunc_reports); } // Truncation distribution priors -void truncation_lp(real[] truncation_mean, real[] truncation_sd, - real[] trunc_mean_mean, real[] trunc_mean_sd, - real[] trunc_sd_mean, real[] trunc_sd_sd) { +void truncation_lp(array[] real truncation_mean, array[] real truncation_sd, + array[] real trunc_mean_mean, array[] real trunc_mean_sd, + array[] real trunc_sd_mean, array[] real trunc_sd_sd) { int truncation = num_elements(truncation_mean); if (truncation) { if (trunc_mean_sd[1] > 0) { @@ -51,8 +51,8 @@ void truncation_lp(real[] truncation_mean, real[] truncation_sd, } } // update log density for reported cases -void report_lp(int[] cases, vector reports, - real[] rep_phi, real phi_mean, real phi_sd, +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) @@ -72,8 +72,8 @@ void report_lp(int[] cases, vector reports, } } // 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) { +vector report_log_lik(array[] int cases, vector reports, + array[] real rep_phi, int model_type, real weight) { int t = num_elements(reports); vector[t] log_lik; @@ -91,9 +91,9 @@ vector report_log_lik(int[] cases, vector reports, return(log_lik); } // sample reported cases from the observation model -int[] report_rng(vector reports, real[] rep_phi, int model_type) { +array[] int report_rng(vector reports, array[] real rep_phi, int model_type) { int t = num_elements(reports); - int sampled_reports[t]; + array[t] int sampled_reports; real sqrt_phi = 1e5; if (model_type) { sqrt_phi = 1 / sqrt(rep_phi[model_type]); diff --git a/inst/stan/functions/pmfs.stan b/inst/stan/functions/pmfs.stan index 9b077b28e..adc27d53e 100644 --- a/inst/stan/functions/pmfs.stan +++ b/inst/stan/functions/pmfs.stan @@ -5,33 +5,33 @@ // @author Sam Abbott // @author Adrian Lison vector discretised_pmf(real mu, real sigma, int n, int dist) { - vector[n] pmf; + vector[n] lpmf; if (sigma > 0) { - vector[n] upper_cdf; + vector[n] upper_lcdf; if (dist == 0) { for (i in 1:n) { - upper_cdf[i] = lognormal_cdf(i, mu, sigma); + upper_lcdf[i] = lognormal_lcdf(i | mu, sigma); } } else if (dist == 1) { real alpha = mu^2 / sigma^2; real beta = mu / sigma^2; for (i in 1:n) { - upper_cdf[i] = gamma_cdf(i, alpha, beta); + upper_lcdf[i] = gamma_lcdf(i | alpha, beta); } } else { reject("Unknown distribution function provided."); } // discretise - pmf[1] = upper_cdf[1]; - pmf[2:n] = upper_cdf[2:n] - upper_cdf[1:(n-1)]; + lpmf[1] = upper_lcdf[1]; + lpmf[2:n] = log_diff_exp(upper_lcdf[2:n], upper_lcdf[1:(n-1)]); // normalize - pmf = pmf / upper_cdf[n]; + lpmf = lpmf - upper_lcdf[n]; } else { // delta function - pmf = rep_vector(0, n); - pmf[n] = 1; + lpmf = rep_vector(negative_infinity(), n); + lpmf[n] = 0; } - return(pmf); + return(exp(lpmf)); } // reverse a mf diff --git a/inst/stan/functions/rt.stan b/inst/stan/functions/rt.stan index 95c0f48fc..a2ba59e25 100644 --- a/inst/stan/functions/rt.stan +++ b/inst/stan/functions/rt.stan @@ -1,6 +1,6 @@ // update a vector of Rts -vector update_Rt(int t, real log_R, vector noise, int[] bps, - real[] bp_effects, int stationary) { +vector update_Rt(int t, real log_R, vector noise, array[] int bps, + array[] real bp_effects, int stationary) { // define control parameters int bp_n = num_elements(bp_effects); int bp_c = 0; @@ -38,8 +38,8 @@ vector update_Rt(int t, real log_R, vector noise, int[] bps, return(R); } // Rt priors -void rt_lp(vector log_R, real[] initial_infections, real[] initial_growth, - real[] bp_effects, real[] bp_sd, int bp_n, int seeding_time, +void rt_lp(vector log_R, array[] real initial_infections, array[] real initial_growth, + array[] real bp_effects, array[] real bp_sd, int bp_n, int seeding_time, real r_logmean, real r_logsd, real prior_infections, real prior_growth) { // prior on R diff --git a/inst/stan/functions/secondary.stan b/inst/stan/functions/secondary.stan index 5c89c4295..5a96a7cf9 100644 --- a/inst/stan/functions/secondary.stan +++ b/inst/stan/functions/secondary.stan @@ -1,5 +1,5 @@ // Calculate secondary reports condition only on primary reports -vector calculate_secondary(vector reports, int[] obs, real[] frac_obs, +vector calculate_secondary(vector reports, array[] int obs, array[] real frac_obs, vector delay_rev_pmf, int cumulative, int historic, int primary_hist_additive, int current, int primary_current_additive, int predict) { diff --git a/inst/stan/simulate_infections.stan b/inst/stan/simulate_infections.stan index a6befed21..85b4452c8 100644 --- a/inst/stan/simulate_infections.stan +++ b/inst/stan/simulate_infections.stan @@ -24,7 +24,7 @@ data { } transformed data { - int delay_type_max[delay_types] = get_delay_type_max( + array[delay_types] int delay_type_max = get_delay_type_max( delay_types, delay_types_p, delay_types_id, delay_types_groups, delay_max, delay_np_pmf_groups ); @@ -34,8 +34,8 @@ generated quantities { // generated quantities matrix[n, t] infections; //latent infections matrix[n, t - seeding_time] reports; // observed cases - int imputed_reports[n, t - seeding_time]; - real r[n, t - seeding_time]; + array[n, t - seeding_time] int imputed_reports; + array[n, t - seeding_time] real r; for (i in 1:n) { // generate infections from Rt trace vector[delay_type_max[gt_id]] gt_rev_pmf; diff --git a/inst/stan/simulate_secondary.stan b/inst/stan/simulate_secondary.stan index c27de4d3f..0370b705f 100644 --- a/inst/stan/simulate_secondary.stan +++ b/inst/stan/simulate_secondary.stan @@ -13,7 +13,7 @@ data { int h; // forecast horizon int all_dates; // should all dates have simulations returned // secondary model specific data - int obs[t - h]; // observed secondary data + array[t - h] int obs; // observed secondary data matrix[n, t] primary; // observed primary data #include data/secondary.stan // delay from infection to report @@ -23,14 +23,14 @@ data { } transformed data { - int delay_type_max[delay_types] = get_delay_type_max( + array[delay_types] int delay_type_max = get_delay_type_max( delay_types, delay_types_p, delay_types_id, delay_types_groups, delay_max, delay_np_pmf_groups ); } generated quantities { - int sim_secondary[n, all_dates ? t : h]; + array[n, all_dates ? t : h] int sim_secondary; for (i in 1:n) { vector[t] secondary; vector[delay_type_max[delay_id]] delay_rev_pmf = get_delay_rev_pmf(