[Pomp-commits] r457 - in pkg: . data inst inst/data-R inst/doc man src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu May 5 18:32:20 CEST 2011
Author: kingaa
Date: 2011-05-05 18:32:19 +0200 (Thu, 05 May 2011)
New Revision: 457
Added:
pkg/data/blowflies.rda
pkg/inst/data-R/blowflies.R
pkg/inst/data-R/blowflies.csv
pkg/man/blowflies.Rd
pkg/src/tsir.c
Modified:
pkg/DESCRIPTION
pkg/data/dacca.rda
pkg/inst/ChangeLog
pkg/inst/data-R/dacca.R
pkg/inst/doc/advanced_topics_in_pomp.pdf
pkg/inst/doc/intro_to_pomp.pdf
pkg/man/dacca.Rd
Log:
- add blowflies example
- fix bug in tsir code
- rearrange 'dacca' example
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2011-05-04 21:33:27 UTC (rev 456)
+++ pkg/DESCRIPTION 2011-05-05 16:32:19 UTC (rev 457)
@@ -2,7 +2,7 @@
Type: Package
Title: Statistical inference for partially observed Markov processes
Version: 0.36-7
-Date: 2011-05-04
+Date: 2011-05-05
Author: Aaron A. King, Edward L. Ionides, Carles Breto, Steve Ellner, Bruce Kendall, Helen Wearing,
Matthew J. Ferrari, Michael Lavine, Daniel C. Reuman
Maintainer: Aaron A. King <kingaa at umich.edu>
Added: pkg/data/blowflies.rda
===================================================================
(Binary files differ)
Property changes on: pkg/data/blowflies.rda
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Modified: pkg/data/dacca.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/ChangeLog
===================================================================
--- pkg/inst/ChangeLog 2011-05-04 21:33:27 UTC (rev 456)
+++ pkg/inst/ChangeLog 2011-05-05 16:32:19 UTC (rev 457)
@@ -1,3 +1,22 @@
+2011-05-04 kingaa
+
+ * [r456] R/aaa.R, R/mif-methods.R, R/pfilter-methods.R,
+ R/pfilter.R, man/mif-methods.Rd, man/pfilter-methods.Rd,
+ man/pfilter.Rd: - 'pred.mean', 'pred.var', and 'filter.mean'
+ methods are now defined for 'pfilterd.pomp' objects (and not just
+ 'mif' objects)
+ * [r455] R/trajectory-pomp.R: - 'trajectory' now gives a more
+ informative error when no skeleton is present
+ * [r454] DESCRIPTION, data/dacca.rda, data/euler.sir.rda,
+ data/gillespie.sir.rda, data/gompertz.rda, data/ou2.rda,
+ data/ricker.rda, data/rw2.rda, data/verhulst.rda, inst/ChangeLog,
+ inst/NEWS, inst/data-R/dacca.R, inst/data-R/euler.sir.R,
+ inst/data-R/gillespie.sir.R, inst/data-R/gompertz.R,
+ inst/data-R/ou2.R, inst/data-R/ricker.R, inst/data-R/rw2.R,
+ inst/data-R/verhulst.R, src/blowfly.c, src/probe_acf.c: - fix bug
+ in error message in 'probe_acf.c'
+ - update data()-loadable examples
+
2011-05-03 kingaa
* [r453] inst/doc/intro_to_pomp.pdf:
Added: pkg/inst/data-R/blowflies.R
===================================================================
--- pkg/inst/data-R/blowflies.R (rev 0)
+++ pkg/inst/data-R/blowflies.R 2011-05-05 16:32:19 UTC (rev 457)
@@ -0,0 +1,121 @@
+## 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
+
+raw.data <- subset(
+ read.csv2("blowflies.csv",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_model_simulator",
+ delta.t=1,
+ ),
+ paramnames=c("log.P","log.N0","log.delta","log.sigma.P","log.sigma.d","tau","log.sigma.y"),
+ statenames=c("N1","R","S","e","eps"),
+ obsnames=c("y"),
+ measurement.model=y~nbinom(mu=N1,size=exp(-2*log.sigma.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),
+ 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(
+ data=subset(raw.data[c("day","y")],day>14&day<400),
+ times="day",
+ t0=14,
+ rprocess=discrete.time.sim(
+ step.fun="_blowfly_model_simulator",
+ delta.t=2,
+ ),
+ paramnames=c("log.P","log.N0","log.delta","log.sigma.P","log.sigma.d","tau","log.sigma.y"),
+ statenames=c("N1","R","S","e","eps"),
+ obsnames=c("y"),
+ measurement.model=y~nbinom(mu=N1,size=exp(-2*log.sigma.y)),
+ 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) <- c(
+ log.P = 1.189 ,
+ log.delta = -1.828 ,
+ log.N0 = 6.522 ,
+ log.sigma.P = 0.301 ,
+ log.sigma.d = -0.292 ,
+ log.sigma.y = -3.625 ,
+ tau = 14
+ )
+
+## mle from search to date
+coef(blowflies2) <- c(
+ log.P = 1.005 ,
+ log.delta = -1.75 ,
+ log.N0 = 6.685 ,
+ log.sigma.P = 0.366 ,
+ log.sigma.d = -0.274 ,
+ log.sigma.y = -4.524 ,
+ tau = 7
+ )
+
+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["log.sigma.P"] <- log(0.00001)
+ params.1.skel["log.sigma.d"] <- log(0.00001)
+ params.1.skel["log.sigma.y"] <- log(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')
+
+}
+
+save(blowflies1,blowflies2,file="blowflies.rda")
Added: pkg/inst/data-R/blowflies.csv
===================================================================
--- pkg/inst/data-R/blowflies.csv (rev 0)
+++ pkg/inst/data-R/blowflies.csv 2011-05-05 16:32:19 UTC (rev 457)
@@ -0,0 +1,859 @@
+"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
Modified: pkg/inst/data-R/dacca.R
===================================================================
--- pkg/inst/data-R/dacca.R 2011-05-04 21:33:27 UTC (rev 456)
+++ pkg/inst/data-R/dacca.R 2011-05-05 16:32:19 UTC (rev 457)
@@ -101,8 +101,6 @@
trend=tcovar-mean(tcovar)
)
-dacca <- list()
-
pomp(
data=cholera,
times='time',
@@ -155,13 +153,13 @@
states[comp.names] <- round(covars['pop']*frac/sum(frac))
states
}
- ) -> dacca$po
+ ) -> dacca
## Parameter transformations
## Parameters are fit in the transformed space.
## Positive parameters are log transformed, probabilities are logit transformed.
## Initial conditions are parameterized in terms of log-transformed fractions in each compartment.
-dacca$transform <- function (params, dir = c("forward","inverse")) {
+dacca.transform <- function (params, dir = c("forward","inverse")) {
dir <- match.arg(dir)
r <- length(dim(params))
nm <- if (r>0) rownames(params) else names(params)
@@ -191,5 +189,5 @@
)
}
-coef(dacca$po) <- dacca$transform(mle,dir="forward")
-save(dacca,file="dacca.rda")
+coef(dacca) <- dacca.transform(mle,dir="forward")
+save(dacca,dacca.transform,file="dacca.rda")
Modified: pkg/inst/doc/advanced_topics_in_pomp.pdf
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/intro_to_pomp.pdf
===================================================================
(Binary files differ)
Added: pkg/man/blowflies.Rd
===================================================================
--- pkg/man/blowflies.Rd (rev 0)
+++ pkg/man/blowflies.Rd 2011-05-05 16:32:19 UTC (rev 457)
@@ -0,0 +1,17 @@
+\name{blowflies}
+\alias{blowflies}
+\alias{blowflies1}
+\alias{blowflies2}
+\docType{data}
+\title{Model for Nicholson's blowflies.}
+\description{
+ \code{blowfly1} and \code{blowfly2} are \code{pomp} objects encoding stochastic delay-difference models.
+}
+\usage{data(blowflies)}
+\examples{
+data(blowflies)
+plot(blowflies1)
+plot(blowflies2)
+}
+\seealso{\code{\link{pomp-class}} and the vignettes}
+\keyword{datasets}
Modified: pkg/man/dacca.Rd
===================================================================
--- pkg/man/dacca.Rd 2011-05-04 21:33:27 UTC (rev 456)
+++ pkg/man/dacca.Rd 2011-05-05 16:32:19 UTC (rev 457)
@@ -1,9 +1,10 @@
\name{dacca}
\alias{dacca}
+\alias{dacca.transform}
\docType{data}
\title{Model of cholera transmission for historic Bengal.}
\description{
- \code{dacca} contains census and cholera mortality data from the Dacca district of the former British province of Bengal over the years 1891 to 1940 together with a stochastic differential equation transmission model.
+ \code{dacca} is a \code{pomp} object containing census and cholera mortality data from the Dacca district of the former British province of Bengal over the years 1891 to 1940 together with a stochastic differential equation transmission model.
The model is that of King et al. (2008).
It also has the MLE for the SIRS model with seasonal reservoir.
@@ -13,23 +14,16 @@
data(dacca)
}
\details{
- \code{dacca} is a list consisting of
- \describe{
- \item{po}{
- A \code{pomp} object containing the model, data, and MLE parameters.
- }
- \item{transform}{
- A function for transforming the parameters.
- }
- }
+ \code{dacca} is a \code{pomp} object containing the model, data, and MLE parameters.
+ \code{dacca.transform} is a function that transforms the parameters.
}
\examples{
data(dacca)
-plot(dacca$po)
-theta <- dacca$transform(coef(dacca$po),dir="inverse") #MLEs
-plot(simulate(dacca$po))
+plot(dacca)
+theta <- dacca.transform(coef(dacca),dir="inverse") #MLEs
+plot(simulate(dacca))
theta["eps"] <- 1
-plot(simulate(dacca$po,params=dacca$transform(theta,dir="forward")))
+plot(simulate(dacca,params=dacca.transform(theta,dir="forward")))
}
\references{
King, A. A., Ionides, E. L., Pascual, M., and Bouma, M. J.
Added: pkg/src/tsir.c
===================================================================
--- pkg/src/tsir.c (rev 0)
+++ pkg/src/tsir.c 2011-05-05 16:32:19 UTC (rev 457)
@@ -0,0 +1,118 @@
+// dear emacs, please treat this as -*- C++ -*-
+
+#include <Rmath.h>
+
+#include "pomp.h"
+
+#define LOG_BETA (p[parindex[0]]) // transmission rates
+#define LOG_M (p[parindex[1]]) // Poisson import rate
+#define LOG_ALPHA (p[parindex[2]]) // mixing exponent
+#define LOG_RHO (p[parindex[3]]) // under-reporting
+#define PERIOD (p[parindex[4]]) // period of seasonality
+#define DEGREE (p[parindex[5]]) // degree of B-splines
+#define NBASIS (p[parindex[6]]) // number of B-spline basis functions
+#define LOG_SIGMA (p[parindex[7]]) // noise intensity
+#define LOG_OD (p[parindex[8]]) // measurement overdispersion parameter
+
+#define BIRTHS (covar[covindex[0]]) // numbers of births
+
+#define SUS (x[stateindex[0]]) // susceptibles
+#define INF (x[stateindex[1]]) // infectives
+#define THETA (x[stateindex[2]]) // imports
+
+#define REPORTS (y[obsindex[0]]) // reported cases
+
+#define NEWSUS (f[stateindex[0]]) // susceptibles
+#define NEWINF (f[stateindex[1]]) // infectives
+#define NEWTHETA (f[stateindex[2]]) // imports
+
+void _tsir_dmeasure (double *lik, double *y, double *x, double *p, int give_log,
+ int *obsindex, int *stateindex, int *parindex, int *covindex,
+ int ncovars, double *covars, double t)
+{
+ double mu, size, prob;
+ size = exp(-LOG_OD); // od = 1/size
+ mu = exp(LOG_RHO)*INF; // mean
+ prob = size/(size+mu);
+ // *lik = dbinom(REPORTS,nearbyint(INF),exp(LOG_RHO),give_log);
+ *lik = dnbinom(REPORTS,size,prob,give_log);
+}
+
+void _tsir_rmeasure (double *y, double *x, double *p,
+ int *obsindex, int *stateindex, int *parindex, int *covindex,
+ int ncovars, double *covars, double t) {
+ double mu, size, prob;
+ size = exp(-LOG_OD); // od = 1/size
+ mu = exp(LOG_RHO)*INF; // mean
+ prob = size/(size+mu);
+ // REPORTS = rbinom(nearbyint(INF),exp(LOG_RHO));
+ REPORTS = rnbinom(size,prob);
+}
+
+// Ricker model with log-normal process noise
+void _tsir_simulator (double *x, const double *p,
+ const int *stateindex, const int *parindex, const int *covindex,
+ int covdim, const double *covar,
+ double t, double dt)
+{
+ double beta;
+ double em;
+ double alpha;
+ double sigma, eps;
+ double prob, lambda;
+ int nbasis;
+ int degree;
+
+ nbasis = (int) NBASIS;
+ degree = (int) DEGREE;
+
+ em = exp(LOG_M);
+ alpha = exp(LOG_ALPHA);
+ sigma = exp(LOG_SIGMA);
+
+ {
+ double seasonality[nbasis];
+ periodic_bspline_basis_eval(t,PERIOD,degree,nbasis,&seasonality[0]);
+ beta = exp(dot_product(nbasis,&seasonality[0],&LOG_BETA));
+ }
+
+ if (sigma!=0)
+ eps = exp(rnorm(0,sigma));
+ else
+ eps = 1;
+
+ THETA = rpois(em); // imports
+ lambda = beta*pow(INF+THETA,alpha)*eps; // force of infection
+ INF = rbinom(nearbyint(SUS),1-exp(-lambda));
+ SUS += BIRTHS-INF; // susceptible balance
+}
+
+void _tsir_skeleton (double *f, double *x, const double *p,
+ const int *stateindex, const int *parindex, const int *covindex,
+ int covdim, const double *covar, double t)
+{
+ double beta;
+ double em;
+ double alpha;
+ double lambda;
+ int nbasis;
+ int degree;
+
+ nbasis = (int) NBASIS;
+ degree = (int) DEGREE;
+
+ em = exp(LOG_M);
+ alpha = exp(LOG_ALPHA);
+
+ {
+ double seasonality[nbasis];
+ periodic_bspline_basis_eval(t,PERIOD,degree,nbasis,&seasonality[0]);
+ beta = exp(dot_product(nbasis,&seasonality[0],&LOG_BETA));
+ }
+
+ lambda = beta*pow(INF+em,alpha); // force of infection
+ NEWINF = SUS*(1-exp(-lambda));
+ NEWSUS = SUS+BIRTHS-NEWINF; // susceptible balance
+ NEWTHETA = em; // imports
+
+}
Property changes on: pkg/src/tsir.c
___________________________________________________________________
Added: svn:eol-style
+ native
More information about the pomp-commits
mailing list