Skip to content

Commit 7e47dfd

Browse files
committed
Guard growth-rate functions against insufficient data explicitly
Co-authored-by: sbfnk <sebastian.funk@lshtm.ac.uk>
1 parent 27bd7dc commit 7e47dfd

1 file changed

Lines changed: 18 additions & 8 deletions

File tree

inst/stan/functions/generated_quantities.stan

Lines changed: 18 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -92,6 +92,11 @@ vector calculate_growth_infections(vector infections, int seeding_time) {
9292
int t = num_elements(infections);
9393
int ot = t - seeding_time;
9494
int start = 1 + seeding_time;
95+
if (ot <= 1) {
96+
reject(
97+
"seeding_time must be >1 time step shorter than the infections vector."
98+
);
99+
}
95100
vector[t] log_inf = log(infections);
96101
vector[ot - 1] growth = log_inf[(1+start):t] - log_inf[start:(t - 1)];
97102
return(growth);
@@ -118,21 +123,26 @@ vector calculate_growth_infness(vector infections, int seeding_time,
118123
int t = num_elements(infections);
119124
int ot = t - seeding_time;
120125
int start = 1 + seeding_time;
121-
if (ot <= 1) {
122-
reject("seeding_time must >1 time step shorter than infections vector.");
123-
}
126+
// mean generation time, will always be >= 1
127+
int gt_length = num_elements(gt_rev_pmf);
128+
int mean_gen = to_int(round( // round weighted mean to nearest int
129+
dot_product(reverse(linspaced_vector(gt_length, 1, gt_length)), gt_rev_pmf)
130+
));
131+
// a shifted growth rate needs more modeled points than the mean generation
132+
// time, otherwise the shift below would index out of range
133+
if (ot <= mean_gen) {
134+
reject(
135+
"infections vector must be longer than seeding_time by more than the ",
136+
"mean generation time."
137+
);
138+
}
124139
// infectiousness
125140
vector[ot] infness_log = rep_vector(1e-5, ot);
126141
for (s in 1:ot) {
127142
infness_log[s] += log(update_infectiousness(
128143
infections, gt_rev_pmf, seeding_time, s
129144
));
130145
}
131-
// mean generation time, will always be >= 1
132-
int gt_length = num_elements(gt_rev_pmf);
133-
int mean_gen = to_int(round( // round weighted mean to nearest int
134-
dot_product(reverse(linspaced_vector(gt_length, 1, gt_length)), gt_rev_pmf)
135-
));
136146
// growth rate
137147
vector[ot - 1] growth = infness_log[2:ot] - infness_log[1:(ot - 1)];
138148
// shift by mean_gen (most recent growth rates remain undefined)

0 commit comments

Comments
 (0)