[Pomp-commits] r800 - in pkg/pomp: . R inst inst/examples
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Jan 9 17:53:21 CET 2013
Author: kingaa
Date: 2013-01-09 17:53:21 +0100 (Wed, 09 Jan 2013)
New Revision: 800
Added:
pkg/pomp/R/example.R
pkg/pomp/inst/examples/blowflies.R
pkg/pomp/inst/examples/dacca.R
pkg/pomp/inst/examples/ewmeas.csv
pkg/pomp/inst/examples/gompertz.R
pkg/pomp/inst/examples/ou2.R
pkg/pomp/inst/examples/ricker.R
pkg/pomp/inst/examples/rw2.R
pkg/pomp/inst/examples/sir.R
pkg/pomp/inst/examples/verhulst.R
Removed:
pkg/pomp/inst/data-R/
Modified:
pkg/pomp/DESCRIPTION
pkg/pomp/NAMESPACE
Log:
- remove data()-loadable pomp objects in favor of the new 'pompExample' mechanism.
Modified: pkg/pomp/DESCRIPTION
===================================================================
--- pkg/pomp/DESCRIPTION 2013-01-09 16:05:51 UTC (rev 799)
+++ pkg/pomp/DESCRIPTION 2013-01-09 16:53:21 UTC (rev 800)
@@ -22,4 +22,4 @@
pmcmc.R pmcmc-methods.R compare-pmcmc.R
nlf-funcs.R nlf-guts.R nlf-objfun.R nlf.R
probe.R probe-match.R basic-probes.R spect.R spect-match.R
- builder.R
+ builder.R example.R
Modified: pkg/pomp/NAMESPACE
===================================================================
--- pkg/pomp/NAMESPACE 2013-01-09 16:05:51 UTC (rev 799)
+++ pkg/pomp/NAMESPACE 2013-01-09 16:53:21 UTC (rev 800)
@@ -92,5 +92,6 @@
sannbox,
spect.match,
traj.match.objfun,
- pompBuilder
+ pompBuilder,
+ pompExample
)
Added: pkg/pomp/R/example.R
===================================================================
--- pkg/pomp/R/example.R (rev 0)
+++ pkg/pomp/R/example.R 2013-01-09 16:53:21 UTC (rev 800)
@@ -0,0 +1,15 @@
+pompExample <- function (example) {
+ if (missing(example)) {
+ avlbl <- list.files(
+ path=system.file("examples",package="pomp"),
+ pattern=".+?R$"
+ )
+ avlbl <- gsub("\\.R$","",avlbl)
+ cat("available pomp examples:\n",sQuote(avlbl),"\n")
+ } else {
+ file <- system.file(file.path("examples",paste(example,".R",sep="")),package="pomp")
+ objs <- source(file,local=TRUE)
+ cat("newly created pomp objects:\n",objs$value,"\n")
+ }
+ invisible(NULL)
+}
Added: pkg/pomp/inst/examples/blowflies.R
===================================================================
--- pkg/pomp/inst/examples/blowflies.R (rev 0)
+++ pkg/pomp/inst/examples/blowflies.R 2013-01-09 16:53:21 UTC (rev 800)
@@ -0,0 +1,977 @@
+## 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(textConnection(blowfly.data),comment.char="#"),
+ set==4
+ )
+
+pomp(
+ data=subset(raw.data[c("day","y")],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"),
+ obsnames=c("y"),
+ 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),
+ parameter.inv.transform=function(params,...) {
+ log(params)
+ },
+ parameter.transform=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),
+ 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')
+
+}
+
+assign("blowflies1",blowflies1,envir=.GlobalEnv)
+assign("blowflies2",blowflies2,envir=.GlobalEnv)
+
+c("blowflies1","blowflies2")
Added: pkg/pomp/inst/examples/dacca.R
===================================================================
--- pkg/pomp/inst/examples/dacca.R (rev 0)
+++ pkg/pomp/inst/examples/dacca.R 2013-01-09 16:53:21 UTC (rev 800)
@@ -0,0 +1,159 @@
+## 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.
+
+require(pomp)
+
+mle <- c(
+ gamma=20.8,eps=19.1,rho=0,
+ delta=0.02, deltaI=0.06, clin=1, alpha=1,
+ beta.trend=-0.00498,
+ log.beta1=0.747, log.beta2=6.38, log.beta3=-3.44, log.beta4=4.23, log.beta5=3.33, log.beta6=4.55,
+ log.omega1=log(0.184), log.omega2=log(0.0786), log.omega3=log(0.0584), log.omega4=log(0.00917), log.omega5=log(0.000208), log.omega6=log(0.0124),
+ sd.beta=3.13, tau=0.23,
+ S.0=0.621, I.0=0.378, Rs.0=0, R1.0=0.000843, R2.0=0.000972, R3.0=1.16e-07,
+ nbasis=6, nrstage=3
+ )
+
+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
+nbasis <- as.integer(mle["nbasis"])
+nrstage <- as.integer(mle["nrstage"])
+
+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,
+ seas=periodic.bspline.basis(tcovar-1/12,nbasis=nbasis,degree=3,period=1),
+ 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)
+ )
+
+pomp(
+ data=cholera,
+ times='time',
+ t0=t0,
+ nrstage = nrstage,
+ rprocess = euler.sim(
+ step.fun = "_cholmodel_one",
+ PACKAGE="pomp",
+ delta.t=1/240
+ ),
+ dmeasure = "_cholmodel_norm_dmeasure",
+ rmeasure="_cholmodel_norm_rmeasure",
+ PACKAGE = "pomp",
+ covar=covartable,
+ tcovar='time',
+ zeronames = c("M","count"),
+ obsnames = "cholera.deaths",
+ statenames = c("S","I","Rs","R1","M","W","count"),
+ paramnames = c("tau","gamma","eps","delta","deltaI",
+ "log.omega1","sd.beta","beta.trend","log.beta1",
+ "alpha","rho","clin","nbasis","nrstage"),
+ covarnames = c("pop","dpopdt","seas.1","trend"),
+ all.state.names=c("S","I","Rs",paste("R",1:nrstage,sep=''),"M","W","count"),
+ comp.names=c("S","I","Rs",paste("R",1:nrstage,sep='')),
+ comp.ic.names=c("S.0","I.0","Rs.0",paste("R",1:nrstage,".0",sep='')),
+ log.trans=c( # parameters to log transform
+ "gamma","eps","rho","delta","deltaI","alpha",
+ "tau","sd.beta",
+ "S.0","I.0","Rs.0",paste("R",1:nrstage,".0",sep='')
+ ),
+ logit.trans="clin", # parameters to logit transform
+ parameter.transform=function (params, log.trans, logit.trans, comp.ic.names, ...) {
+ params[logit.trans] <- plogis(params[logit.trans])
+ params[log.trans] <- exp(params[log.trans])
+ params[comp.ic.names] <- params[comp.ic.names]/sum(params[comp.ic.names])
+ params
+ },
+ parameter.inv.transform=function (params, log.trans, logit.trans, comp.names, ...) {
+ params[logit.trans] <- qlogis(params[logit.trans])
+ params[log.trans] <- log(params[log.trans])
+ params
+ },
+ initializer = function (params, t0, covars, nrstage, comp.ic.names, comp.names, all.state.names, ...) {
+ states <- numeric(length(all.state.names))
+ names(states) <- all.state.names
+ ## translate fractions into initial conditions
+ frac <- params[comp.ic.names]
+ states[comp.names] <- round(covars['pop']*frac/sum(frac))
+ states
+ }
+ ) -> dacca
+
+coef(dacca) <- mle
+
+assign("dacca",dacca,envir=.GlobalEnv)
+c("dacca")
Added: pkg/pomp/inst/examples/ewmeas.csv
===================================================================
--- pkg/pomp/inst/examples/ewmeas.csv (rev 0)
+++ pkg/pomp/inst/examples/ewmeas.csv 2013-01-09 16:53:21 UTC (rev 800)
@@ -0,0 +1,992 @@
+"time";"reports"
+1948,0274;3762
+1948,0466;3301
+1948,0658;3454
+1948,0849;4222
+1948,1041;5604
+1948,1233;6218
+1948,1425;7216
+1948,1616;7370
+1948,1781;8828
+1948,1973;8521
+1948,2164;10056
+1948,2356;9365
+1948,2548;10901
+1948,274;11285
+1948,2932;9749
+1948,3123;9135
+1948,3315;10210
+1948,3507;10210
+1948,3699;11975
+1948,389;11745
+1948,4082;13511
+1948,4274;10901
+1948,4466;12283
+1948,4658;10670
+1948,4849;10670
+1948,5041;9058
+1948,5233;8905
+1948,5425;9519
+1948,5616;8521
+1948,5808;8291
+1948,6;6909
+1948,6192;5911
+1948,6384;4606
+1948,6575;3762
+1948,6767;2994
+1948,6959;2610
+1948,7151;2610
+1948,7342;2687
+1948,7534;3531
+1948,7726;4069
+1948,7918;4529
+1948,811;5297
+1948,8301;6218
+1948,8493;6218
+1948,8685;6986
+1948,8877;6986
+1948,9068;8598
+1948,926;8598
+1948,9452;10133
+1948,9644;9826
+1948,9836;8982
+1949,0027;11208
+1949,0219;13279
+1949,0411;11053
+1949,0603;10311
+1949,0795;11721
+1949,0986;13872
+1949,1178;18472
+1949,137;19584
+1949,1562;20104
+1949,1753;18694
+1949,1945;17359
+1949,2137;16320
+1949,2329;15356
+1949,2521;15875
+1949,2712;14762
+1949,2904;13130
+1949,3096;13575
+1949,3288;11498
+1949,3479;9273
+1949,3671;8383
+1949,3863;8531
+1949,4055;9644
+1949,4247;9792
+1949,4438;9273
+1949,463;9199
+1949,4822;6899
+1949,5014;6973
+1949,5205;5564
+1949,5397;5415
+1949,5589;5044
+1949,5781;4303
+1949,5973;3709
+1949,6164;3264
+1949,6356;2151
+1949,6548;1780
+1949,674;1187
+1949,6932;890
+1949,7123;668
+1949,7315;668
+1949,7507;816
+1949,7699;964
+1949,789;964
+1949,8082;1039
+1949,8274;1039
+1949,8466;1113
+1949,8658;1039
+1949,8849;1335
+1949,9041;1484
+1949,9233;1706
+1949,9425;1632
+1949,9616;1780
+1949,9808;1706
+1950;2003
+1950,0192;2993
+1950,0384;2993
+1950,0575;2993
+1950,0767;3265
+1950,0959;4535
+1950,1151;5442
+1950,1342;5804
+1950,1534;6439
+1950,1726;5986
+1950,1918;7074
+1950,211;1723
+1950,2301;7346
+1950,2493;8072
+1950,2685;7074
+1950,2877;10158
+1950,3068;8797
+1950,326;8162
+1950,3452;8072
+1950,3644;9341
+1950,3836;11337
+1950,4027;12062
+1950,4219;12969
+1950,4411;13785
+1950,4603;11337
+1950,4795;11699
+1950,4986;11427
+1950,5178;11518
+1950,537;11971
+1950,5562;11246
+1950,5753;10611
+1950,5945;10248
+1950,6137;8253
+1950,6329;7074
+1950,6521;5532
+1950,6712;4353
+1950,6904;3265
+1950,7096;3174
+1950,7288;3446
+1950,7479;3809
+1950,7671;4988
+1950,7863;6439
+1950,8055;7527
+1950,8247;8525
+1950,8438;10067
+1950,863;11337
+1950,8822;13060
+1950,9014;14601
+1950,9205;15780
+1950,9397;18864
+1950,9589;18320
+1950,9781;18501
+1950,9973;22401
+1951,0164;21682
+1951,0356;18720
+1951,0548;15995
+1951,074;15995
+1951,0932;16350
+1951,1123;21682
+1951,1315;24644
+1951,1507;27843
+1951,1699;29264
+1951,189;29146
+1951,2082;31871
+1951,2274;27843
+1951,2466;31278
+1951,2658;28790
+1951,2849;22629
+1951,3041;19075
+1951,3233;18009
+1951,3425;18364
+1951,3616;19904
+1951,3808;18483
+1951,4;19194
+1951,4192;15639
+1951,4384;14810
+1951,4575;12203
+1951,4767;11374
+1951,4959;10071
+1951,5151;8530
+1951,5342;7701
+1951,5534;6279
+1951,5726;5213
+1951,5918;4976
+1951,611;3673
+1951,6301;3199
+1951,6493;2370
+1951,6685;1896
+1951,6877;1422
+1951,7068;1185
+1951,726;1066
+1951,7452;1303
+1951,7644;1422
+1951,7836;1659
+1951,8027;1777
+1951,8219;2014
+1951,8411;2014
+1951,8603;1896
+1951,8795;2133
+1951,8986;2014
+1951,9178;2370
+1951,937;2370
+1951,9562;2488
+1951,9753;2370
+1951,9945;2370
+1952,0137;3445
+1952,0329;2846
+1952,0521;2771
+1952,0712;2846
+1952,0904;3370
+1952,1096;4419
+1952,1288;4943
+1952,1479;5692
+1952,1644;6141
+1952,1836;6815
+1952,2027;7040
+1952,2219;6815
+1952,2411;6890
+1952,2603;6740
+1952,2795;5392
+1952,2986;6740
+1952,3178;6815
+1952,337;4718
+1952,3562;4419
+1952,3753;5692
+1952,3945;6216
+1952,4137;6591
+1952,4329;6291
+1952,4521;7265
+1952,4712;6216
+1952,4904;6740
+1952,5096;6666
+1952,5288;7489
+1952,5479;7115
+1952,5671;7265
+1952,5863;7265
+1952,6055;6965
+1952,6247;6066
+1952,6438;5168
+1952,663;3595
+1952,6822;2921
+1952,7014;2397
+1952,7205;2546
+1952,7397;3220
+1952,7589;4643
+1952,7781;5692
+1952,7973;7714
+1952,8164;9062
+1952,8356;11084
+1952,8548;11159
+1952,874;11908
+1952,8932;11983
+1952,9123;14679
+1952,9315;16027
+1952,9507;17600
+1952,9699;19398
+1952,989;14080
+1953,0082;31125
+1953,0274;25257
+1953,0466;21798
+1953,0658;22322
+1953,0849;27143
+1953,1041;32068
+1953,1233;33012
+1953,1425;33536
+1953,1616;31544
+1953,1808;28191
+1953,2;27457
+1953,2192;27772
+1953,2384;26305
+1953,2575;21589
+1953,2767;23894
+1953,2959;19493
+1953,3151;14043
+1953,3342;11004
+1953,3534;9956
+1953,3726;10270
+1953,3918;10480
+1953,411;9222
+1953,4301;8908
+1953,4493;7650
+1953,4685;7755
+1953,4877;7022
+1953,5068;6707
+1953,526;5659
+1953,5452;4926
+1953,5644;4506
+1953,5836;3668
+1953,6027;2830
+1953,6219;2201
+1953,6411;1572
+1953,6603;1048
+1953,6795;838
+1953,6986;629
+1953,7178;629
+1953,737;524
+1953,7562;629
+1953,7753;734
+1953,7945;734
+1953,8137;838
+1953,8329;838
+1953,8521;1048
+1953,8712;838
+1953,8904;838
+1953,9096;943
+1953,9288;943
+1953,9479;943
+1953,9671;838
+1953,9863;629
+1954,0055;1362
+1954,0247;1131
+1954,0438;1159
+1954,063;1102
+1954,0822;1187
+1954,1014;1498
+1954,1205;1837
+1954,1397;1950
+1954,1589;2092
+1954,1781;2431
+1954,1973;2205
+1954,2164;2318
+1954,2356;2063
+1954,2548;2318
+1954,274;2289
+1954,2932;1922
+1954,3123;2205
+1954,3315;1978
+1954,3507;1724
+1954,3699;1696
+1954,389;1950
+1954,4082;2035
+1954,4274;2826
+1954,4466;2798
+1954,4658;3109
+1954,4849;2826
+1954,5041;3081
+1954,5233;2685
+1954,5425;3081
+1954,5616;2968
+1954,5808;3392
+1954,6;2883
+1954,6192;2826
+1954,6384;2261
+1954,6575;1696
+1954,6767;1583
+1954,6959;1074
+1954,7151;1272
+1954,7342;1413
+1954,7534;1837
+1954,7726;2346
+1954,7918;2798
+1954,811;3448
+1954,8301;3787
+1954,8493;4635
+1954,8685;4211
+1954,8877;4720
+1954,9068;5285
+1954,926;5992
+1954,9452;6275
+1954,9644;6586
+1954,9836;5822
+1955,0027;8310
+1955,0219;9339
+1955,0411;7871
+1955,0603;7071
+1955,0795;10673
+1955,0986;12807
+1955,1178;18411
+1955,137;20012
+1955,1562;23747
+1955,1753;26549
+1955,1945;25481
+1955,2137;29217
+1955,2329;28683
+1955,2521;32419
+1955,2712;26949
+1955,2904;33086
+1955,3096;28683
+1955,3288;22546
+1955,3479;19478
+1955,3671;20412
+1955,3863;23747
+1955,4055;26148
+1955,4247;24814
+1955,4438;27749
+1955,463;21879
+1955,4822;20545
+1955,5014;17744
+1955,5205;18411
+1955,5397;16409
+1955,5589;14542
+1955,5781;12140
+1955,5973;10006
+1955,6164;8138
+1955,6356;5603
+1955,6548;3735
+1955,674;2668
+1955,6932;1601
+1955,7123;1067
+1955,7315;1067
+1955,7507;1201
+1955,7699;1468
+1955,789;1734
+1955,8082;2135
+1955,8274;2268
+1955,8466;2268
+1955,8658;2401
+1955,8849;2401
+1955,9041;2535
+1955,9233;2802
+1955,9425;2802
+1955,9616;2935
+1955,9808;2535
+1956;2802
+1956,0192;3272
+1956,0384;2407
+1956,0575;1975
+1956,0767;2161
+1956,0959;2469
+1956,1151;2932
+1956,1342;3179
+1956,1534;3210
+1956,1699;2840
+1956,189;2685
+1956,2082;2932
+1956,2274;2840
+1956,2466;2716
+1956,2658;3395
+1956,2849;3025
+1956,3041;2130
+1956,3233;2037
+1956,3425;2191
+1956,3616;3117
+1956,3808;3457
+1956,4;3272
+1956,4192;3889
+1956,4384;3241
+1956,4575;3426
+1956,4767;3272
+1956,4959;3796
+1956,5151;3488
+1956,5342;3333
+1956,5534;3333
+1956,5726;3333
+1956,5918;3272
+1956,611;2716
+1956,6301;2438
+1956,6493;1883
+1956,6685;1389
+1956,6877;1111
+1956,7068;957
+1956,726;1111
+1956,7452;1420
+1956,7644;2006
+1956,7836;2161
+1956,8027;2840
+1956,8219;3056
+1956,8411;3117
+1956,8603;3148
+1956,8795;3549
+1956,8986;4105
+1956,9178;4568
+1956,937;6050
+1956,9562;6574
+1956,9753;6729
+1956,9945;6975
+1957,0137;10966
+1957,0329;9260
+1957,0521;8164
+1957,0712;9504
+1957,0904;13037
+1957,1096;16814
+1957,1288;19739
+1957,1479;21323
+1957,1671;23029
+1957,1863;24856
+1957,2055;24856
+1957,2247;27171
+1957,2438;25222
+1957,263;25344
+1957,2822;27293
+1957,3014;21810
+1957,3205;27049
+1957,3397;23150
+1957,3589;18277
+1957,3781;16693
+1957,3973;19495
+1957,4164;21688
+1957,4356;23394
+1957,4548;21810
+1957,474;22297
+1957,4932;17424
+1957,5123;15474
+1957,5315;13403
+1957,5507;12794
+1957,5699;11575
+1957,589;9869
+1957,6082;7920
+1957,6274;6580
+1957,6466;4265
+1957,6658;3046
+1957,6849;2071
+1957,7041;1218
+1957,7233;1097
+1957,7425;975
+1957,7616;1340
+1957,7808;1218
+1957,8;1340
+1957,8192;1462
+1957,8384;1706
+1957,8575;1706
+1957,8767;1706
+1957,8959;1950
+1957,9151;1828
+1957,9342;2193
+1957,9534;2437
+1957,9726;2437
+1957,9918;2071
+1958,011;3660
+1958,0301;3001
+1958,0493;2538
+1958,0685;2428
+1958,0877;3608
+1958,1068;3898
+1958,126;4709
+1958,1452;4798
+1958,1644;4554
+1958,1836;4709
+1958,2027;4168
+1958,2219;4812
+1958,2411;4456
+1958,2603;4359
+1958,2795;5265
+1958,2986;5294
+1958,3178;4074
+1958,337;3843
+1958,3562;3717
+1958,3753;4736
+1958,3945;5494
+1958,4137;5747
+1958,4329;6677
+1958,4521;5804
+1958,4712;6362
+1958,4904;6257
+1958,5096;6696
+1958,5288;7356
+1958,5479;6934
+1958,5671;6850
+1958,5863;6554
+1958,6055;5824
+1958,6247;5301
+1958,6438;3944
+1958,663;2906
+1958,6822;2056
+1958,7014;1669
+1958,7205;1529
+1958,7397;1688
+1958,7589;2066
+1958,7781;2888
+1958,7973;3886
+1958,8164;4394
+1958,8356;5209
+1958,8548;5525
+1958,874;6019
+1958,8932;6352
+1958,9123;7183
+1958,9315;8723
+1958,9507;9814
+1958,9699;10539
+1958,989;8368
+1959,0082;17815
+1959,0274;14841
+1959,0466;12422
+1959,0658;13435
+1959,0849;17703
+1959,1041;21782
+1959,1233;24696
+1959,1425;25932
+1959,1616;24506
+1959,1808;23081
+1959,2;22336
+1959,2192;22897
+1959,2384;19163
+1959,2575;24763
+1959,2767;22426
+1959,2959;17344
+1959,3151;14185
+1959,3342;14046
+1959,3534;15422
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/pomp -r 800
More information about the pomp-commits
mailing list