inst/examples/rw2.R

In pomp: Statistical Inference for Partially Observed Markov Processes

require(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")