pomp: Constructor of the basic POMP object

In pomp: Statistical Inference for Partially Observed Markov Processes

Description

This function constructs a pomp object, encoding a partially-observed Markov process model together with a uni- or multivariate time series. One implements the model by specifying its components, each of which can be written as R functions or, for much greater computational efficiency, using C code. The preferred way to specify most components (as detailed below) is through the use of Csnippets, snippets of C that are compiled and linked into a running R session.

Usage

## S4 method for signature 'data.frame'
pomp(data, times, t0, ..., rprocess, dprocess, rmeasure, dmeasure,
       measurement.model,
       skeleton, skeleton.type = c("map","vectorfield"), skelmap.delta.t = 1,
       initializer, rprior, dprior, params, covar, tcovar,
       obsnames, statenames, paramnames, covarnames, zeronames,
       PACKAGE, fromEstimationScale, toEstimationScale, globals)
## S4 method for signature 'numeric'
pomp(data, times, t0, ..., rprocess, dprocess, rmeasure, dmeasure,
       measurement.model,
       skeleton, skeleton.type = c("map","vectorfield"), skelmap.delta.t = 1,
       initializer, rprior, dprior, params, covar, tcovar,
       obsnames, statenames, paramnames, covarnames, zeronames,
       PACKAGE, fromEstimationScale, toEstimationScale, globals)
## S4 method for signature 'matrix'
pomp(data, times, t0, ..., rprocess, dprocess, rmeasure, dmeasure,
       measurement.model,
       skeleton, skeleton.type = c("map","vectorfield"), skelmap.delta.t = 1,
       initializer, rprior, dprior, params, covar, tcovar,
       obsnames, statenames, paramnames, covarnames, zeronames,
       PACKAGE, fromEstimationScale, toEstimationScale, globals)
## S4 method for signature 'pomp'
pomp(data, times, t0, ..., rprocess, dprocess, rmeasure, dmeasure,
       measurement.model, skeleton, skeleton.type, skelmap.delta.t,
       initializer, rprior, dprior, params, covar, tcovar,
       obsnames, statenames, paramnames, covarnames, zeronames,
       PACKAGE, fromEstimationScale, toEstimationScale, globals)

Arguments

data, times	The time series data and times at which observations are made. data can be specified as a vector, a matrix, a data-frame. Alternatively, a pomp object can be supplied in the data argument. If data is a numeric vector, times must be a numeric vector of the same length. If data is a matrix, it should have dimensions nobs x ntimes, where nobs is the number of observed variables and ntimes is the number of times at which observations were made (i.e., each column is a distinct observation of the nobs variables). In this case, times must be given as a numeric vector (of length ntimes). If data is a data-frame, times must name the column of observation times. Note that, in this case, data will be internally coerced to an array with storage-mode ‘double’. times must be numeric and strictly increasing.
t0	The zero-time, at which the stochastic dynamical system is to be initialized. This must be no later than the time of the first observation, i.e., t0 <= times[1].
rprocess	optional function; a function of prototype rprocess(xstart,times,params,\dots) that simulates from the unobserved process. The form of this function is given below. pomp provides a number of plugins that construct appropriate rprocess arguments corresponding to several common stochastic simulation algorithms. See below for more details.
dprocess	optional function; a function of prototype dprocess(x,times,params,log,\dots) that evaluates the likelihood of a sequence of consecutive state transitions. The form of this function is given below. It is not typically necessary (or even feasible) to define dprocess. This functionality is provided to support future algorithm development.
rmeasure	optional; the measurement model simulator. This can be specified in one of four ways: as a function of prototype rmeasure(x,t,params,\dots) that makes a draw from the observation process given states x, time t, and parameters params. as the name of a native (compiled) routine of type pomp_measure_model_simulator as defined in the header file ‘pomp.h’. (To view the header file, execute file.show(system.file("include/pomp.h",package="pomp")) in an R session.) using the formula-based measurement.model facility (see below). as a snippet of C code (via Csnippet) that draws from the observation process as above. The last is typically the preferred option, as it results in much faster code execution.
dmeasure	optional; the measurement model probability density function. This can be specified in one of four ways: as a function of prototype dmeasure(y,x,t,params,log,\dots) that computes the p.d.f. of y given x, t, and params. as the name of a native (compiled) routine of type pomp_measure_model_density as defined in the header file ‘pomp.h’. (To view the header file, execute file.show(system.file("include/pomp.h",package="pomp")) in an R session.) using the formula-based measurement.model facility (see below). as a snippet of C code (via Csnippet) that computes the p.d.f. as above. The last is typically the preferred option, as it results in much faster code execution. As might be expected, if log=TRUE, this function should return the log likelihood.
measurement.model	optional; a formula or list of formulae, specifying the measurement model. These formulae are parsed internally to generate rmeasure and dmeasure functions. If measurement.model is given it overrides any specification of rmeasure or dmeasure. NB: This is a convenience function, primarily designed to facilitate exploration; it will typically be possible to acclerate measurement model computations by writing dmeasure and/or rmeasure using Csnippets.
skeleton, skeleton.type, skelmap.delta.t	The function skeleton specifies the deterministic skeleton of the unobserved Markov process. If we are dealing with a discrete-time Markov process, its deterministic skeleton is a map: indicate this by specifying skeleton.type="map". In this case, the default assumption is that time advances 1 unit per iteration of the map; to change this, set skelmap.delta.t to the appropriate time-step. If we are dealing with a continuous-time Markov process, its deterministic skeleton is a vectorfield: indicate this by specifying skeleton.type="vectorfield". The skeleton function can be specified in one of three ways: as an R function of prototype skeleton(x,t,params,\dots) that evaluates the deterministic skeleton at state x and time t given the parameters params, as the name of a native (compiled) routine of type pomp_skeleton as defined in the header file ‘pomp.h’. (To view the header file, execute file.show(system.file("include/pomp.h",package="pomp")) in an R session.) as a snippet of C code (via Csnippet) that performs this evaluation. The latter is typically the preferred option, for reasons of computational efficiency.
initializer	The initializer gives the parameterization of the initial state of the unobserved Markov process. Specifically, given a vector of parameters, params and an initial time, t0, the initializer determines the state vector at time t0. By default, any parameters in params, the names of which end in “.0”, are assumed to be initial values of states. To initialize the unobserved state process, these are simply copied over as initial conditions. The names of the resulting state variables are obtained by dropping the “.0” suffix. A custom initializer can be supplied here in one of two formats: as an R function of prototype initializer(params,t0,\dots) that yields initial conditions for the state process when given a vector, params, of parameters. as a snippet of C code (via Csnippet). See the Csnippet help for rules on writing a custom initializer.
rprior	optional; function drawing a sample from a prior distribution on parameters. This can be specified in one of three ways: as an R function of prototype rprior(params,\dots) that makes a draw from the prior distribution given params, as the name of a native (compiled) routine of type pomp_rprior as defined in the header file ‘pomp.h’, or (To view the header file, execute file.show(system.file("include/pomp.h",package="pomp")) in an R session.) as a snippet of C code (via Csnippet). As above, the latter is typically preferable.
dprior	optional; function evaluating the prior distribution. This can be specified in one of three ways: as an R function of prototype dprior(params,log=FALSE,\dots) that evaluates the prior probability density, as the name of a native (compiled) routine of type pomp_dprior as defined in the header file ‘pomp.h’, or (To view the header file, execute file.show(system.file("include/pomp.h",package="pomp")) in an R session.) as a snippet of C code (via Csnippet). As above, the latter is typically preferable.
params	optional named numeric vector of parameters. This will be coerced internally to storage mode double.
covar, tcovar	An optional matrix or data frame of covariates: covar is the table of covariates (one column per variable); tcovar the corresponding times (one entry per row of covar). covar can be specified as either a matrix or a data frame. In either case the columns are taken to be distinct covariates. If covar is a data frame, tcovar can be either the name or the index of the time variable. If a covariate table is supplied, then the value of each of the covariates is interpolated as needed. The resulting interpolated values are passed to the corresponding functions as a numeric vector named covars; see below.
obsnames, statenames, paramnames, covarnames	Optional character vectors specifying the names of observables, state variables, parameters, and covariates, respectively. These are only used in the event that one or more of the basic functions (rprocess, dprocess, rmeasure, dmeasure, skeleton, rprior, dprior) are defined using Csnippet or native routines.
zeronames	optional character vector specifying the names of accumulator variables (see below).
PACKAGE	An optional string giving the name of the dynamically loaded library in which any native routines are to be found.
fromEstimationScale, toEstimationScale	Optional functions specifying parameter transformations. Many algorithms for parameter estimation search an unconstrained space of parameters. When working with such an algorithm and a model for which the parameters are constrained, it can be useful to transform parameters. toEstimationScale and fromEstimationScale are transformations from the model scale to the estimation scale, and vice versa, respectively. These functions must have arguments params and .... toEstimationScale should transform parameters from the scale that rprocess, dprocess, rmeasure, dmeasure, and skeleton use internally to the scale used in estimation. fromEstimationScale should be the inverse of toEstimationScale. The parameter transformations can be defined (as above) using either R functions, native routines, or Csnippets. Note that it is the user's responsibility to make sure that these transformations are mutually inverse. If obj is the constructed pomp object, and coef(obj) is non-empty, a simple check of this property is x <- coef(obj,transform=TRUE) obj1 <- obj coef(obj1,transform=TRUE) <- x identical(coef(obj),coef(obj1)) identical(coef(obj1,transform=TRUE),x). By default, both functions are the identity transformation. See the demos, demo(package="pomp"), pompExample, and the tutorials on the package website for examples.
globals	optional character; C code that will be included in the source for (and therefore hard-coded into) the shared-object library created when the call to pomp uses Csnippets. If no Csnippets are used, globals has no effect.
...	Any additional arguments given to pomp will be stored in the pomp object and passed as arguments to each of the basic functions whenever they are evaluated.

Value

pomp returns an object of class pomp. If data is an object of class pomp, then by default the returned pomp object is identical to data. If additional arguments are given, these override the defaults.

Important note

It is not typically necessary (or even feasible) to define all of the components rprocess, dprocess, rmeasure, dmeasure, and skeleton in any given problem. Each algorithm makes use of only a subset of these components. Any algorithm requiring a component that has not been defined will return an informative error.

The state process model

Specification of process-model codes rprocess and/or dprocess in most cases is facilitated by pomp's so-called plugins, which have been developed to handle common use-cases. Currently, if one's process model evolves in discrete time or one is willing to make such an approximation (e.g., via an Euler approximation), then the euler.sim, discrete.time.sim, and onestep.sim plugins for rprocess and onestep.dens plugin for dprocess are available. In addition, for exact simulation of certain continuous-time Markov chains, an implementation of Gillespie's algorithm is available (see gillespie.sim). To learn more about the use of plugins, consult the help documentation (plugins) and the tutorials on the package website. Several of the demos and examples make use of these as well.

In specific cases, it may be possible to obtain increased computational efficiency by writing custom versions of rprocess and/or dprocess instead of using the plugins. If such custom versions are desired, the following describes how these functions should be written.

rprocess

