[Pomp-commits] r15 - pkg/inst/examples

[noreply at r-forge.r-project.org]

Author: kingaa
Date: 2008-07-25 12:02:49 +0200 (Fri, 25 Jul 2008)
New Revision: 15

Modified:
   pkg/inst/examples/sir.R
Log:
fix up R-function version

Modified: pkg/inst/examples/sir.R
===================================================================
--- pkg/inst/examples/sir.R	2008-07-25 09:49:54 UTC (rev 14)
+++ pkg/inst/examples/sir.R	2008-07-25 10:02:49 UTC (rev 15)
@@ -27,18 +27,13 @@
 
 ## set up the pomp object
 po <- pomp(
-           times=seq(1/52,4,by=1/52),
+           times=1/52*seq.int(length=4*52),
            data=rbind(measles=numeric(52*4)),
            t0=0,
            tcovar=tbasis,
            covar=basis,
            delta.t=1/52/20,
-           statenames=c("S","I","R","cases","W","B","dW"),
-           paramnames=c("gamma","mu","iota","beta1","beta.sd","pop"),
-           covarnames=c("seas1"),
            zeronames=c("cases"),
-           measurement.model=measles~binom(size=cases,prob=exp(rho)),
-           rprocess=euler.simulate,
            step.fun=function(t,x,params,covars,delta.t,...) {
              params <- exp(params)
              with(
@@ -55,18 +50,22 @@
                                reulermultinom(n=1,size=R,rate=c(mu),dt=delta.t)
                                )
                     c(
-                      S+trans[1]-trans[2]-trans[3],
-                      I+trans[2]-trans[4]-trans[5],
-                      R+trans[4]-trans[6],
-                      cases+trans[4],
-                      if (beta.sd>0) W+(dW-delta.t)/beta.sd else W,
-                      trans,
-                      dW
+                      S=S+trans[1]-trans[2]-trans[3],
+                      I=I+trans[2]-trans[4]-trans[5],
+                      R=R+trans[4]-trans[6],
+                      cases=cases+trans[4],
+                      W=if (beta.sd>0) W+(dW-delta.t)/beta.sd else W,
+                      B=trans[1],
+                      SI=trans[2],
+                      SD=trans[3],
+                      IR=trans[4],
+                      ID=trans[5],
+                      RD=trans[6],
+                      dW=dW
                       )
                   }
                   )
            },
-           dprocess=euler.density,
            dens.fun=function(t1,t2,params,x1,x2,covars,...) {
              params <- exp(params)
              with(
@@ -112,6 +111,9 @@
                   }
                   )
            },
+           rprocess=euler.simulate,
+           dprocess=euler.density,
+           measurement.model=measles~binom(size=cases,prob=exp(rho)),
            initializer=function(params,t0,...){
              p <- exp(params)
              with(
@@ -129,13 +131,6 @@
            }
            )
 
-set.seed(3049953)
-## simulate from the model
-tic <- Sys.time()
-x <- simulate(po,params=log(params),nsim=3)
-toc <- Sys.time()
-print(toc-tic)
-
 # alternatively, one can define the computationally intensive bits using native routines:
 ## the C codes "sir_euler_simulator" and "sir_euler_density" are included in the "examples" directory (file "sir.c")
 po <- pomp(
@@ -149,13 +144,13 @@
            paramnames=c("gamma","mu","iota","beta1","beta.sd","pop"),
            covarnames=c("seas1"),
            zeronames=c("cases"),
-           measurement.model=measles~binom(size=cases,prob=exp(rho)),
-           rprocess=euler.simulate,
            step.fun="sir_euler_simulator",
-           dprocess=euler.density,
+           rprocess=euler.simulate,
            dens.fun="sir_euler_density",
+           dprocess=euler.density,
            skeleton="sir_ODE",
            PACKAGE="pomp",
+           measurement.model=measles~binom(size=cases,prob=exp(rho)),
            initializer=function(params,t0,...){
              p <- exp(params)
              with(