7 return floor(R_unif_index(
n));
11 int n,
int from,
int to) {
15 if (!ISNA(color[i]) && nearbyint(color[i]) == from)
n--;
19 assert(nearbyint(color[i]) == from);
23 #define Beta (__p[__parindex[0]])
24 #define sigma (__p[__parindex[1]])
25 #define gamma (__p[__parindex[2]])
26 #define psi (__p[__parindex[3]])
27 #define omega (__p[__parindex[4]])
28 #define S0 (__p[__parindex[5]])
29 #define E0 (__p[__parindex[6]])
30 #define I0 (__p[__parindex[7]])
31 #define R0 (__p[__parindex[8]])
32 #define N (__p[__parindex[9]])
33 #define S (__x[__stateindex[0]])
34 #define E (__x[__stateindex[1]])
35 #define I (__x[__stateindex[2]])
36 #define R (__x[__stateindex[3]])
37 #define ll (__x[__stateindex[4]])
38 #define node (__x[__stateindex[5]])
39 #define ellE (__x[__stateindex[6]])
40 #define ellI (__x[__stateindex[7]])
41 #define COLOR (__x[__stateindex[8]])
44 event_rates(__x,__p,t, \
45 __stateindex,__parindex,__covindex, \
46 __covars,rate,logpi,&penalty) \
53 const int *__stateindex,
54 const int *__parindex,
55 const int *__covindex,
56 const double *__covars,
61 double event_rate = 0;
67 pi = (
I > 0) ? 1-
ellI/
I : 0;
69 event_rate += (*rate = alpha*pi); rate++;
70 *logpi = log(pi); logpi++;
73 event_rate += (*rate = alpha*pi); rate++;
74 *logpi = log(pi)-log(
ellI); logpi++;
76 event_rate += (*rate = alpha*pi); rate++;
77 *logpi = log(pi)-log(
ellI); logpi++;
83 event_rate += (*rate = alpha*(1-pi)); rate++;
84 *logpi = log(1-pi); logpi++;
86 event_rate += (*rate = alpha*pi); rate++;
87 *logpi = log(pi)-log(
ellE); logpi++;
92 event_rate += (*rate = alpha); rate++;
100 event_rate += (*rate =
omega*
R); rate++;
104 assert(R_FINITE(event_rate));
119 const int *__stateindex,
120 const int *__parindex,
121 const int *__covindex,
122 const double *__covars
125 S = nearbyint(
S0*adj);
126 E = nearbyint(
E0*adj);
127 I = nearbyint(
I0*adj);
128 R = nearbyint(
R0*adj);
143 const int *__stateindex,
144 const int *__parindex,
145 const int *__covindex,
146 const double *__covars,
150 double tstep = 0.0, tmax = t + dt;
151 double *color = &
COLOR;
159 int parent = (int) nearbyint(
node);
164 assert(parent<=nnode);
167 int parlin = lineage[parent];
168 int parcol = color[parlin];
169 assert(parlin >= 0 && parlin <
nsample);
172 switch (nodetype[parent]) {
178 assert(sat[parent]==1);
179 int c = child[index[parent]];
180 assert(lineage[parent]==lineage[c]);
183 if (unif_rand() < x) {
184 color[lineage[c]] = 0;
188 color[lineage[c]] = 1;
195 if (unif_rand() < 0.5) {
196 color[lineage[c]] = 0;
200 color[lineage[c]] = 1;
216 if (sat[parent] == 1) {
217 int c = child[index[parent]];
218 color[lineage[c]] = 1;
220 }
else if (sat[parent] == 0) {
227 color[parlin] = R_NaReal;
239 assert(sat[parent]==2);
240 ll += (
S > 0 &&
I > 0) ? log(
Beta*
S/
N/(
E+1)) : R_NegInf;
244 int c1 = child[index[parent]];
245 int c2 = child[index[parent]+1];
247 assert(lineage[c1] != lineage[c2]);
248 assert(lineage[c1] != parlin || lineage[c2] != parlin);
249 assert(lineage[c1] == parlin || lineage[c2] == parlin);
250 if (unif_rand() < 0.5) {
251 color[lineage[c1]] = 0;
252 color[lineage[c2]] = 1;
254 color[lineage[c1]] = 1;
255 color[lineage[c2]] = 0;
267 double event_rate = 0;
271 tstep = exp_rand()/event_rate;
273 while (t + tstep < tmax) {
275 ll -= penalty*tstep + logpi[event];
332 tstep = exp_rand()/event_rate;
341 # define lik (__lik[0])
351 const int *__obsindex,
352 const int *__stateindex,
353 const int *__parindex,
354 const int *__covindex,
355 const double *__covars,
359 lik = (give_log) ?
ll : exp(
ll);
get_userdata_int_t * get_userdata_int
static int rcateg(double erate, double *rate, int nrate)
static void change_color(double *color, int nsample, int n, int from, int to)
static double event_rates(double *__x, const double *__p, double t, const int *__stateindex, const int *__parindex, const int *__covindex, const double *__covars, double *rate, double *logpi, double *penalty)
static int random_choice(double n)
void seirs_gill(double *__x, const double *__p, const int *__stateindex, const int *__parindex, const int *__covindex, const double *__covars, double t, double dt)
void seirs_dmeas(double *__lik, const double *__y, const double *__x, const double *__p, int give_log, const int *__obsindex, const int *__stateindex, const int *__parindex, const int *__covindex, const double *__covars, double t)
Measurement model likelihood (dmeasure).
void seirs_rinit(double *__x, const double *__p, double t0, const int *__stateindex, const int *__parindex, const int *__covindex, const double *__covars)