If the plugins are not used rprocess must be an R function with at least the following arguments: xstart, times, params, and .... It can also take additional arguments. It is guaranteed that these will be filled with the corresponding elements the user has included as additional arguments in the construction of the pomp object.

In calls to rprocess, xstart can be assumed to be an nvar x nrep matrix; its rows correspond to components of the state vector and columns correspond to independent realizations of the process. params will similarly be an npar x nrep matrix with rows corresponding to parameters and columns corresponding to independent realizations. Note that the columns of params correspond to those of xstart; in particular, they will agree in number. Both xstart and params are guaranteed to have rownames.

rprocess must return a rank-3 array with rownames. Suppose x is the array returned. Then dim(x)=c(nvars,nrep,ntimes), where ntimes is the length of the vector times. x[,j,k] is the value of the state process in the j-th realization at time times[k]. In particular, x[,,1] must be identical to xstart. The rownames of x must correspond to those of xstart.

dprocess

If the plugins are not used, dprocess must have at least the following arguments: x, times, params, log, and .... It may take additional arguments: again, these will be filled with the corresponding elements the user defines when the pomp object is constructed.

In calls to dprocess, x may be assumed to be an nvars x nrep x ntimes array, where these terms have the same meanings as above. params will be a matrix with rows corresponding to individual parameters and columns corresponding to independent realizations. The columns of params correspond to those of x; in particular, they will agree in number. Both x and params are guaranteed to have rownames.

dprocess must return a nrep x ntimes-1 matrix. Suppose d is the array returned. d[j,k] is the probability density of the transition from state x[,j,k-1] at time times[k-1] to state x[,j,k] at time times[k]. If log=TRUE, then the log of the pdf must be returned.

In writing this function, you may assume that the transitions are consecutive. It should be clear that, but for this assumption, it will in general be impossible to write the transition probabilities explicitly. In such cases, algorithms that make no use of dprocess, which are said to have the “plug and play” property, are useful. Most of the algorithms in pomp have this property. In particular, at present, no methods in pomp make use of dprocess.

The observation process model

The following is a guide to writing the measurement model components as R functions. For a description on how to write these components using Csnippets, see the tutorials on the package website.

rmeasure

if provided, must take at least the arguments x, t, params, and .... It may take additional arguments, which will be filled with user-specified data as above. x will be a named numeric vector of length nvars (which has the same meaning as above). t will be a scalar quantity, the time at which the measurement is made. params will be a named numeric vector of length npars. The rmeasure function may take additional arguments which will be filled with user-specified data as above.

rmeasure must return a named numeric vector of length nobs, the number of observable variables.

dmeasure

if provided, must take at least the arguments y, x, t, params, log, and .... y will be a named numeric vector of length nobs containing (actual or simulated) values of the observed variables; x will be a named numeric vector of length nvar containing state variables; params will be a named numeric vector containing parameters; and t will be a scalar, the corresponding observation time. The dmeasure function may take additional arguments which will be filled with user-specified data as above.

dmeasure must return a single numeric value, the probability density of y given x at time t. If log=TRUE, then dmeasure should return the log of the probability density.

The deterministic skeleton

The following describes how to specify the deterministic skeleton as an R function. For a description on how to write this component using Csnippets, see the tutorials on the package website and the Csnippet help.

If skeleton if provided, must have at least the arguments x, t, params, and .... x is a numeric vector containing the coordinates of a point in state space at which evaluation of the skeleton is desired. t is a numeric value giving the time at which evaluation of the skeleton is desired. Of course, these will be irrelevant in the case of an autonomous skeleton. params is a numeric vector holding the parameters. skeleton may take additional arguments, which will be filled, as above, with user-specified data.

skeleton must return a numeric vector of the same length as x, which contains the value vectorfield (if the dynamical system is continuous) or the value of the map (if the dynamical system is discrete), at the point x at time t.

The state-process initializer

if provided, must have at least the arguments params, t0, and .... params will be a named numeric vector of parameters. t0 will be the time at which initial conditions are desired. initializer must return a named numeric vector of initial states.

Covariates

If the pomp object contains covariates (via the covar argument; see above), then whenever any of the R functions described above are called, they will each be supplied with an additional argument covars. This will be a named numeric vector containing the (interpolated) values of the covariates at the time t. In particular, covars will have one value for each column of the covariate table.

Accumulator variables

In formulating models, one often wishes to define a state variable that will accumulate some quantity over the interval between successive observations. pomp provides a facility to make such features more convenient. Specifically, variables named in the pomp's zeronames argument will be set to zero immediately following each observation. See euler.sir and the tutorials on the package website for examples.

Warning

Some error checking is done by pomp, but complete error checking is impossible. If the user-specified functions do not conform to the above specifications, then the results may be invalid. In particular, if both rmeasure and dmeasure are specified, the user should verify that these two functions correspond to the same probability distribution. If skeleton is specified, the user is responsible for verifying that it corresponds to a deterministic skeleton of the model. Each pomp-package algorithm uses some subset of the five basic functions (rprocess, dprocess, rmeasure, dmeasure, skeleton). If an algorithm requires a component that has not been specified, an informative error will be generated.

Author(s)

Aaron A. King kingaa at umich dot edu

See Also

pomp methods, pomp low-level interface, process model plugins

Examples

## Not run: 
pompExample()
pomp.home <- system.file("examples",package="pomp")
pomp.examples <- list.files(pomp.home)
file.show(
          file.path(pomp.home,pomp.examples),
          header=paste("======",pomp.examples,"=======")
         )

## End(Not run)

Example output

examples in /usr/lib/R/site-library/pomp/examples:
[1] "bbs"           "blowflies"     "dacca"         "euler.sir"    
[5] "gillespie.sir" "gompertz"      "ou2"           "ricker"       
[9] "rw2"          
====== bbs.R =======

library(pomp)

