Nothing
inst/examples/blowflies.R
In pomp: Statistical Inference for Partially Observed Markov Processes
## 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
require(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"),
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)
}
) -> 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
)
test <- FALSE
if(test){
sim1 <- simulate(blowflies1,nsim=1)
plot(obs(sim1)['y',],ty='l')
lines(obs(blowflies1)['y',],lty="dashed")
states(sim1)[,1]
sim2 <- simulate(blowflies2,nsim=1)
plot(obs(sim2)['y',],ty='l')
lines(obs(blowflies2)['y',],lty="dashed")
states(sim2)[,1]
## check that it matches the deterministic skeleton when noise is small
params.1.skel <- coef(blowflies1)
params.1.skel["sigma.P"] <- 0.00001
params.1.skel["sigma.d"] <- 0.00001
params.1.skel["sigma.y"] <- 0.00001
simulate(blowflies1,params=params.1.skel,nsim=1,seed=73691676L) -> b1.skel
plot(obs(blowflies1)['y',],ty='l',lty="dashed")
lines(obs(b1.skel)['y',],ty='l')
}
c("blowflies1","blowflies2")
Try the pomp package in your browser
Any scripts or data that you put into this service are public.
pomp documentation built on May 2, 2019, 4:09 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
## 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
require(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"),
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)
}
) -> 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
)
test <- FALSE
if(test){
sim1 <- simulate(blowflies1,nsim=1)
plot(obs(sim1)['y',],ty='l')
lines(obs(blowflies1)['y',],lty="dashed")
states(sim1)[,1]
sim2 <- simulate(blowflies2,nsim=1)
plot(obs(sim2)['y',],ty='l')
lines(obs(blowflies2)['y',],lty="dashed")
states(sim2)[,1]
## check that it matches the deterministic skeleton when noise is small
params.1.skel <- coef(blowflies1)
params.1.skel["sigma.P"] <- 0.00001
params.1.skel["sigma.d"] <- 0.00001
params.1.skel["sigma.y"] <- 0.00001
simulate(blowflies1,params=params.1.skel,nsim=1,seed=73691676L) -> b1.skel
plot(obs(blowflies1)['y',],ty='l',lty="dashed")
lines(obs(b1.skel)['y',],ty='l')
}
c("blowflies1","blowflies2")
Try the pomp package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.