[Pomp-commits] r46 - pkg/tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Aug 15 12:33:57 CEST 2008
Author: kingaa
Date: 2008-08-15 12:33:57 +0200 (Fri, 15 Aug 2008)
New Revision: 46
Modified:
pkg/tests/sir.R
pkg/tests/sir.Rout.save
Log:
improve test of SIR model
Modified: pkg/tests/sir.R
===================================================================
--- pkg/tests/sir.R 2008-08-15 10:32:57 UTC (rev 45)
+++ pkg/tests/sir.R 2008-08-15 10:33:57 UTC (rev 46)
@@ -129,59 +129,61 @@
x <- simulate(po,params=log(params),nsim=3)
toc <- Sys.time()
print(toc-tic)
+
pdf(file='sir.pdf')
+
plot(x[[1]],variables=c("S","I","R","cases","W"))
-t <- seq(0,4/52,by=1/52/25)
-X <- simulate(po,params=log(params),nsim=10,states=TRUE,obs=TRUE,times=t)
+t1 <- seq(0,4/52,by=1/52/25)
+X1 <- simulate(po,params=log(params),nsim=10,states=TRUE,obs=TRUE,times=t1)
-f <- dprocess(
- po,
- x=X$states[,,31:40],
- times=t[31:40],
- params=matrix(
- log(params),
- nrow=length(params),
- ncol=10,
- dimnames=list(names(params),NULL)
- ),
- log=TRUE
- )
-print(apply(f,1,sum),digits=4)
+t2 <- seq(0,2,by=1/52)
+X2 <- simulate(po,params=log(params),nsim=1,states=TRUE,obs=TRUE,times=t2)
-g <- dmeasure(
- po,
- y=rbind(measles=X$obs[,7,]),
- x=X$states,
- times=t,
- params=matrix(
- log(params),
- nrow=length(params),
- ncol=10,
- dimnames=list(names(params),NULL)
- ),
- log=TRUE
- )
-print(apply(g,1,sum),digits=4)
-
-t <- seq(0,2,by=1/52)
-X <- simulate(po,params=log(params),nsim=1,states=TRUE,obs=TRUE,times=t)
-
-h <- skeleton(
- po,
- x=X$states[,1,55:70,drop=FALSE],
- t=t[55:70],
- params=as.matrix(log(params))
- )
-print(h[c("S","I","R"),,],digits=4)
-
-t <- seq(0,20,by=1/52)
+t3 <- seq(0,20,by=1/52)
tic <- Sys.time()
-X <- trajectory(po,params=log(params),times=t,hmax=1/52)
+X3 <- trajectory(po,params=log(params),times=t3,hmax=1/52)
toc <- Sys.time()
print(toc-tic)
-plot(t,X['I',1,],type='l')
+plot(t3,X3['I',1,],type='l')
+f1 <- dprocess(
+ po,
+ x=X1$states[,,31:40],
+ times=t1[31:40],
+ params=matrix(
+ log(params),
+ nrow=length(params),
+ ncol=10,
+ dimnames=list(names(params),NULL)
+ ),
+ log=TRUE
+ )
+print(apply(f1,1,sum),digits=4)
+
+g1 <- dmeasure(
+ po,
+ y=rbind(measles=X1$obs[,7,]),
+ x=X1$states,
+ times=t1,
+ params=matrix(
+ log(params),
+ nrow=length(params),
+ ncol=10,
+ dimnames=list(names(params),NULL)
+ ),
+ log=TRUE
+ )
+print(apply(g1,1,sum),digits=4)
+
+h1 <- skeleton(
+ po,
+ x=X2$states[,1,55:70,drop=FALSE],
+ t=t2[55:70],
+ params=as.matrix(log(params))
+ )
+print(h1[c("S","I","R"),,],digits=4)
+
po <- pomp(
times=seq(1/52,4,by=1/52),
data=rbind(measles=numeric(52*4)),
@@ -225,53 +227,52 @@
print(toc-tic)
plot(x[[1]],variables=c("S","I","R","cases","W"))
-t <- seq(0,4/52,by=1/52/25)
-X <- simulate(po,params=log(params),nsim=10,states=TRUE,obs=TRUE,times=t)
+t3 <- seq(0,20,by=1/52)
+tic <- Sys.time()
+X4 <- trajectory(po,params=log(params),times=t3,hmax=1/52)
+toc <- Sys.time()
+print(toc-tic)
+plot(t3,X4['I',1,],type='l')
-f <- dprocess(
- po,
- x=X$states[,,31:40],
- times=t[31:40],
- params=matrix(
- log(params),
- nrow=length(params),
- ncol=10,
- dimnames=list(names(params),NULL)
- ),
- log=TRUE
- )
-print(apply(f,1,sum),digits=4)
+f2 <- dprocess(
+ po,
+ x=X1$states[,,31:40],
+ times=t1[31:40],
+ params=matrix(
+ log(params),
+ nrow=length(params),
+ ncol=10,
+ dimnames=list(names(params),NULL)
+ ),
+ log=TRUE
+ )
+print(apply(f2,1,sum),digits=4)
-g <- dmeasure(
- po,
- y=rbind(measles=X$obs[,7,]),
- x=X$states,
- times=t,
- params=matrix(
- log(params),
- nrow=length(params),
- ncol=10,
- dimnames=list(names(params),NULL)
- ),
- log=TRUE
- )
-print(apply(g,1,sum),digits=4)
+g2 <- dmeasure(
+ po,
+ y=rbind(measles=X1$obs[,7,]),
+ x=X1$states,
+ times=t1,
+ params=matrix(
+ log(params),
+ nrow=length(params),
+ ncol=10,
+ dimnames=list(names(params),NULL)
+ ),
+ log=TRUE
+ )
+print(apply(g2,1,sum),digits=4)
-t <- seq(0,2,by=1/52)
-X <- simulate(po,params=log(params),nsim=1,states=TRUE,obs=TRUE,times=t)
+h2 <- skeleton(
+ po,
+ x=X2$states[,1,55:70,drop=FALSE],
+ t=t2[55:70],
+ params=as.matrix(log(params))
+ )
+print(h2[c("S","I","R"),,],digits=4)
-h <- skeleton(
- po,
- x=X$states[,1,55:70,drop=FALSE],
- t=t[55:70],
- params=as.matrix(log(params))
- )
-print(h[c("S","I","R"),,],digits=4)
+print(max(abs(f2-f1),na.rm=T),digits=4)
+print(max(abs(g2-g1),na.rm=T),digits=4)
+print(max(abs(h2-h1),na.rm=T),digits=4)
-t <- seq(0,20,by=1/52)
-tic <- Sys.time()
-X <- trajectory(po,params=log(params),times=t,hmax=1/52)
-toc <- Sys.time()
-print(toc-tic)
-plot(t,X['I',1,],type='l')
dev.off()
Modified: pkg/tests/sir.Rout.save
===================================================================
--- pkg/tests/sir.Rout.save 2008-08-15 10:32:57 UTC (rev 45)
+++ pkg/tests/sir.Rout.save 2008-08-15 10:33:57 UTC (rev 46)
@@ -147,71 +147,73 @@
> x <- simulate(po,params=log(params),nsim=3)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 4.001635 secs
+Time difference of 3.935402 secs
+>
> pdf(file='sir.pdf')
+>
> plot(x[[1]],variables=c("S","I","R","cases","W"))
>
-> t <- seq(0,4/52,by=1/52/25)
-> X <- simulate(po,params=log(params),nsim=10,states=TRUE,obs=TRUE,times=t)
+> t1 <- seq(0,4/52,by=1/52/25)
+> X1 <- simulate(po,params=log(params),nsim=10,states=TRUE,obs=TRUE,times=t1)
>
-> f <- dprocess(
-+ po,
-+ x=X$states[,,31:40],
-+ times=t[31:40],
-+ params=matrix(
-+ log(params),
-+ nrow=length(params),
-+ ncol=10,
-+ dimnames=list(names(params),NULL)
-+ ),
-+ log=TRUE
-+ )
-> print(apply(f,1,sum),digits=4)
- [1] -56.62 -41.79 -51.58 -47.57 -47.92 -44.74 -42.84 -46.91 -55.07 -47.48
+> t2 <- seq(0,2,by=1/52)
+> X2 <- simulate(po,params=log(params),nsim=1,states=TRUE,obs=TRUE,times=t2)
>
-> g <- dmeasure(
-+ po,
-+ y=rbind(measles=X$obs[,7,]),
-+ x=X$states,
-+ times=t,
-+ params=matrix(
-+ log(params),
-+ nrow=length(params),
-+ ncol=10,
-+ dimnames=list(names(params),NULL)
-+ ),
-+ log=TRUE
-+ )
-> print(apply(g,1,sum),digits=4)
- [1] -Inf -Inf -Inf -Inf -Inf -Inf -260.9 -Inf -Inf -Inf
->
-> t <- seq(0,2,by=1/52)
-> X <- simulate(po,params=log(params),nsim=1,states=TRUE,obs=TRUE,times=t)
->
-> h <- skeleton(
-+ po,
-+ x=X$states[,1,55:70,drop=FALSE],
-+ t=t[55:70],
-+ params=as.matrix(log(params))
-+ )
-> print(h[c("S","I","R"),,],digits=4)
- [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
-S 39923.7 39832.3 39214.4 39275.9 38864 38820.4 37892 36798 36093
-I 184.5 234.1 420.2 445.7 608 662.3 1027 1481 1806
-R -40115.5 -40072.7 -39640.2 -39728.2 -39479 -39488.1 -38927 -38287 -37908
- [,10] [,11] [,12] [,13] [,14] [,15] [,16]
-S 33608 28740 22452 16799 5834 -8382 -27621.9
-I 2847 4893 7554 9945 14452 20004 27020.1
-R -36462 -33640 -30013 -26751 -20293 -11629 593.4
->
-> t <- seq(0,20,by=1/52)
+> t3 <- seq(0,20,by=1/52)
> tic <- Sys.time()
-> X <- trajectory(po,params=log(params),times=t,hmax=1/52)
+> X3 <- trajectory(po,params=log(params),times=t3,hmax=1/52)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 9.436516 secs
-> plot(t,X['I',1,],type='l')
+Time difference of 7.444703 secs
+> plot(t3,X3['I',1,],type='l')
>
+> f1 <- dprocess(
++ po,
++ x=X1$states[,,31:40],
++ times=t1[31:40],
++ params=matrix(
++ log(params),
++ nrow=length(params),
++ ncol=10,
++ dimnames=list(names(params),NULL)
++ ),
++ log=TRUE
++ )
+> print(apply(f1,1,sum),digits=4)
+ [1] -51.57 -43.19 -50.14 -39.65 -41.63 -49.34 -39.60 -48.19 -42.60 -48.26
+>
+> g1 <- dmeasure(
++ po,
++ y=rbind(measles=X1$obs[,7,]),
++ x=X1$states,
++ times=t1,
++ params=matrix(
++ log(params),
++ nrow=length(params),
++ ncol=10,
++ dimnames=list(names(params),NULL)
++ ),
++ log=TRUE
++ )
+> print(apply(g1,1,sum),digits=4)
+ [1] -Inf -Inf -Inf -Inf -Inf -Inf -263.6 -408.1 -Inf -468.1
+>
+> h1 <- skeleton(
++ po,
++ x=X2$states[,1,55:70,drop=FALSE],
++ t=t2[55:70],
++ params=as.matrix(log(params))
++ )
+> print(h1[c("S","I","R"),,],digits=4)
+ [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
+S 39362.0 39111.6 39195.1 38583.3 38374.5 36367 34769 33057 29765
+I 272.3 376.1 407.8 615.6 731.6 1423 2046 2755 4083
+R -39640.5 -39495.0 -39609.9 -39205.0 -39112.3 -37796 -36819 -35816 -33852
+ [,10] [,11] [,12] [,13] [,14] [,15] [,16]
+S 24577 17622 11393 319.5 -14707 -33518 -57305
+I 6165 8975 11503 15834.2 21418 27852 34849
+R -30746 -26601 -22900 -16156.5 -6714 5662 22453
+>
> po <- pomp(
+ times=seq(1/52,4,by=1/52),
+ data=rbind(measles=numeric(52*4)),
@@ -253,69 +255,71 @@
> x <- simulate(po,params=log(params),nsim=3)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 0.3902411 secs
+Time difference of 0.3896399 secs
> plot(x[[1]],variables=c("S","I","R","cases","W"))
>
-> t <- seq(0,4/52,by=1/52/25)
-> X <- simulate(po,params=log(params),nsim=10,states=TRUE,obs=TRUE,times=t)
+> t3 <- seq(0,20,by=1/52)
+> tic <- Sys.time()
+> X4 <- trajectory(po,params=log(params),times=t3,hmax=1/52)
+> toc <- Sys.time()
+> print(toc-tic)
+Time difference of 5.439281 secs
+> plot(t3,X4['I',1,],type='l')
>
-> f <- dprocess(
-+ po,
-+ x=X$states[,,31:40],
-+ times=t[31:40],
-+ params=matrix(
-+ log(params),
-+ nrow=length(params),
-+ ncol=10,
-+ dimnames=list(names(params),NULL)
-+ ),
-+ log=TRUE
-+ )
-> print(apply(f,1,sum),digits=4)
- [1] -54.89 -50.09 -47.17 -35.39 -44.83 -44.58 -47.16 -35.21 -46.25 -46.85
+> f2 <- dprocess(
++ po,
++ x=X1$states[,,31:40],
++ times=t1[31:40],
++ params=matrix(
++ log(params),
++ nrow=length(params),
++ ncol=10,
++ dimnames=list(names(params),NULL)
++ ),
++ log=TRUE
++ )
+> print(apply(f2,1,sum),digits=4)
+ [1] -51.57 -43.19 -50.14 -39.65 -41.63 -49.34 -39.60 -48.19 -42.60 -48.26
>
-> g <- dmeasure(
-+ po,
-+ y=rbind(measles=X$obs[,7,]),
-+ x=X$states,
-+ times=t,
-+ params=matrix(
-+ log(params),
-+ nrow=length(params),
-+ ncol=10,
-+ dimnames=list(names(params),NULL)
-+ ),
-+ log=TRUE
-+ )
-> print(apply(g,1,sum),digits=4)
- [1] -Inf -Inf -Inf -Inf -Inf -Inf -261.2 -Inf -Inf -Inf
+> g2 <- dmeasure(
++ po,
++ y=rbind(measles=X1$obs[,7,]),
++ x=X1$states,
++ times=t1,
++ params=matrix(
++ log(params),
++ nrow=length(params),
++ ncol=10,
++ dimnames=list(names(params),NULL)
++ ),
++ log=TRUE
++ )
+> print(apply(g2,1,sum),digits=4)
+ [1] -Inf -Inf -Inf -Inf -Inf -Inf -263.6 -408.1 -Inf -468.1
>
-> t <- seq(0,2,by=1/52)
-> X <- simulate(po,params=log(params),nsim=1,states=TRUE,obs=TRUE,times=t)
+> h2 <- skeleton(
++ po,
++ x=X2$states[,1,55:70,drop=FALSE],
++ t=t2[55:70],
++ params=as.matrix(log(params))
++ )
+> print(h2[c("S","I","R"),,],digits=4)
+ [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
+S 39362.0 39111.6 39195.1 38583.3 38374.5 36367 34769 33057 29765
+I 272.3 376.1 407.8 615.6 731.6 1423 2046 2755 4083
+R -39640.5 -39495.0 -39609.9 -39205.0 -39112.3 -37796 -36819 -35816 -33852
+ [,10] [,11] [,12] [,13] [,14] [,15] [,16]
+S 24577 17622 11393 319.5 -14707 -33518 -57305
+I 6165 8975 11503 15834.2 21418 27852 34849
+R -30746 -26601 -22900 -16156.5 -6714 5662 22453
>
-> h <- skeleton(
-+ po,
-+ x=X$states[,1,55:70,drop=FALSE],
-+ t=t[55:70],
-+ params=as.matrix(log(params))
-+ )
-> print(h[c("S","I","R"),,],digits=4)
- [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
-S 37319.0 37046.3 35259 33635 31996 30625 27251 25043 16279 7989
-I 644.9 823.6 1374 1969 2632 3261 4593 5593 8989 12303
-R -37964.4 -37871.3 -36633 -35605 -34628 -33886 -31845 -30638 -25270 -20294
- [,11] [,12] [,13] [,14] [,15] [,16]
-S -4107 -21296 -45772 -70414 -100996 -132252
-I 17033 23456 31943 39251 46266 49391
-R -12929 -2160 13830 31164 54730 82862
+> print(max(abs(f2-f1),na.rm=T),digits=4)
+[1] 7.994e-15
+> print(max(abs(g2-g1),na.rm=T),digits=4)
+[1] 0
+> print(max(abs(h2-h1),na.rm=T),digits=4)
+[1] 8.731e-11
>
-> t <- seq(0,20,by=1/52)
-> tic <- Sys.time()
-> X <- trajectory(po,params=log(params),times=t,hmax=1/52)
-> toc <- Sys.time()
-> print(toc-tic)
-Time difference of 7.409947 secs
-> plot(t,X['I',1,],type='l')
> dev.off()
null device
1
More information about the pomp-commits
mailing list