flu <- read.csv2(text="
day;reports
1;3
2;8
3;28
4;76
5;222
6;293
7;257
8;237
9;192
10;126
11;70
12;28
13;12
14;5
")

pomp(
     data=flu,
     times="day",
     t0=0,
     params=c(
       gamma=1/3,mu=0.0,iota=0.0,
       beta=1.4,
       beta.sd=0,
       pop=1400,
       rho=0.9,sigma=3.6,
       S_0=0.999,I_0=0.001,R_0=0
       ),
     rprocess=euler.sim(
       step.fun="_sir_euler_simulator",
       delta.t=1/12,
       PACKAGE="pomp"
       ),
     skeleton=vectorfield("_sir_ODE"),
     measurement.model=reports~nbinom(mu=rho*cases,size=1/sigma^2),
     PACKAGE="pomp",
     statenames=c("S","I","R","cases","W"),
     paramnames=c(
       "gamma","mu","iota",
       "beta","beta.sd","pop","rho",
       "S_0","I_0","R_0"
       ),
     zeronames=c("cases"),
     nbasis=1L,
     degree=0L,
     period=1.0,
     logvar=c(
       "beta","gamma","mu","iota","sigma","beta.sd",
       "S_0","I_0","R_0"
       ),
     logitvar="rho",
     toEstimationScale=function (params, logvar, logitvar, ...) {
       params[logvar] <- log(params[logvar])
       params[logitvar] <- qlogis(params[logitvar])
       params
     },
     fromEstimationScale=function (params, logvar, logitvar, ...) {
       params[logvar] <- exp(params[logvar])
       params[logitvar] <- plogis(params[logitvar])
       params
     },
     initializer="_sir_init"
     ) -> bbs

c("bbs")

====== blowflies.R =======

## blowfly model, with general dt
## here, set up for dt=1 and dt=2
## dt is hard-coded, and initial values are customized for each dt

library(pomp)

## following xia and tong, the delay is treated as fixed at 14 days
## xia and tong claim to be using tau=8 bidays, but on inspection 
## their Euler method is really tau=7 bidays

blowfly.data <- '"day";"y";"set"
40;3721;1
41;3373;1
42;2880;1
43;1805;1
44;1195;1
45;557;1
46;267;1
47;239;1
48;182;1
49;270;1
50;300;1
51;330;1
52;360;1
53;390;1
54;420;1
55;362;1
56;363;1
57;423;1
58;423;1
59;483;1
60;425;1
61;397;1
62;573;1
63;864;1
64;1533;1
65;2377;1
66;3366;1
67;4472;1
68;5403;1
69;7613;1
70;7498;1
71;7440;1
72;7529;1
73;7820;1
74;7996;1
75;8316;1
76;5061;1
77;5150;1
78;5296;1
79;4337;1
80;3728;1
81;3060;1
82;2654;1
83;2074;1
84;1697;1
85;1290;1
86;971;1
87;769;1
88;683;1
89;655;1
90;626;1
91;511;1
92;425;1
93;426;1
94;310;1
95;253;1
96;196;1
97;168;1
98;197;1
99;286;1
100;781;1
101;927;1
102;1073;1
103;1249;1
104;1366;1
105;1570;1
106;1949;1
107;2473;1
108;3578;1
110;4800;1
111;5440;1
112;4715;1
113;4192;1
114;3844;1
115;3177;1
116;3439;1
117;3731;1
118;4168;1
119;3617;1
120;3123;1
121;2223;1
122;1759;1
123;800;1
124;656;1
125;483;1
126;367;1
127;223;1
128;224;1
129;138;1
130;51;1
131;0;1
132;82;1
133;170;1
134;287;1
135;172;1
136;80;1
137;115;1
138;0;1
139;379;1
140;961;1
141;2095;1
142;2678;1
143;3231;1
144;3116;1
145;2914;1
146;2740;1
147;2944;1
148;3526;1
149;4225;1
150;5127;1
151;6116;1
152;7222;1
153;6554;1
154;5944;1
155;4405;1
156;3010;1
157;2255;1
158;1413;1
159;862;1
160;601;1
161;544;1
162;516;1
163;488;1
164;488;1
165;489;1
166;490;1
167;462;1
168;463;1
169;464;1
170;464;1
171;465;1
172;496;1
173;584;1
174;701;1
175;818;1
176;1400;1
177;1953;1
178;2710;1
179;3234;1
180;3729;1
181;3875;1
182;4022;1
183;3761;1
184;3559;1
185;3996;1
186;4462;1
187;4928;1
188;5248;1
189;5889;1
190;6616;1
191;7548;1
192;9642;1
193;8596;1
194;4585;1
195;2812;1
196;1883;1
197;1477;1
198;984;1
199;607;1
200;404;1
201;434;1
202;377;1
203;349;1
204;321;1
205;292;1
206;235;1
207;207;1
208;121;1
209;92;1
210;122;1
211;152;1
212;328;1
213;648;1
214;1347;1
215;1784;1
216;2337;1
217;2803;1
218;3269;1
219;3735;1
220;4376;1
221;4260;1
222;4174;1
223;4116;1
224;4321;1
225;4496;1
226;4701;1
227;4411;1
228;4034;1
229;3657;1
230;3890;1
231;4269;1
232;4939;1
233;5754;1
234;4940;1
235;4040;1
236;3227;1
237;2356;1
238;1427;1
239;1108;1
240;905;1
241;1052;1
242;1110;1
243;1169;1
244;1257;1
245;1258;1
246;1230;1
247;1173;1
248;1203;1
249;1116;1
250;1059;1
251;1002;1
252;1090;1
253;1237;1
254;1615;1
255;2168;1
256;2721;1
257;3420;1
258;4381;1
259;5254;1
260;6185;1
261;7697;1
262;7465;1
263;6478;1
264;6188;1
265;5898;1
266;5522;1
267;5057;1
268;4331;1
269;4042;1
270;3839;1
271;3607;1
272;3289;1
273;2446;1
274;1749;1
275;1023;1
276;646;1
277;444;1
278;242;1
279;214;1
280;69;1
281;70;1
282;100;1
283;101;1
284;130;1
285;248;1
286;336;1
287;482;1
288;629;1
289;862;1
290;1066;1
291;1329;1
292;1765;1
293;2202;1
294;2406;1
295;4733;1
296;4851;1
297;4502;1
298;4213;1
299;3982;1
300;3895;1
301;3663;1
302;4013;1
303;4508;1
304;5352;1
305;6079;1
306;6546;1
307;5588;1
308;4774;1
309;3932;1
310;3119;1
311;2073;1
312;1638;1
313;1348;1
314;797;1
315;536;1
240;439;2
242;244;2
244;146;2
246;78;2
248;39;2
250;10;2
252;19;2
254;156;2
256;1550;2
258;2262;2
260;2301;2
262;2067;2
264;1521;2
266;1024;2
268;1326;2
270;1238;2
272;1033;2
274;809;2
276;556;2
278;322;2
280;136;2
282;49;2
284;29;2
286;10;2
288;19;2
290;19;2
292;1521;2
294;2662;2
296;2798;2
298;2506;2
300;1774;2
302;1248;2
304;809;2
306;829;2
308;526;2
310;331;2
312;214;2
314;97;2
316;39;2
318;49;2
320;585;2
322;1521;2
324;3013;2
326;3519;2
328;3149;2
330;2554;2
332;1979;2
334;2525;2
336;2174;2
338;1638;2
340;1121;2
342;643;2
344;370;2
346;234;2
348;117;2
350;49;2
352;19;2
354;19;2
356;58;2
358;575;2
360;1716;2
362;1950;2
364;1813;2
366;1482;2
368;975;2
370;595;2
372;487;2
374;536;2
376;419;2
378;312;2
380;234;2
382;107;2
384;49;2
386;1316;2
388;1930;2
390;1813;2
392;1813;2
394;2398;2
396;3256;2
398;3334;2
400;2964;2
402;1872;2
404;1131;2
406;731;2
408;458;2
410;253;2
412;107;2
414;49;2
416;10;2
418;29;2
420;273;2
422;1199;2
424;1716;2
426;1716;2
428;1521;2
430;1199;2
432;731;2
434;419;2
436;819;2
438;643;2
440;526;2
442;234;2
444;97;2
446;721;2
448;1024;2
450;936;2
452;829;2
454;916;2
456;1833;2
458;2662;2
460;2623;2
240;33;3
242;16;3
244;6;3
246;22;3
248;378;3
250;860;3
252;1084;3
254;1249;3
256;1172;3
258;909;3
260;718;3
262;543;3
264;367;3
266;269;3
268;258;3
270;181;3
272;110;3
274;61;3
276;17;3
278;1;3
280;6;3
282;297;3
284;565;3
286;505;3
288;549;3
290;691;3
292;822;3
294;899;3
296;1195;3
298;1047;3
300;845;3
302;631;3
304;423;3
306;221;3
308;117;3
310;46;3
312;18;3
314;13;3
316;18;3
318;111;3
320;511;3
322;719;3
324;785;3
326;708;3
328;610;3
330;511;3
332;610;3
334;709;3
336;599;3
338;495;3
340;369;3
342;232;3
344;145;3
346;68;3
348;36;3
350;8;3
352;3;3
354;3;3
356;145;3
358;726;3
360;1328;3
362;1503;3
364;1459;3
366;1251;3
368;1016;3
370;748;3
372;638;3
374;491;3
376;332;3
378;190;3
380;69;3
382;25;3
384;9;3
386;4;3
388;37;3
390;80;3
392;425;3
394;847;3
396;1082;3
398;1050;3
400;803;3
402;749;3
404;935;3
406;891;3
408;716;3
410;514;3
412;344;3
414;229;3
416;141;3
418;87;3
420;65;3
422;21;3
424;26;3
426;142;3
428;727;3
430;919;3
432;979;3
434;925;3
436;782;3
438;536;3
440;350;3
442;213;3
444;115;3
446;65;3
448;22;3
450;11;3
452;0;3
454;5;3
456;334;3
458;498;3
460;761;3
0;948;4
2;942;4
4;911;4
6;858;4
8;801;4
10;676;4
12;504;4
14;397;4
16;248;4
18;146;4
20;1801;4
22;6235;4
24;5974;4
26;8921;4
28;6610;4
30;5973;4
32;5673;4
34;3875;4
36;2361;4
38;1352;4
40;1226;4
42;912;4
44;521;4
46;363;4
48;229;4
50;142;4
52;82;4
54;542;4
56;939;4
58;2431;4
60;3687;4
62;4543;4
64;4535;4
66;5441;4
68;4412;4
70;3022;4
72;2656;4
74;1967;4
76;1295;4
78;915;4
80;551;4
82;313;4
84;167;4
86;95;4
88;93;4
90;60;4
92;68;4
94;5259;4
96;6673;4
98;5441;4
100;3987;4
102;2952;4
104;3648;4
106;4222;4
108;3889;4
110;2295;4
112;1509;4
114;928;4
116;739;4
118;566;4
120;383;4
122;274;4
124;192;4
126;226;4
128;519;4
130;1224;4
132;2236;4
134;3818;4
136;6208;4
138;5996;4
140;5789;4
142;6652;4
144;7939;4
146;4868;4
148;3952;4
150;2712;4
152;1734;4
154;1224;4
156;703;4
158;508;4
160;366;4
162;279;4
164;243;4
166;343;4
168;761;4
170;1025;4
172;1221;4
174;1600;4
176;2267;4
178;3290;4
180;3471;4
182;3637;4
184;3703;4
186;4876;4
188;5364;4
190;4890;4
192;3029;4
194;1950;4
196;1225;4
198;1076;4
200;905;4
202;772;4
204;628;4
206;473;4
208;539;4
210;825;4
212;1702;4
214;2868;4
216;4473;4
218;5221;4
220;6592;4
222;6400;4
224;4752;4
226;3521;4
228;2719;4
230;1931;4
232;1500;4
234;1082;4
236;849;4
238;774;4
240;864;4
242;1308;4
244;1624;4
246;2224;4
248;2423;4
250;2959;4
252;3547;4
254;7237;4
256;5218;4
258;5311;4
260;4273;4
262;3270;4
264;2281;4
266;1549;4
268;1091;4
270;796;4
272;610;4
274;445;4
276;894;4
278;1454;4
280;2262;4
282;2363;4
284;3847;4
286;3876;4
288;3935;4
290;3479;4
292;3415;4
294;3861;4
296;3571;4
298;3113;4
300;2319;4
302;1630;4
304;1297;4
306;861;4
308;761;4
310;659;4
312;701;4
314;762;4
316;1188;4
318;1778;4
320;2428;4
322;3806;4
324;4519;4
326;5646;4
328;4851;4
330;5374;4
332;4713;4
334;7367;4
336;7236;4
338;5245;4
340;3636;4
342;2417;4
344;1258;4
346;766;4
348;479;4
350;402;4
352;248;4
354;254;4
356;604;4
358;1346;4
360;2342;4
362;3328;4
364;3599;4
366;4081;4
368;7643;4
370;7919;4
372;6098;4
374;6896;4
376;5634;4
378;5134;4
380;4188;4
382;3469;4
384;2442;4
386;1931;4
388;1790;4
390;1722;4
392;1488;4
394;1416;4
396;1369;4
398;1666;4
400;2627;4
402;2840;4
404;4044;4
406;4929;4
408;5111;4
410;3152;4
412;4462;4
414;4082;4
416;3026;4
418;1589;4
420;2075;4
422;1829;4
424;1388;4
426;1149;4
428;968;4
430;1170;4
432;1465;4
434;1676;4
436;3075;4
438;3815;4
440;4639;4
442;4424;4
444;2784;4
446;5860;4
448;5781;4
450;4897;4
452;3920;4
454;3835;4
456;3618;4
458;3050;4
460;3772;4
462;3517;4
464;3350;4
466;3018;4
468;2625;4
470;2412;4
472;2221;4
474;2619;4
476;3203;4
478;2706;4
480;2717;4
482;2175;4
484;1628;4
486;2388;4
488;3677;4
490;3156;4
492;4272;4
494;3771;4
496;4955;4
498;5584;4
500;3891;4
502;3501;4
504;4436;4
506;4369;4
508;3394;4
510;3869;4
512;2922;4
514;1843;4
516;2837;4
518;4690;4
520;5119;4
522;5839;4
524;5389;4
526;4993;4
528;4446;4
530;4851;4
532;4243;4
534;4620;4
536;4849;4
538;3664;4
540;3016;4
542;2881;4
544;3821;4
546;4300;4
548;4168;4
550;5446;4
552;5917;4
554;8579;4
556;7533;4
558;6884;4
560;4127;4
562;5546;4
564;6313;4
566;6650;4
568;6304;4
570;4842;4
572;4352;4
574;3215;4
576;2652;4
578;2330;4
580;3123;4
582;3955;4
584;4494;4
586;4780;4
588;5753;4
590;5555;4
592;5712;4
594;4786;4
596;4066;4
598;2891;4
600;3270;4
602;4404;4
604;4398;4
606;4112;4
608;4401;4
610;5779;4
612;6597;4
614;8091;4
616;11282;4
618;12446;4
620;13712;4
622;11017;4
624;14683;4
626;7258;4
628;6195;4
630;5962;4
632;4213;4
634;2775;4
636;1781;4
638;936;4
640;898;4
642;1160;4
644;3158;4
646;3386;4
648;4547;4
650;4823;4
652;4970;4
654;4940;4
656;5793;4
658;7836;4
660;4457;4
662;6901;4
664;8191;4
666;6766;4
668;5165;4
670;2919;4
672;3415;4
674;3431;4
676;3162;4
678;2525;4
680;2290;4
682;1955;4
684;1936;4
686;2384;4
688;4666;4
690;7219;4
692;8306;4
694;8027;4
696;7010;4
698;8149;4
700;8949;4
702;6105;4
704;5324;4
706;5766;4
708;6214;4
710;7007;4
712;8154;4
714;9049;4
716;6883;4
718;8103;4
720;6803;4
'

raw.data <- subset(
                   read.csv2(text=blowfly.data,comment.char="#"),
                   set==4,
                   select=-set
                   )

pomp(
     data=subset(raw.data,day>14&day<400),
     times="day",
     t0=14,
     rprocess=discrete.time.sim(
       step.fun="_blowfly_simulator_one",
       delta.t=1,
       PACKAGE="pomp"
       ),
     rmeasure="_blowfly_rmeasure",
     dmeasure="_blowfly_dmeasure",
     PACKAGE="pomp",
     paramnames=c("P","N0","delta","sigma.P","sigma.d","sigma.y"),
     statenames=c("N1","R","S","e","eps"),
     y.init=with( ## initial data
       raw.data,
       approx(x=day,y=y,xout=seq(from=0,to=14,by=1),rule=2)$y
       ),
#     y.init=c(948, 948, 942, 930, 911, 885, 858, 833.7, 801, 748.3, 676, 589.8, 504, 434.9, 397),
     toEstimationScale=function(params,...) {
       log(params)
     },
     fromEstimationScale=function(params,...) {
       exp(params)
     },
     initializer=function (params, t0, y.init, ...) {
       ntau <- length(y.init)
       n <- y.init[ntau:1]
       names(n) <- paste("N",seq_len(ntau),sep="")
       c(n,R=0,S=0,e=0,eps=0)
     }
     ) -> blowflies1

pomp(
     blowflies1,
     rprocess=discrete.time.sim(
       step.fun="_blowfly_simulator_two",
       delta.t=2,
       PACKAGE="pomp"
       ),
     y.init=with( ## initial data
       raw.data,
       approx(x=day,y=y,xout=seq(from=0,to=14,by=2),rule=2)$y
       ),
     #y.init=c(948, 942, 911, 858, 801, 676, 504, 397),
     paramnames=c("P","N0","delta","sigma.P","sigma.d","sigma.y"),
     statenames=c("N1","R","S","e","eps"),
     ) -> blowflies2

## mle from search to date
coef(blowflies1,transform=TRUE) <- c(
                  P = 1.189, 
                  delta = -1.828, 
                  N0 = 6.522, 
                  sigma.P = 0.301, 
                  sigma.d = -0.292, 
                  sigma.y = -3.625
                  )

## mle from search to date
coef(blowflies2,transform=TRUE) <- c(
                  P = 1.005, 
                  delta = -1.75, 
                  N0 = 6.685, 
                  sigma.P = 0.366, 
                  sigma.d = -0.274, 
                  sigma.y = -4.524
                  )

c("blowflies1","blowflies2")

====== dacca.R =======

## pomp object encoding the "SIRS model with seasonal reservoir" of
##   King, A. A., Ionides, E. L., Pascual, M., & Bouma, M. J.
##   Inapparent infections and cholera dynamics.
##   Nature 454:877-880 (2008)
## Data are cholera deaths and decadal census figures from the Dacca district of Bengal province, 1891-1941.
##
## Data courtesy of Menno J. Bouma, London School of Tropical Medicine & Hygiene
##
## Native codes are in the package source.

library(pomp)

nrstage <- 3L
nbasis <- 6L

mle <- c(
         gamma=20.8,eps=19.1,rho=0,
         delta=0.02, deltaI=0.06, clin=1, alpha=1,
         beta_trend=-0.00498,
         logbeta1=0.747, logbeta2=6.38, logbeta3=-3.44, logbeta4=4.23, logbeta5=3.33, logbeta6=4.55,
         logomega1=log(0.184), logomega2=log(0.0786), logomega3=log(0.0584), logomega4=log(0.00917), logomega5=log(0.000208), logomega6=log(0.0124),
         sd_beta=3.13, tau=0.23,
         S_0=0.621, I_0=0.378, Y_0=0, R1_0=0.000843, R2_0=0.000972, R3_0=1.16e-07
         )

census <- data.frame(
                     year = c(1891L, 1901L, 1911L, 1921L, 1931L, 1941L),
                     census = c(2420656L, 2649522L, 2960402L, 3125967L, 3432577L, 4222142L)
                     )

cholera <- data.frame(
                      time=seq(from=1891+1/12,to=1941,by=1/12),
                      cholera.deaths = c(
                        2641L, 939L, 905L, 1219L, 368L, 78L, 29L, 12L, 30L, 44L, 270L, 1149L, 
                        633L, 501L, 855L, 1271L, 666L, 101L, 62L, 23L, 20L, 28L, 461L, 
                        892L, 751L, 170L, 253L, 906L, 700L, 98L, 57L, 72L, 471L, 4217L, 
                        5168L, 4747L, 2380L, 852L, 1166L, 2122L, 576L, 60L, 53L, 62L, 
                        241L, 403L, 551L, 739L, 862L, 348L, 490L, 5596L, 1180L, 142L, 
                        41L, 28L, 39L, 748L, 3934L, 3562L, 587L, 311L, 1639L, 1903L, 
                        601L, 110L, 32L, 19L, 82L, 420L, 1014L, 1073L, 416L, 168L, 909L, 
                        1355L, 447L, 59L, 13L, 21L, 43L, 109L, 338L, 470L, 489L, 394L, 
                        483L, 842L, 356L, 29L, 17L, 16L, 57L, 110L, 488L, 1727L, 1253L, 
                        359L, 245L, 549L, 215L, 9L, 7L, 31L, 236L, 279L, 819L, 1728L, 
                        1942L, 1251L, 3521L, 3412L, 290L, 46L, 35L, 14L, 79L, 852L, 2951L, 
                        2656L, 607L, 172L, 325L, 2191L, 584L, 58L, 38L, 8L, 22L, 50L, 
                        380L, 2059L, 938L, 389L, 767L, 1882L, 286L, 94L, 61L, 10L, 106L, 
                        281L, 357L, 1388L, 810L, 306L, 381L, 1308L, 702L, 87L, 9L, 14L, 
                        36L, 46L, 553L, 1302L, 618L, 147L, 414L, 768L, 373L, 39L, 10L, 
                        36L, 151L, 1130L, 3437L, 4041L, 1415L, 207L, 92L, 128L, 147L, 
                        32L, 7L, 59L, 426L, 2644L, 2891L, 4249L, 2291L, 797L, 680L, 1036L, 
                        404L, 41L, 19L, 12L, 10L, 121L, 931L, 2158L, 1886L, 803L, 397L, 
                        613L, 132L, 48L, 17L, 22L, 26L, 34L, 344L, 657L, 117L, 75L, 443L, 
                        972L, 646L, 107L, 18L, 6L, 9L, 5L, 12L, 142L, 133L, 189L, 1715L, 
                        3115L, 1412L, 182L, 50L, 37L, 77L, 475L, 1730L, 1489L, 620L, 
                        190L, 571L, 1558L, 440L, 27L, 7L, 14L, 93L, 1462L, 2467L, 1703L, 
                        1262L, 458L, 453L, 717L, 232L, 26L, 16L, 18L, 9L, 78L, 353L, 
                        897L, 777L, 404L, 799L, 2067L, 613L, 98L, 19L, 26L, 47L, 171L, 
                        767L, 1896L, 887L, 325L, 816L, 1653L, 355L, 85L, 54L, 88L, 609L, 
                        882L, 1363L, 2178L, 580L, 396L, 1493L, 2154L, 683L, 78L, 19L, 
                        10L, 27L, 88L, 1178L, 1862L, 611L, 478L, 2697L, 3395L, 520L, 
                        67L, 41L, 36L, 209L, 559L, 971L, 2144L, 1099L, 494L, 586L, 508L, 
                        269L, 27L, 19L, 21L, 12L, 22L, 333L, 676L, 487L, 262L, 535L, 
                        979L, 170L, 25L, 9L, 19L, 13L, 45L, 229L, 673L, 432L, 107L, 373L, 
                        1126L, 339L, 19L, 11L, 3L, 15L, 101L, 539L, 709L, 200L, 208L, 
                        926L, 1783L, 831L, 103L, 37L, 17L, 33L, 179L, 426L, 795L, 481L, 
                        491L, 773L, 936L, 325L, 101L, 22L, 25L, 24L, 88L, 633L, 513L, 
                        298L, 93L, 687L, 1750L, 356L, 33L, 2L, 18L, 70L, 648L, 2471L, 
                        1270L, 616L, 193L, 706L, 1372L, 668L, 107L, 58L, 21L, 23L, 93L, 
                        318L, 867L, 332L, 118L, 437L, 2233L, 491L, 27L, 7L, 21L, 96L, 
                        360L, 783L, 1492L, 550L, 176L, 633L, 922L, 267L, 91L, 42L, 4L, 
                        10L, 7L, 43L, 377L, 563L, 284L, 298L, 625L, 131L, 35L, 12L, 8L, 
                        9L, 83L, 502L, 551L, 256L, 198L, 664L, 1701L, 425L, 76L, 17L, 
                        9L, 16L, 5L, 141L, 806L, 1603L, 587L, 530L, 771L, 511L, 97L, 
                        35L, 39L, 156L, 1097L, 1233L, 1418L, 1125L, 420L, 1592L, 4169L, 
                        1535L, 371L, 139L, 55L, 85L, 538L, 1676L, 1435L, 804L, 370L, 
                        477L, 394L, 306L, 132L, 84L, 87L, 53L, 391L, 1541L, 1859L, 894L, 
                        326L, 853L, 1891L, 1009L, 131L, 77L, 63L, 66L, 33L, 178L, 1003L, 
                        1051L, 488L, 911L, 1806L, 837L, 280L, 132L, 76L, 381L, 1328L, 
                        2639L, 2164L, 1082L, 326L, 254L, 258L, 119L, 106L, 93L, 29L, 
                        17L, 17L, 17L, 46L, 79L, 135L, 1290L, 2240L, 561L, 116L, 24L, 
                        15L, 33L, 18L, 16L, 38L, 26L, 45L, 151L, 168L, 57L, 32L, 29L, 
                        27L, 20L, 106L, 1522L, 2013L, 434L, 205L, 528L, 634L, 195L, 45L, 
                        33L, 19L, 20L, 46L, 107L, 725L, 572L, 183L, 2199L, 4018L, 428L, 
                        67L, 31L, 8L, 44L, 484L, 1324L, 2054L, 467L, 216L, 673L, 887L, 
                        353L, 73L, 46L, 15L, 20L, 27L, 25L, 38L, 158L, 312L, 1226L, 1021L, 
                        222L, 90L, 31L, 93L, 368L, 657L, 2208L, 2178L, 702L, 157L, 317L, 
                        146L, 63L, 27L, 22L, 23L, 28L, 225L, 483L, 319L, 120L, 59L, 274L, 
                        282L, 155L, 31L, 16L, 15L, 12L, 14L, 14L, 42L
                        )
                      )

covar.dt <- 0.01

t0 <- with(cholera,2*time[1]-time[2])
tcovar <- seq(from=t0,to=max(cholera$time)+2/12,by=covar.dt)
covartable <- data.frame(
                         time=tcovar,
                         periodic.bspline.basis(tcovar-1/12,nbasis=nbasis,degree=3,period=1,names="seas%d"),
                         pop=predict(smooth.spline(x=census$year,y=census$census),x=tcovar)$y,
                         dpopdt=predict(smooth.spline(x=census$year,y=census$census),x=tcovar,deriv=1)$y,
                         trend=tcovar-mean(tcovar)
                         )

to_est <- Csnippet("
  Ttau = log(tau);
  Tgamma = log(gamma);
  Teps = log(eps);
  Tdelta = log(delta);
  TdeltaI = log(deltaI);
  Tsd_beta = log(sd_beta);
  Talpha = log(alpha);
  Trho = log(rho);
  Tclin = logit(clin);
  to_log_barycentric(&TS_0,&S_0,nrstage+3);
")
 
from_est <- Csnippet("
  Ttau = exp(tau);
  Tgamma = exp(gamma);
  Teps = exp(eps);
  Tdelta = exp(delta);
  TdeltaI = exp(deltaI);
  Tsd_beta = exp(sd_beta);
  Talpha = exp(alpha);
  Trho = exp(rho);
  Tclin = expit(clin);
  from_log_barycentric(&TS_0,&S_0,nrstage+3);
")

rinit <- Csnippet("
  int k;
  double sum = S_0+I_0+Y_0;
  double *R = &R1;
  const double *R0 = &R1_0;
  for (k = 0; k < nrstage; k++) sum += R0[k];
  S = nearbyint(pop*S_0/sum);
  I = nearbyint(pop*I_0/sum);
  Y = nearbyint(pop*Y_0/sum);
  for (k = 0; k < nrstage; k++) R[k] = nearbyint(pop*R0[k]/sum);
  W = 0;
  deaths = 0;
  count = 0;
")

norm_rmeas <- Csnippet("
  double v, tol = 1.0e-18;
  v = deaths*tau;
  if ((count > 0) || (!(R_FINITE(v)))) {
    cholera_deaths = R_NaReal;
  } else {
    cholera_deaths = rnorm(deaths,v+tol);
  }
")

norm_dmeas <- Csnippet("
  double v, tol = 1.0e-18;
  v = deaths*tau;
  if ((count>0.0) || (!(R_FINITE(v)))) {
    lik = tol;
  } else {
    lik = dnorm(cholera_deaths,deaths,v+tol,0)+tol;
  }
  if (give_log) lik = log(lik);
")

## two-path SIRS cholera model using SDEs
## exponent (alpha) on I/n
## only "severe" infections are infectious
## truncation is not used
## instead, particles with negative states are killed

cholmodel_one <- Csnippet("
  double births;
  double infections;
  double sdeaths;
  double ideaths;
  double ydeaths;
  double rdeaths[nrstage];
  double disease;
  double wanings;
  double passages[nrstage+1];
  double effI;
  double neps;
  double beta;
  double omega;
  double dw;
  double *pt;
  int j;

  if (count != 0.0) return;

  neps = eps*nrstage;

  beta = exp(dot_product(nbasis,&seas1,&logbeta1)+beta_trend*trend);
  omega = exp(dot_product(nbasis,&seas1,&logomega1));

  dw = rnorm(0,sqrt(dt));	// white noise

  effI = pow(I/pop,alpha);
  births = dpopdt + delta*pop;	// births

  passages[0] = gamma*I;	// recovery
  ideaths = delta*I;	        // natural i deaths
  disease = deltaI*I;	        // disease death
  ydeaths = delta*Y;     	// natural rs deaths
  wanings = rho*Y;		// loss of immunity

  for (pt = &R1, j = 0; j < nrstage; j++, pt++) {
    rdeaths[j] = *pt*delta;	// natural R deaths
    passages[j+1] = *pt*neps;	// passage to the next immunity class
  }

  infections = (omega+(beta+sd_beta*dw/dt)*effI)*S; // infection
  sdeaths = delta*S;	        // natural S deaths

  S += (births - infections - sdeaths + passages[nrstage] + wanings)*dt;
  I += (clin*infections - disease - ideaths - passages[0])*dt;
  Y += ((1-clin)*infections - ydeaths - wanings)*dt;
  for (pt = &R1, j = 0; j < nrstage; j++, pt++) 
    *pt += (passages[j] - passages[j+1] - rdeaths[j])*dt;
  deaths += disease*dt;		// cumulative deaths due to disease
  W += dw;

  // check for violations of positivity constraints
  // nonzero 'count' variable signals violation
  if (S < 0.0) {
    S = 0.0; I = 0.0; Y = 0.0; 
    count += 1; 
  }
  if (I < 0.0) {
    I = 0.0; S = 0.0; 
    count += 1e3; 
  }
  if (Y < 0.0) { 
    Y = 0.0; S = 0.0; 
    count += 1e6; 
  }
  if (deaths < 0.0) { 
    deaths = 0.0; 
    count += 1e9; 
  }
  for (pt = &R1, j = 0; j < nrstage-1; j++, pt++) {
    if (*pt < 0.0) {
      *pt = 0.0; *(pt+1) = 0.0;
      count += 1e12; 
    }
  }
  if (*pt < 0.0) {
    *pt = 0.0; S = 0.0;
    count += 1e12;
  }
")

pomp(
     data=cholera,
     times='time',
     t0=t0,
     params=mle,
     globals = sprintf("int nrstage = %d, nbasis = %d;",nrstage,nbasis),
     rprocess = euler.sim(
       step.fun = cholmodel_one,
       delta.t=1/240
       ),
     dmeasure = norm_dmeas,
     rmeasure=norm_rmeas,
     covar=covartable,
     tcovar='time',
     zeronames = c("deaths","count"),
     statenames = c("S","I","Y",sprintf("R%d",seq_len(nrstage)),"deaths","W","count"),
     paramnames = c("tau","gamma","eps","delta","deltaI",
       "logomega1","sd_beta","beta_trend","logbeta1",
       "alpha","rho","clin","S_0","I_0","Y_0","R1_0"),
     fromEstimationScale=from_est,
     toEstimationScale=to_est,
     initializer=rinit
     ) -> dacca

c("dacca")

====== euler.sir.R =======

library(pomp)

dat <- '"time";"reports"
0,0192307692307692;617
0,0384615384615385;638
0,0576923076923077;603
0,0769230769230769;655
0,0961538461538462;585
0,115384615384615;677
0,134615384615385;686
0,153846153846154;674
0,173076923076923;710
0,192307692307692;798
0,211538461538462;803
0,230769230769231;795
0,25;892
0,269230769230769;925
0,288461538461538;1000
0,307692307692308;1002
0,326923076923077;1088
0,346153846153846;1012
0,365384615384615;1043
0,384615384615385;1095
0,403846153846154;998
0,423076923076923;1045
0,442307692307692;1034
0,461538461538462;904
0,480769230769231;938
0,5;854
0,519230769230769;831
0,538461538461539;734
0,557692307692308;678
0,576923076923077;715
0,596153846153846;593
0,615384615384615;581
0,634615384615385;549
0,653846153846154;509
0,673076923076923;419
0,692307692307692;389
0,711538461538462;323
0,730769230769231;305
0,75;258
0,769230769230769;255
0,788461538461539;233
0,807692307692308;208
0,826923076923077;197
0,846153846153846;168
0,865384615384616;171
0,884615384615385;160
0,903846153846154;138
0,923076923076923;138
0,942307692307692;134
0,961538461538462;100
0,980769230769231;110
1;126
1,01923076923077;120
1,03846153846154;115
1,05769230769231;125
1,07692307692308;109
1,09615384615385;98
1,11538461538462;112
1,13461538461538;136
1,15384615384615;123
1,17307692307692;123
1,19230769230769;113
1,21153846153846;155
1,23076923076923;132
1,25;158
1,26923076923077;139
1,28846153846154;139
1,30769230769231;161
1,32692307692308;171
1,34615384615385;156
1,36538461538462;199
1,38461538461538;183
1,40384615384615;221
1,42307692307692;200
1,44230769230769;221
1,46153846153846;227
1,48076923076923;200
1,5;233
1,51923076923077;221
1,53846153846154;225
1,55769230769231;228
1,57692307692308;223
1,59615384615385;215
1,61538461538462;233
1,63461538461538;227
1,65384615384615;220
1,67307692307692;208
1,69230769230769;199
1,71153846153846;243
1,73076923076923;182
1,75;197
1,76923076923077;183
1,78846153846154;176
1,80769230769231;191
1,82692307692308;173
1,84615384615385;172
1,86538461538462;171
1,88461538461538;224
1,90384615384615;196
1,92307692307692;212
1,94230769230769;210
1,96153846153846;204
1,98076923076923;189
2;224
2,01923076923077;216
2,03846153846154;198
2,05769230769231;261
2,07692307692308;267
2,09615384615385;292
2,11538461538462;284
2,13461538461538;311
2,15384615384615;368
2,17307692307692;374
2,19230769230769;473
2,21153846153846;482
2,23076923076923;506
2,25;626
2,26923076923077;644
2,28846153846154;675
2,30769230769231;826
2,32692307692308;917
2,34615384615385;1014
2,36538461538462;1083
2,38461538461538;1174
2,40384615384615;1265
2,42307692307692;1244
2,44230769230769;1390
2,46153846153846;1468
2,48076923076923;1507
2,5;1468
2,51923076923077;1455
2,53846153846154;1470
2,55769230769231;1333
2,57692307692308;1294
2,59615384615385;1207
2,61538461538462;1126
2,63461538461538;1025
2,65384615384615;979
2,67307692307692;942
2,69230769230769;828
2,71153846153846;738
2,73076923076923;707
2,75;594
2,76923076923077;532
2,78846153846154;525
2,80769230769231;493
2,82692307692308;452
2,84615384615385;424
2,86538461538462;369
2,88461538461538;355
2,90384615384615;330
2,92307692307692;336
2,94230769230769;270
2,96153846153846;254
2,98076923076923;250
3;238
3,01923076923077;230
3,03846153846154;238
3,05769230769231;241
3,07692307692308;242
3,09615384615385;249
3,11538461538462;218
3,13461538461538;240
3,15384615384615;266
3,17307692307692;237
3,19230769230769;233
3,21153846153846;213
3,23076923076923;224
3,25;244
3,26923076923077;249
3,28846153846154;280
3,30769230769231;262
3,32692307692308;292
3,34615384615385;282
3,36538461538462;281
3,38461538461538;298
3,40384615384615;280
3,42307692307692;343
3,44230769230769;296
3,46153846153846;281
3,48076923076923;297
3,5;323
3,51923076923077;287
3,53846153846154;276
3,55769230769231;259
3,57692307692308;238
3,59615384615385;259
3,61538461538462;273
3,63461538461538;206
3,65384615384615;240
3,67307692307692;219
3,69230769230769;227
3,71153846153846;206
3,73076923076923;196
3,75;191
3,76923076923077;180
3,78846153846154;166
3,80769230769231;165
3,82692307692308;180
3,84615384615385;157
3,86538461538462;152
3,88461538461538;154
3,90384615384615;150
3,92307692307692;157
3,94230769230769;147
3,96153846153846;152
3,98076923076923;149
4;168
'

pomp(
     data=read.csv2(text=dat),
     times="time",
     t0=0,
     params=c(
       gamma=26,mu=0.02,iota=0.01,
       beta1=400,beta2=480,beta3=320,
       beta.sd=1e-3,
       pop=2.1e6,
       rho=0.6,
       S_0=26/400,I_0=0.001,R_0=1-26/400
       ),
     rprocess=euler.sim(
       step.fun="_sir_euler_simulator",
       delta.t=1/52/20,
       PACKAGE="pomp"
       ),
     skeleton=vectorfield("_sir_ODE"),
     rmeasure="_sir_binom_rmeasure",
     dmeasure="_sir_binom_dmeasure",
     PACKAGE="pomp",
     statenames=c("S","I","R","cases","W"),
     paramnames=c(
       "gamma","mu","iota",
       "beta1","beta.sd","pop","rho",
       "S_0","I_0","R_0"
       ),
     zeronames=c("cases"),
     fromEstimationScale="_sir_par_trans",
     toEstimationScale="_sir_par_untrans",
     nbasis=3L,
     degree=3L,
     period=1.0,
     initializer="_sir_init"
     ) -> euler.sir

## the following was originally used to generate the data
## simulate(po,nsim=1,seed=329343545L) -> euler.sir

c("euler.sir")

====== gillespie.sir.R =======

library(pomp)

dat <- '"time";"reports"
0;0
0,0192307692307692;5
0,0384615384615385;5
0,0576923076923077;2
0,0769230769230769;2
0,0961538461538462;3
0,115384615384615;2
0,134615384615385;2
0,153846153846154;1
0,173076923076923;3
0,192307692307692;2
0,211538461538462;3
0,230769230769231;1
0,25;1
0,269230769230769;2
0,288461538461538;1
0,307692307692308;2
0,326923076923077;3
0,346153846153846;1
0,365384615384615;0
0,384615384615385;0
0,403846153846154;2
0,423076923076923;1
0,442307692307692;0
0,461538461538462;1
0,480769230769231;2
0,5;0
0,519230769230769;0
0,538461538461539;0
0,557692307692308;1
0,576923076923077;2
0,596153846153846;0
0,615384615384615;1
0,634615384615385;2
0,653846153846154;0
0,673076923076923;2
0,692307692307692;2
0,711538461538462;3
0,730769230769231;1
0,75;1
0,769230769230769;0
0,788461538461539;1
0,807692307692308;1
0,826923076923077;0
0,846153846153846;3
0,865384615384615;2
0,884615384615385;10
0,903846153846154;3
0,923076923076923;4
0,942307692307692;2
0,961538461538462;2
0,980769230769231;11
1;4
1,01923076923077;8
1,03846153846154;4
1,05769230769231;1
1,07692307692308;7
1,09615384615385;7
1,11538461538462;5
1,13461538461538;3
1,15384615384615;2
1,17307692307692;7
1,19230769230769;9
1,21153846153846;9
1,23076923076923;5
1,25;3
1,26923076923077;4
1,28846153846154;12
1,30769230769231;3
1,32692307692308;11
1,34615384615385;2
1,36538461538462;10
1,38461538461538;9
1,40384615384615;5
1,42307692307692;10
1,44230769230769;10
1,46153846153846;8
1,48076923076923;7
1,5;7
1,51923076923077;5
1,53846153846154;10
1,55769230769231;10
1,57692307692308;13
1,59615384615385;14
1,61538461538462;22
1,63461538461538;18
1,65384615384615;22
1,67307692307692;16
1,69230769230769;27
1,71153846153846;24
1,73076923076923;24
1,75;36
1,76923076923077;42
1,78846153846154;35
1,80769230769231;43
1,82692307692308;50
1,84615384615385;62
1,86538461538462;47
1,88461538461538;58
1,90384615384615;50
1,92307692307692;57
1,94230769230769;51
1,96153846153846;50
1,98076923076923;54
2;61
2,01923076923077;64
2,03846153846154;58
2,05769230769231;61
2,07692307692308;80
2,09615384615385;48
2,11538461538462;61
2,13461538461538;59
2,15384615384615;59
2,17307692307692;57
2,19230769230769;39
2,21153846153846;42
2,23076923076923;45
2,25;55
2,26923076923077;26
2,28846153846154;40
2,30769230769231;50
2,32692307692308;40
2,34615384615385;49
2,36538461538462;44
2,38461538461538;57
2,40384615384615;44
2,42307692307692;56
2,44230769230769;40
2,46153846153846;54
2,48076923076923;59
2,5;67
2,51923076923077;52
2,53846153846154;51
2,55769230769231;57
2,57692307692308;66
2,59615384615385;69
2,61538461538462;61
2,63461538461538;63
2,65384615384615;55
2,67307692307692;59
2,69230769230769;59
2,71153846153846;75
2,73076923076923;62
2,75;66
2,76923076923077;80
2,78846153846154;53
2,80769230769231;69
2,82692307692308;53
2,84615384615385;56
2,86538461538462;50
2,88461538461538;45
2,90384615384615;41
2,92307692307692;40
2,94230769230769;32
2,96153846153846;39
2,98076923076923;35
3;25
3,01923076923077;32
3,03846153846154;15
3,05769230769231;20
3,07692307692308;19
3,09615384615385;16
3,11538461538462;11
3,13461538461538;14
3,15384615384615;16
3,17307692307692;11
3,19230769230769;8
3,21153846153846;6
3,23076923076923;3
3,25;8
3,26923076923077;8
3,28846153846154;2
3,30769230769231;3
3,32692307692308;4
3,34615384615385;8
3,36538461538462;7
3,38461538461539;4
3,40384615384615;9
3,42307692307692;7
3,44230769230769;8
3,46153846153846;5
3,48076923076923;5
3,5;7
3,51923076923077;1
3,53846153846154;6
3,55769230769231;7
3,57692307692308;6
3,59615384615385;5
3,61538461538462;8
3,63461538461539;10
3,65384615384615;1
3,67307692307692;0
3,69230769230769;5
3,71153846153846;8
3,73076923076923;1
3,75;8
3,76923076923077;8
3,78846153846154;5
3,80769230769231;7
3,82692307692308;3
3,84615384615385;7
3,86538461538462;4
3,88461538461539;6
3,90384615384615;7
3,92307692307692;6
3,94230769230769;7
3,96153846153846;7
3,98076923076923;5
4;5
4,01923076923077;7
4,03846153846154;2
4,05769230769231;8
4,07692307692308;3
4,09615384615385;3
4,11538461538462;5
4,13461538461539;1
4,15384615384615;3
4,17307692307692;5
4,19230769230769;5
4,21153846153846;2
4,23076923076923;3
4,25;6
4,26923076923077;1
4,28846153846154;2
4,30769230769231;2
4,32692307692308;5
4,34615384615385;3
4,36538461538462;3
4,38461538461539;7
4,40384615384615;2
4,42307692307692;2
4,44230769230769;2
4,46153846153846;3
4,48076923076923;4
4,5;4
4,51923076923077;7
4,53846153846154;1
4,55769230769231;5
4,57692307692308;12
4,59615384615385;5
4,61538461538462;8
4,63461538461539;3
4,65384615384615;17
4,67307692307692;14
4,69230769230769;15
4,71153846153846;12
4,73076923076923;19
4,75;19
4,76923076923077;17
4,78846153846154;23
4,80769230769231;28
4,82692307692308;23
4,84615384615385;32
4,86538461538462;24
4,88461538461539;35
4,90384615384615;31
4,92307692307692;23
4,94230769230769;41
4,96153846153846;27
4,98076923076923;28
5;26
5,01923076923077;29
5,03846153846154;38
5,05769230769231;20
5,07692307692308;17
5,09615384615385;29
5,11538461538462;31
5,13461538461539;25
5,15384615384615;31
5,17307692307692;26
5,19230769230769;23
5,21153846153846;36
5,23076923076923;25
5,25;25
5,26923076923077;35
5,28846153846154;37
5,30769230769231;32
5,32692307692308;30
5,34615384615385;41
5,36538461538462;45
5,38461538461539;61
5,40384615384615;40
5,42307692307692;43
5,44230769230769;46
5,46153846153846;49
5,48076923076923;59
5,5;48
5,51923076923077;69
5,53846153846154;61
5,55769230769231;68
5,57692307692308;79
5,59615384615385;91
5,61538461538462;93
5,63461538461539;77
5,65384615384615;87
5,67307692307692;107
5,69230769230769;103
5,71153846153846;119
5,73076923076923;111
5,75;110
5,76923076923077;91
5,78846153846154;120
5,80769230769231;106
5,82692307692308;93
5,84615384615385;69
5,86538461538462;90
5,88461538461539;100
5,90384615384615;73
5,92307692307692;72
5,94230769230769;67
5,96153846153846;63
5,98076923076923;41
6;43
6,01923076923077;33
6,03846153846154;40
6,05769230769231;29
6,07692307692308;25
6,09615384615385;22
6,11538461538462;19
6,13461538461539;15
6,15384615384615;18
6,17307692307692;13
6,19230769230769;15
6,21153846153846;14
6,23076923076923;10
6,25;8
6,26923076923077;13
6,28846153846154;4
6,30769230769231;10
6,32692307692308;3
6,34615384615385;9
6,36538461538462;10
6,38461538461539;9
6,40384615384615;4
6,42307692307692;4
6,44230769230769;3
6,46153846153846;2
6,48076923076923;6
6,5;4
6,51923076923077;5
6,53846153846154;3
6,55769230769231;8
6,57692307692308;3
6,59615384615385;6
6,61538461538462;5
6,63461538461539;6
6,65384615384615;1
6,67307692307692;8
6,69230769230769;12
6,71153846153846;10
6,73076923076923;3
6,75;3
6,76923076923077;6
6,78846153846154;5
6,80769230769231;10
6,82692307692308;6
6,84615384615385;6
6,86538461538462;5
6,88461538461539;11
6,90384615384615;6
6,92307692307692;14
6,94230769230769;7
6,96153846153846;6
6,98076923076923;5
7;4
7,01923076923077;7
7,03846153846154;3
7,05769230769231;4
7,07692307692308;3
7,09615384615385;3
7,11538461538462;2
7,13461538461539;3
7,15384615384615;7
7,17307692307692;2
7,19230769230769;2
7,21153846153846;3
7,23076923076923;3
7,25;5
7,26923076923077;2
7,28846153846154;3
7,30769230769231;7
7,32692307692308;2
7,34615384615385;6
7,36538461538462;2
7,38461538461539;5
7,40384615384615;5
7,42307692307692;6
7,44230769230769;5
7,46153846153846;8
7,48076923076923;9
7,5;5
7,51923076923077;5
7,53846153846154;15
7,55769230769231;9
7,57692307692308;18
7,59615384615385;17
7,61538461538462;10
7,63461538461539;16
7,65384615384615;15
7,67307692307692;19
7,69230769230769;14
7,71153846153846;23
7,73076923076923;28
7,75;32
7,76923076923077;32
7,78846153846154;31
7,80769230769231;40
7,82692307692308;37
7,84615384615385;41
7,86538461538462;37
7,88461538461539;40
7,90384615384615;38
7,92307692307692;44
7,94230769230769;35
7,96153846153846;34
7,98076923076923;26
8;49
8,01923076923077;32
8,03846153846154;45
8,05769230769231;38
8,07692307692308;31
8,09615384615385;26
8,11538461538461;27
8,13461538461539;26
8,15384615384616;26
8,17307692307692;22
8,19230769230769;33
8,21153846153846;36
8,23076923076923;27
8,25;36
8,26923076923077;39
8,28846153846154;36
8,30769230769231;30
8,32692307692308;39
8,34615384615385;28
8,36538461538461;31
8,38461538461539;40
8,40384615384616;47
8,42307692307692;41
8,44230769230769;42
8,46153846153846;42
8,48076923076923;50
8,5;60
8,51923076923077;53
8,53846153846154;64
8,55769230769231;59
8,57692307692308;65
8,59615384615385;61
8,61538461538462;69
8,63461538461539;66
8,65384615384616;71
8,67307692307692;71
8,69230769230769;72
8,71153846153846;94
8,73076923076923;79
8,75;79
8,76923076923077;78
8,78846153846154;78
8,80769230769231;92
8,82692307692308;64
8,84615384615385;70
8,86538461538462;50
8,88461538461539;63
8,90384615384616;58
8,92307692307692;53
8,94230769230769;48
8,96153846153846;51
8,98076923076923;28
9;36
9,01923076923077;40
9,03846153846154;20
9,05769230769231;16
9,07692307692308;19
9,09615384615385;16
9,11538461538462;20
9,13461538461539;16
9,15384615384616;11
9,17307692307692;9
9,19230769230769;15
9,21153846153846;15
9,23076923076923;7
9,25;12
9,26923076923077;9
9,28846153846154;5
9,30769230769231;4
9,32692307692308;3
9,34615384615385;4
9,36538461538462;4
9,38461538461539;5
9,40384615384616;3
9,42307692307692;2
9,44230769230769;2
9,46153846153846;3
9,48076923076923;3
9,5;5
9,51923076923077;3
9,53846153846154;3
9,55769230769231;6
9,57692307692308;6
9,59615384615385;5
9,61538461538462;2
9,63461538461539;5
9,65384615384616;5
9,67307692307692;3
9,69230769230769;3
9,71153846153846;3
9,73076923076923;4
9,75;3
9,76923076923077;9
9,78846153846154;4
9,80769230769231;2
9,82692307692308;5
9,84615384615385;4
9,86538461538462;2
9,88461538461539;4
9,90384615384616;5
9,92307692307692;3
9,94230769230769;2
9,96153846153846;3
9,98076923076923;2
10;6
'

pomp(
     data=read.csv2(text=dat),
     times="time",
     t0=0,
     params=c(
       gamma=24,mu=1/70,iota=0.1,
       beta1=330,beta2=410,beta3=490,
       rho=0.1,
       S_0=0.05,I_0=1e-4,R_0=0.95,
       pop=1000000,
       beta.sd=0
       ),
     rprocess=gillespie.sim(
       rate.fun="_sir_rates",
       PACKAGE="pomp",
       v=cbind(
         birth=c(1,0,0,1,0),
         sdeath=c(-1,0,0,-1,0),
         infection=c(-1,1,0,0,0),
         ideath=c(0,-1,0,-1,0),
         recovery=c(0,-1,1,0,1),
         rdeath=c(0,0,-1,-1,0)
         )
       ),
     skeleton=vectorfield("_sir_ODE"),
     measurement.model=reports~binom(size=cases,prob=rho),
     PACKAGE="pomp",
     statenames=c("S","I","R","N","cases"),
     paramnames=c(
       "gamma","mu","iota",
       "beta1","beta.sd","pop","rho",
       "S_0","I_0","R_0"
       ),
     zeronames=c("cases"),
     comp.names=c("S","I","R"),
     ic.names=c("S_0","I_0","R_0"),
     fromEstimationScale="_sir_par_trans",
     toEstimationScale="_sir_par_untrans",
     nbasis=3L,
     degree=3L,
     period=1.0,
     initializer=function(params, t0, comp.names, ic.names, ...) {
       x0 <- setNames(numeric(5),c("S","I","R","N","cases"))
       fracs <- params[ic.names]
       x0["N"] <- params["pop"]
       x0[comp.names] <- round(params["pop"]*fracs/sum(fracs))
       x0
     }
     ) -> gillespie.sir

## originally, the data were created via:
## simulate(po,nsim=1,seed=1165270654L) -> gillespie.sir

c("gillespie.sir")

====== gompertz.R =======

library(pomp)

dat <- '"time";"X";"Y"
0;1;1,07952841598157
1;0,950636618475229;1,11861594539405
2;0,83244916029991;0,828286900987495
3;0,93103099656842;1,13600466557311
4;0,911439385512777;0,806676551188204
5;0,837105441016674;0,873084434929483
6;1,00268501398128;1,04124835021771
7;1,06579807727341;1,05979968381669
8;1,04671681272393;1,14368595456944
9;1,01143788899437;1,17581280964431
10;1,12153364435953;1,16125775736251
11;1,23274218227861;1,22891596556583
12;1,32655395564463;1,38234823254117
13;1,35069512981151;1,13856582462482
14;1,13875163542966;1,10278598215362
15;1,26066897039704;1,29730413792106
16;1,24543791707035;1,33112792763266
17;1,23134905281316;1,20013316479374
18;1,15533668996193;1,1748481850297
19;1,19423810951996;1,37677723538777
20;1,20635310746995;1,30467220974374
21;1,04307227953634;0,933048055379677
22;1,04107912990884;1,08563062740614
23;1,00318246871677;1,02180297594675
24;0,970849938478495;1,25099936514459
25;0,79243164319502;0,79761814762639
26;0,771762127953167;0,68253288426557
27;0,760988977405389;0,785573183621746
28;0,836517935205031;0,880032437134315
29;0,7939171601806;0,817662333089958
30;0,723293224297652;0,682204778379272
31;0,775999242794153;0,609158268599648
32;0,801813842772822;0,953075868457542
33;0,76994574267228;0,64000611215012
34;0,746021978630474;0,824778708853423
35;0,81336335287903;0,75477831785326
36;0,923882876057083;1,05020086682025
37;1,02216853265321;1,07915643829565
38;0,953397077328701;0,910377207219948
39;1,01222135734455;1,07196221241156
40;0,978288274555936;0,927402543398725
41;1,07783628249536;1,13562270683309
42;1,18185012207163;1,15138433126514
43;0,952520534180798;0,979766621130238
44;0,84439364958955;0,932583376950161
45;0,85025855989352;0,874190322172554
46;0,868591242831582;0,917317195376153
47;0,96734181848679;0,810217908784557
48;0,90424208054782;0,91751282983596
49;0,77998160754507;0,777229472252738
50;0,807245881695201;0,984536558782125
51;0,796994124555762;0,829289645018853
52;0,811536735801658;0,738549113698673
53;0,817376297833073;0,751664773418347
54;0,911274275731693;0,901320737117393
55;1,09671858179842;1,02114400728858
56;1,17782304259497;1,15480952716325
57;1,48974254173377;1,52171146422867
58;1,59363218853134;1,60303998339523
59;1,33894250041079;1,3673888909755
60;1,24282152770653;1,07557075960246
61;1,21473587944885;1,26679772028961
62;1,31898738946143;1,22392737161959
63;1,45195233453426;1,40639369164472
64;1,58364437767931;1,61885414154031
65;1,60727226980948;1,70623639556045
66;1,55648312498095;1,87962959035064
67;1,66568956879442;1,73633383185204
68;1,71151559987914;1,59157235244771
69;1,53634236355877;1,58129601393277
70;1,62285253887348;1,66169160756609
71;1,27898023105463;1,35952077984963
72;1,1818761514106;1,22549784715231
73;1,09002630947639;1,2025122303878
74;1,2232879353834;1,25096891846816
75;1,43006968443063;1,41018069268411
76;1,41189860111872;1,43675715990225
77;1,55201995283258;1,44229999201853
78;1,45818570565819;1,38883706742128
79;1,22002519977696;1,25058856058262
80;1,23865238275549;1,23288019435894
81;1,48621496132753;1,58848483557968
82;1,36483886467245;1,23628960315986
83;1,51299070383272;1,60393599831722
84;1,25759262244331;1,1541357227336
85;1,19746015105645;1,35609775043505
86;1,05742547057021;0,855756862528805
87;1,22304508672589;1,38718382523001
88;1,25379301734834;1,24797305664899
89;1,13949050640685;1,35307215916429
90;0,817425730431819;0,970652163066452
91;0,802233823759296;0,847914335870894
92;0,857912063302802;0,897827069767532
93;1,05787680391846;1,16338321886402
94;1,05448994160666;1,04849780806396
95;1,11773425002073;1,02168552045925
96;1,03321705886828;1,23012631227719
97;0,985172093658399;1,06931864474828
98;0,942026595214108;1,14999282942163
99;1,16257505590149;0,972309774206453
100;1,28071364632546;1,23495259845453
'

pomp(
     data=subset(read.csv2(text=dat),select=-X),
     times="time",
     t0=0,
     params=c(K=1,r=0.1,sigma=0.1,tau=0.1,X.0=1),
     rprocess=discrete.time.sim(
       step.fun="_gompertz_simulator"
       ),
     rmeasure="_gompertz_normal_rmeasure",
     dmeasure="_gompertz_normal_dmeasure",
     skeleton=map("_gompertz_skeleton",delta.t=1),
     paramnames=c("r","K","sigma","tau"),
     statenames=c("X"),
     fromEstimationScale=function(params,...){
       exp(params)
     },
     toEstimationScale=function(params,...){
       log(params)
     }
     ) -> gompertz

## the following was used to create the data included
## simulate(po,nsim=1,seed=299438676L) -> gompertz

c("gompertz")

====== ou2.R =======

library(pomp)

dat <- '"time";"y1";"y2"
1;-1,70081613181228;4,24746468495763
2;0,242170113354698;7,03941580692464
3;3,76758938785925;2,97520908716299
4;5,68588808149477;-1,47040054732921
5;3,77452601600709;-4,43380857201258
6;-3,08021312625511;-3,91887931755601
7;-2,62344860290795;-1,03202457662666
8;-2,86592128836057;4,35432873921529
9;-0,57760946681527;3,39758289565346
10;2,74097620676726;4,91071465265353
11;-0,95441253676346;10,1937011763325
12;-0,991464112428821;6,61823476151473
13;3,84297778313153;7,16145392465002
14;11,0789623105051;4,16444989466932
15;8,59819025333584;-4,90334695500249
16;2,19964347503517;-10,2511439107283
17;-1,13492397511025;-10,3624312929729
18;-2,56413579798821;-8,16885880550387
19;-4,91536849799727;-7,2203116145846
20;-4,33899642531645;-2,31022153289937
21;-1,35937286561288;2,99698279033669
22;-4,0242840268505;3,69116465405162
23;2,9387644957218;6,15414586126603
24;4,31828658150857;5,69290844809883
25;3,53319543208941;3,46054629112961
26;6,44281783433088;-4,79748856413256
27;0,718035321052913;-7,62976005316743
28;1,28027528075623;-6,25270934959545
29;3,70910008919496;-9,94125507650066
30;1,48104366840107;-10,2850854614516
31;-0,209117126395869;-8,80046618352614
32;-5,81639714549404;-13,2283385097361
33;-10,0225182259708;-6,4113389217938
34;-6,89506104214256;-0,998400718547352
35;-6,08611873904796;1,29299352993825
36;-7,74305796307499;5,90263448168222
37;0,0287081408447805;5,88274335246936
38;9,62512958922631;2,54860024892302
39;4,61772012246597;-4,51699357410536
40;-1,40028864684968;-3,66192759013006
41;-6,10811040670547;-4,38003912739696
42;-0,658058036920092;-1,7236046423518
43;-3,49326574531689;-4,14959269754278
44;2,54826871576915;-1,40873093710545
45;-0,640114486048052;0,896862978414845
46;-1,67566665036179;0,379384529694427
47;0,488141779118929;2,07741847310818
48;0,975197260459797;2,83822178936221
49;-1,0207098675894;0,112186038274661
50;0,0342598900248652;3,50605179482128
51;-2,39842556182189;2,92854885312383
52;5,19015239107048;-0,0089221935515531
53;4,14310779510487;0,268899079252409
54;2,63260618646816;1,60371893328454
55;-1,55607266713321;-0,764379706108408
56;-0,545766983243217;-0,837246196106889
57;1,34912551213756;-4,15285169037281
58;0,970330264719094;0,379978231358855
59;-5,04652345030695;-0,444477177942712
60;-3,97209430320528;1,10490455647656
61;-3,20325965055922;5,69632695844117
62;1,03436272221434;7,86129746544832
63;0,667898750837872;10,4068175981576
64;3,50430639591665;9,2816749036714
65;0,925079731276049;5,05429708955105
66;1,46558192088153;0,403851869094214
67;2,69317536547979;-4,51756265008713
68;-2,41608073604264;-8,3370959564069
69;-4,53015031919341;-3,20513929382432
70;-11,0133845324347;-1,33315487824171
71;-10,7508370993112;5,17621264446314
72;-10,5218612264298;8,61966772995883
73;-8,57143649619155;11,4982952428018
74;-3,99562214384817;12,5937325399812
75;-2,71455145450133;15,5796601067951
76;0,512617531996642;16,8599826830534
77;8,77766626349685;11,8289064956647
78;9,92868196527894;5,40891154021026
79;4,66011726941136;1,15575075091348
80;2,67821752937628;-0,585565498904033
81;1,15068730079675;-5,26069159041533
82;-2,203829735812;-7,73574749590483
83;-1,36113119160702;-6,51696428184672
84;-4,49088592785424;-6,88875950760792
85;-6,95435086585582;-5,73210381080353
86;-7,54286719246903;-2,04743677554387
87;-7,25975156626116;-1,10509448997233
88;-5,31555444374504;2,94363686064014
89;1,9728018375925;4,08030586218171
90;2,87436582964618;1,2087743319357
91;2,96600944673318;1,85680168956727
92;4,63179378892858;3,07085679112279
93;5,62360583343836;-2,89483235988436
94;1,87661311794902;-3,00789937429844
95;-0,0144619140532691;-2,67204400181558
96;1,94214418106964;-6,14711561527305
97;1,98468512734267;-5,61898548532659
98;-1,41012196458721;-11,086585198673
99;-7,29715227580992;-4,07855448780968
100;-4,81538906352058;-1,00365875150418
'

pomp( 
     data=read.csv2(text=dat),
     times="time",
     t0=0,
     rprocess=discrete.time.sim("_ou2_step",PACKAGE="pomp"),
     dprocess=onestep.dens("_ou2_pdf",PACKAGE="pomp"),
     dmeasure = "_ou2_dmeasure",
     rmeasure = "_ou2_rmeasure",
     skeleton = map("_ou2_skel",delta.t=1),
     PACKAGE="pomp",
     paramnames = c(
       "alpha.1","alpha.2","alpha.3","alpha.4",
       "sigma.1","sigma.2","sigma.3",
       "tau"
       ),
     statenames = c("x1","x2"),
     params=c(
       alpha.1=0.8, alpha.2=-0.5, alpha.3=0.3, alpha.4=0.9,
       sigma.1=3, sigma.2=-0.5, sigma.3=2,
       tau=1, 
       x1.0=-3, x2.0=4
       )
     ) -> ou2

c("ou2")

====== ricker.R =======

library(pomp)

dat <- '"time";"y"
0;68
1;2
2;87
3;0
4;12
5;174
6;0
7;0
8;1
9;57
10;11
11;178
12;0
13;1
14;0
15;34
16;72
17;3
18;101
19;0
20;8
21;156
22;0
23;0
24;3
25;93
26;0
27;17
28;121
29;0
30;0
31;19
32;107
33;0
34;4
35;127
36;0
37;1
38;47
39;8
40;117
41;0
42;3
43;82
44;2
45;39
46;70
47;11
48;275
49;0
50;0
'

pomp(
     data=read.csv2(text=dat),
     times="time",
     t0=0,
     params=c(r=exp(3.8),sigma=0.3,phi=10,c=1,N.0=7,e.0=0), # originally used to generate the data
     rprocess=discrete.time.sim(
       step.fun="_ricker_simulator"
       ),
     rmeasure="_ricker_poisson_rmeasure",
     dmeasure="_ricker_poisson_dmeasure",
     skeleton=map("_ricker_skeleton",delta.t=1),
     paramnames=c("r","sigma","phi","c"),
     statenames=c("N","e"),
     toEstimationScale=function(params,...) {
       params[c("r","sigma","phi","c","N.0")] <- log(params[c("r","sigma","phi","c","N.0")])
       params
     },
     fromEstimationScale=function(params,...) {
       params[c("r","sigma","phi","c","N.0")] <- exp(params[c("r","sigma","phi","c","N.0")])
       params
     },
     PACKAGE="pomp"
     ) -> ricker

c("ricker")

====== rw2.R =======

library(pomp)

dat <- '"time";"y1";"y2"
1;-1,88442138814419;1,16694679794963
2;0,397345132104944;2,33332499487602
3;-4,43903823096924;3,98978800111369
4;-3,07795898556872;6,68358950652736
5;-2,5084377724861;4,37016589397865
6;-2,66457136038371;10,2825383170911
7;-1,44355492377755;13,9831162784522
8;0,968437851613671;12,9665236550447
9;1,57185098023574;12,7659614606791
10;1,88961638094342;10,0954449344935
11;-1,12051626823462;10,7719406327611
12;0,162635095977655;15,025468449779
13;2,83442380381637;14,726164694454
14;0,847618205728183;13,5546613154093
15;0,105441944282331;11,7850550087575
16;0,296461065693377;16,0964275588307
17;0,808351449922765;12,7573633795118
18;0,180557522256085;14,1494260918142
19;-1,54720386433481;12,640929745979
20;-1,36317527995832;16,1564652150424
21;-2,14975241338301;14,0862494453931
22;-1,73362387899149;9,57894436683076
23;-4,11641486767628;12,3396158954935
24;-4,79684528795191;8,98019410182055
25;-4,65677411351422;4,53773975598888
26;-3,37927023673216;0,778920969568143
27;-2,49299858289081;2,98210668216634
28;-5,07199979564729;0,568214561291144
29;-8,14846233492297;12,4192780139244
30;-8,34021579959152;10,7851893298924
31;-8,0554839106157;8,38050019430963
32;-7,49140128186129;9,13845554655233
33;-9,94560640628807;11,8836891121119
34;-9,45170675671395;8,17857284326955
35;-8,4441895939124;5,9109331963516
36;-9,45982753677811;8,86602299854314
37;-9,24437347396321;16,0305936657508
38;-10,1336328971692;16,0647075572491
39;-11,9028131037299;18,7933411823485
40;-9,72801079504317;15,4101385644413
41;-11,4973632506713;23,630492310884
42;-8,98750967738246;20,5435130645714
43;-8,04940313116889;17,0322851685699
44;-8,45016514464591;15,6776626603499
45;-7,86249185601368;17,4261671143166
46;-9,1108152257913;13,2514704672741
47;-5,74359776629995;15,0391723729132
48;-9,45085485740965;9,84383972567
49;-9,43356939335136;8,58325220775032
50;-8,97952290180562;6,20169409394793
51;-8,44595737562255;5,71784567554888
52;-8,97118172249087;9,97933250202793
53;-8,45151319327985;9,11680772778525
54;-5,26131264548228;15,8933416924517
55;-5,16430864323134;14,170773142009
56;-6,78932972657802;13,4267139677724
57;-3,56789727727952;8,14108380116751
58;-6,74287381245534;11,0166048628484
59;-4,0245289198821;9,58793090014379
60;-5,03772395574636;9,68903665265393
61;-5,77457667405596;9,55816973802239
62;-5,0333451779765;8,98748477219307
63;-5,19083971746316;4,23538959486914
64;-3,96919059964988;1,38887983713013
65;-7,13867015535884;1,11548035325103
66;-4,6845514647567;5,23516815549156
67;-7,50412917294511;6,9527660027643
68;-7,50935496777481;9,44828450753692
69;-10,7888354929516;10,578008043056
70;-12,5924363121098;4,52419273323586
71;-12,4859985294564;-5,47179920501655
72;-11,6652687012524;-11,0285638166536
73;-12,0804923414833;-15,5180994963905
74;-10,8417578661395;-12,2110565771484
75;-9,05271883434386;-14,0944222890911
76;-10,9906152262749;-14,0225161881394
77;-9,19521310176121;-13,9652127577178
78;-10,6418892515783;-9,50194644402169
79;-10,0967357596072;-14,4818357261641
80;-10,2251248195123;-14,213435762591
81;-9,73893509065385;-4,48392863563174
82;-12,140122618462;-5,6261537871297
83;-10,915100326829;-1,34815640031926
84;-7,81183298776309;1,38933495434592
85;-11,9357645137211;-1,94337533715718
86;-10,4648437535608;-5,77201117370161
87;-10,0071269987997;-3,76206574032637
88;-12,7343781232615;-6,37608048073594
89;-7,27244655421996;-9,95764261287881
90;-10,4430703790289;-6,99969220125347
91;-10,4068713238478;-11,4302846737015
92;-9,26191379316831;-11,8452636790998
93;-8,54562417903865;-11,9494361360662
94;-10,3112275427129;-9,1309892339989
95;-12,7175391606616;-11,6769838280081
96;-9,59825698000055;-17,8102613321856
97;-10,3482501430257;-17,7997590083428
98;-12,2620571933296;-19,6661424778843
99;-10,2627440209765;-18,6163616927085
100;-10,1301733745858;-19,3339318994661
'

pomp(
     data=read.csv2(text=dat),
     times="time",
     t0=0,
     params=c(x1.0=0,x2.0=0,s1=1,s2=3,tau=1), # parameters at which data were generated
     rprocess = onestep.sim(
       step.fun = function(x, t, params, delta.t, ...) {
         c(
           x1=rnorm(n=1,mean=x['x1'],sd=params['s1']*delta.t),
           x2=rnorm(n=1,mean=x['x2'],sd=params['s2']*delta.t)
           )
       }
       ),
     dprocess = onestep.dens(
       dens.fun = function (x1, t1, x2, t2, params, ...) {
         sum(
             dnorm(
                   x=x2[c('x1','x2')],
                   mean=x1[c('x1','x2')],
                   sd=params[c('s1','s2')]*(t2-t1),
                   log=TRUE
                   ),
             na.rm=TRUE
             )
       }
       ),
     measurement.model=list(
       y1 ~ norm(mean=x1,sd=tau),
       y2 ~ norm(mean=x2,sd=tau)
       )
     ) -> rw2

c("rw2")

pomp documentation built on May 2, 2019, 4:09 p.m.