-
-
Notifications
You must be signed in to change notification settings - Fork 479
/
sir.stan
58 lines (53 loc) · 1.47 KB
/
sir.stan
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
// Simple SIR model inspired by the presentation in
// http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3380087/pdf/nihms372789.pdf
functions {
vector simple_SIR(real t, vector y, real beta,
// water contact rate
real kappa, // C_{50}
real gamma, // recovery rate
real xi, // bacteria production rate
real delta) {
// bacteria removal rate
vector[4] dydt;
dydt[1] = -beta * y[4] / (y[4] + kappa) * y[1];
dydt[2] = beta * y[4] / (y[4] + kappa) * y[1] - gamma * y[2];
dydt[3] = gamma * y[2];
dydt[4] = xi * y[2] - delta * y[4];
return dydt;
}
}
data {
int<lower=0> N_t;
array[N_t] real t;
vector[4] y0;
array[N_t] int stoi;
array[N_t] real B;
}
transformed data {
real t0 = 0;
real<lower=0> kappa = 1000000;
}
parameters {
real<lower=0> beta;
real<lower=0> gamma;
real<lower=0> xi;
real<lower=0> delta;
}
transformed parameters {
array[N_t] vector<lower=0>[4] y = ode_rk45_tol(simple_SIR, y0, t0, t, 1e-6,
1e-6, 1000, beta, kappa,
gamma, xi, delta);
}
model {
vector[N_t] y_diff;
y_diff[1] = y0[1] - y[1, 1];
for (n in 2 : N_t) {
y_diff[n] = y[n - 1, 1] - y[n, 1];
}
beta ~ cauchy(0, 2.5);
gamma ~ cauchy(0, 1);
xi ~ cauchy(0, 25);
delta ~ cauchy(0, 1);
stoi ~ poisson(y_diff);
B ~ lognormal(log(y[ : , 4]), 0.15);
}