epinowcast
Loading...
Searching...
No Matches
delay_lpmf.stan
Go to the documentation of this file.
1
28real delay_snap_lpmf(array[] int dummy, int start, int end, array[] int obs,
29 array[] int sl, array[] int csl,
30 array[] int nsl, array[] int cnsl, array[] int obs_lookup,
31 array[] vector imp_obs, array[] int sg,
32 array[] int st, array[,] int rdlurd,
33 vector srdlh, matrix refp_lh, array[] int dpmfs,
34 int ref_p, int rep_h, int ref_as_p, array[] real phi,
35 int model_obs, vector refnp_lh, int ref_np,
36 array[] int sdmax, array[] int csdmax) {
37 real tar = 0;
38 // Where am I in the observed data?
39 array[3] int nc = filt_obs_indexes(start, end, cnsl, nsl);
40 // Where am I in the observed data filling in gaps?
41 array[3] int n = filt_obs_indexes(start, end, csl, sl);
42 if (nc[3]) {
43 // Filter observed data and observed data lookup
44 array[nc[3]] int filt_obs = segment(obs, nc[1], nc[3]);
45 array[nc[3]] int filt_obs_lookup = segment(obs_lookup, nc[1], nc[3]);
46 array[nc[3]] int filt_obs_lookup_local;
47 for (i in 1:nc[3]) {
48 filt_obs_lookup_local[i] = filt_obs_lookup[i] - n[1] + 1;
49 }
50
51 // What is going to be used for storage
52 vector[n[3]] log_exp_obs;
53
54 // combine expected final obs and time effects to get expected obs
55 log_exp_obs = expected_obs_from_snaps(
56 start, end, imp_obs, rdlurd, srdlh, refp_lh, dpmfs, ref_p, rep_h, ref_as_p, sl, csl, sg, st, n[3], refnp_lh, ref_np, sdmax, csdmax
57 );
58
59 // observation error model (across all reference dates and groups)
60 profile("model_likelihood_neg_binomial") {
61 tar = obs_lpmf(
62 filt_obs | log_exp_obs[filt_obs_lookup_local], phi, model_obs
63 );
64 }
65 }
66 return(tar);
67}
68
122real delay_group_lpmf(array[] int groups, int start, int end, array[] int obs,
123 array[] int sl, array[] int csl, array[] int nsl,
124 array[] int cnsl, array[] int obs_lookup,
125 array[] vector imp_obs, int t, array[] int sg,
126 array[,] int ts, array[] int st,
127 array[,] int rdlurd, vector srdlh, matrix refp_lh,
128 array[] int dpmfs, int ref_p, int rep_h, int ref_as_p,
129 array[] real phi, int model_obs, int model_miss,
130 int miss_obs, array[] int missing_reference,
131 array[,] int obs_by_report, vector miss_ref_lprop,
132 array[] int sdmax, array[] int csdmax,
133 array[] int miss_st, array[] int miss_cst,
134 vector refnp_lh, int ref_np) {
135 // Where am I?
136 real tar = 0;
137 int i_start = ts[1, start];
138 int i_end = ts[t, end];
139 // Where am I in the observed data?
140 array[3] int nc = filt_obs_indexes(i_start, i_end, cnsl, nsl);
141 // Where am I in the observed data filling in gaps?
142 array[3] int n = filt_obs_indexes(i_start, i_end, csl, sl);
143
144 // What is going to be used for storage
145 vector[n[3]] log_exp_obs;
146 vector[model_miss ? miss_obs : 0] log_exp_obs_miss;
147
148 // Combine expected final obs and time effects to get expected obs
149 // If missing reference module in place calculate all expected obs vs
150 // just those observed and allocate if missing or not.
151 if (model_miss) {
152 array[3] int f = filt_obs_indexes(i_start, i_end, csdmax, sdmax);
153 array[3] int l = filt_obs_indexes(start, end, miss_cst, miss_st);
154 vector[f[3]] log_exp_all;
155
156 // Calculate all expected observations
157 log_exp_all = expected_obs_from_snaps(
158 i_start, i_end, imp_obs, rdlurd, srdlh, refp_lh, dpmfs, ref_p, rep_h, ref_as_p, sdmax, csdmax, sg, st, f[3], refnp_lh, ref_np, sdmax, csdmax
159 );
160
161 // Allocate to just those actually observed
162 log_exp_obs = allocate_observed_obs(
163 i_start, i_end, log_exp_all, sl, csl, sdmax, csdmax
164 );
166 i_start, i_end, log_exp_obs, sl, csl, log1m_exp(miss_ref_lprop)
167 );
168
169 // Allocate expected cases by reporting time
170 if (miss_obs) {
172 i_start, i_end, log_exp_all, sdmax, csdmax, miss_ref_lprop
173 );
174 log_exp_obs_miss = log_expected_by_report(
175 log_exp_all, obs_by_report[l[1]:l[2]]
176 );
177 }
178 }else{
179 log_exp_obs = expected_obs_from_snaps(
180 i_start, i_end, imp_obs, rdlurd, srdlh, refp_lh, dpmfs, ref_p, rep_h, ref_as_p, sl, csl, sg, st, n[3], refnp_lh, ref_np, sdmax, csdmax
181 );
182 }
183 // Observation error model (across all reference dates and groups)
184 profile("model_likelihood_neg_binomial") {
185 if (nc[3]) {
186 // Filter observed data and observed data lookup
187 array[nc[3]] int filt_obs = segment(obs, nc[1], nc[3]);
188 array[nc[3]] int filt_obs_lookup = segment(obs_lookup, nc[1], nc[3]);
189 array[nc[3]] int filt_obs_lookup_local;
190 for (i in 1:nc[3]) {
191 filt_obs_lookup_local[i] = filt_obs_lookup[i] - n[1] + 1;
192 }
193
194 tar = obs_lpmf(
195 filt_obs | log_exp_obs[filt_obs_lookup_local], phi, model_obs
196 );
197 }
198 if (model_miss && miss_obs) {
199 array[3] int l = filt_obs_indexes(start, end, miss_cst, miss_st);
200 array[l[3]] int filt_miss_ref = segment(missing_reference, l[1], l[3]);
201 tar += obs_lpmf(filt_miss_ref | log_exp_obs_miss, phi, model_obs);
202 }
203 }
204 return(tar);
205}
vector allocate_observed_obs(int start, int end, vector obs, array[] int sl, array[] int csl, array[] int sdmax, array[] int csdmax)
vector apply_missing_reference_effects(int start, int end, vector obs, array[] int sl, array[] int csl, vector miss_ref_lprop)
real delay_group_lpmf(array[] int groups, int start, int end, array[] int obs, array[] int sl, array[] int csl, array[] int nsl, array[] int cnsl, array[] int obs_lookup, array[] vector imp_obs, int t, array[] int sg, array[,] int ts, array[] int st, array[,] int rdlurd, vector srdlh, matrix refp_lh, array[] int dpmfs, int ref_p, int rep_h, int ref_as_p, array[] real phi, int model_obs, int model_miss, int miss_obs, array[] int missing_reference, array[,] int obs_by_report, vector miss_ref_lprop, array[] int sdmax, array[] int csdmax, array[] int miss_st, array[] int miss_cst, vector refnp_lh, int ref_np)
real delay_snap_lpmf(array[] int dummy, int start, int end, array[] int obs, array[] int sl, array[] int csl, array[] int nsl, array[] int cnsl, array[] int obs_lookup, array[] vector imp_obs, array[] int sg, array[] int st, array[,] int rdlurd, vector srdlh, matrix refp_lh, array[] int dpmfs, int ref_p, int rep_h, int ref_as_p, array[] real phi, int model_obs, vector refnp_lh, int ref_np, array[] int sdmax, array[] int csdmax)
vector expected_obs_from_snaps(int start, int end, array[] vector imp_obs, array[,] int rdlurd, vector srdlh, matrix refp_lh, array[] int dpmfs, int ref_p, int rep_h, int ref_as_p, array[] int sl, array[] int csl, array[] int sg, array[] int st, int n, vector refnp_lh, int ref_np, array[] int sdmax, array[] int csdmax)
array[] int filt_obs_indexes(int start, int end, array[] int csl, array[] int sl)
vector log_expected_by_report(vector log_exp, array[,] int obs_by_report)
real obs_lpmf(array[] int obs, vector log_exp_obs, array[] real phi, int model_obs)
Definition obs_lpmf.stan:27