Skip to content

Instantly share code, notes, and snippets.

@euxoa
Last active October 30, 2020 12:19
Show Gist options
  • Save euxoa/b3dcbd60918ba2b43caa8030754f44f4 to your computer and use it in GitHub Desktop.
Save euxoa/b3dcbd60918ba2b43caa8030754f44f4 to your computer and use it in GitHub Desktop.
data {
int<lower=0> N; // rows
int<lower=1> D; // delays, first is 0, last D-1
int<lower=0> T; // days, first is 1
int<lower=0> counts[N];
int<lower=1> t[N];
int<lower=0, upper=D-1> dt0[N];
int<lower=0, upper=D-1> dt1[N];
real<lower=0> dir_prior;
}
transformed data {
vector<lower=0>[D] theta;
for (d in 1:D) theta[d] = dir_prior;
}
parameters {
simplex[D] beta;
real<lower=0> sd_change;
real bl; // baseline (intercept) for latent intensity of the epidemic
real k; // intercept for the growth of the epidemic
vector[T] u_nova; // increments (innovations) of growth over time, standardized
real<lower=0> sd_overdisp[4];
vector[6] weekdays_6;
// matrix[T, D] u_dispersion; // deviations due to overdispersion, log-scale, standardized
matrix[T, D] log_lambda;
}
transformed parameters {
vector[T] log_growth;
vector[T] alpha;
vector[7] weekdays = append_row(weekdays_6, -sum(weekdays_6));
log_growth = k + cumulative_sum(sd_change * u_nova);
alpha = bl + cumulative_sum(log_growth);
}
model {
//matrix[T, D] log_lambda;
for (i in 1:T)
for (d in 1:D)
log_lambda[i, d] ~ normal(alpha[i] + weekdays[1 + i%7] + log(beta[d]),
sd_overdisp[d>4 ? 4 : d]);
for (i in 1:N) {
int d0 = dt0[i] + 1;
int d1 = dt1[i] + 1;
if (d0==d1)
counts[i] ~ poisson_log(log_lambda[t[i], d0]); // redundant but more efficient maybe
else
counts[i] ~ poisson_log(log_sum_exp(log_lambda[t[i], d0:d1]));
}
// for (i in 1:T) u_dispersion[i] ~ normal(0, 1);
u_nova ~ normal(0, 1); // student_t(4, 0, 1);
sd_change ~ normal(0, 0.1);
bl ~ normal(0, 10);
k ~ normal(0, .2);
beta ~ dirichlet(theta);
sd_overdisp ~ lognormal(0, .5);
weekdays ~ normal(0, 1);
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment