[Pomp-commits] r473 - pkg/tests

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

Author: kingaa
Date: 2011-05-11 23:39:37 +0200 (Wed, 11 May 2011)
New Revision: 473

Modified:
   pkg/tests/ou2-forecast.R
   pkg/tests/ou2-forecast.Rout.save
Log:
- update this test to use new 'save.params' option


Modified: pkg/tests/ou2-forecast.R
===================================================================
--- pkg/tests/ou2-forecast.R	2011-05-11 21:38:40 UTC (rev 472)
+++ pkg/tests/ou2-forecast.R	2011-05-11 21:39:37 UTC (rev 473)
@@ -3,17 +3,35 @@
 set.seed(921625222L)
 
 data(ou2)
-pf <- pfilter(ou2,Np=1000,save.states=TRUE)
+tm <- time(ou2)
+y <- obs(ou2)
+theta <- coef(ou2)
+pf <- pfilter(ou2,Np=1000,save.states=TRUE,save.params=TRUE)
 ll <- cumsum(pf$cond.loglik)
-pp <- matrix(data=coef(ou2),nrow=length(coef(ou2)),ncol=1000,dimnames=list(names(coef(ou2)),NULL))
-for (t in seq(60,90,by=10)) {
-  pp[c("x1.0","x2.0"),] <- pf$last.states[c("x1","x2"),,t]
-  y <- simulate(ou2,params=pp,obs=TRUE,t0=time(ou2)[t],times=time(ou2)[(t+1):(t+10)])
-  mn <- apply(y,c(1,3),mean)
-  sd <- apply(y,c(1,3),sd)
-  z <- (data.array(ou2)[,(t+1):(t+10)]-mn)/sd       ## z score
-  mse <- (data.array(ou2)[,(t+1):(t+10)]-mn)^2+sd^2 ## mean squared error
+pp <- matrix(
+             data=theta,
+             nrow=length(theta),
+             ncol=1000,
+             dimnames=list(names(theta),NULL)
+             )
+z <- array(dim=c(2,9,10))
+mse <- array(dim=c(2,9,10))
+t0 <- seq(from=10,to=90,by=10)
+for (k in seq_along(t0)) {
+  pp[c("x1.0","x2.0"),] <- pf$saved.states[c("x1","x2"),,tm==t0[k]]
+  inds <- which(tm>t0[k]&tm<=t0[k]+10)
+  Y <- simulate(ou2,params=pp,obs=TRUE,t0=t0[k],times=tm[inds])
+  mn <- apply(Y,c(1,3),mean)
+  sd <- apply(Y,c(1,3),sd)
+  bias <- mn-y[,inds]
+  z[,k,] <- bias/sd                               ## z score
+  mse[,k,] <- bias^2+sd^2                         ## mean squared error
 }
+
 fit <- mif(ou2,Nmif=3,rw.sd=c(alpha.1=0.1,alpha.4=0.1),Np=1000,cooling.factor=0.98,var.factor=1,ic.lag=2)
-pf <- pfilter(fit,save.states=TRUE)
+pf <- pfilter(fit,save.states=TRUE,save.params=TRUE)
 
+pdf(file="ou2-forecast.pdf")
+matplot(t(mse[1,,]),xlab="horizon",ylab="MSE",type='l',lty=1,col=hsv(h=seq(from=0,by=0.1,length=9)))
+legend("topleft",bty='n',lty=1,col=hsv(h=seq(from=0,by=0.1,length=9)),legend=t0)
+dev.off()

Modified: pkg/tests/ou2-forecast.Rout.save
===================================================================
--- pkg/tests/ou2-forecast.Rout.save	2011-05-11 21:38:40 UTC (rev 472)
+++ pkg/tests/ou2-forecast.Rout.save	2011-05-11 21:39:37 UTC (rev 473)
@@ -1,7 +1,8 @@
 
-R version 2.11.1 (2010-05-31)
-Copyright (C) 2010 The R Foundation for Statistical Computing
+R version 2.12.2 (2011-02-25)
+Copyright (C) 2011 The R Foundation for Statistical Computing
 ISBN 3-900051-07-0
+Platform: x86_64-unknown-linux-gnu (64-bit)
 
 R is free software and comes with ABSOLUTELY NO WARRANTY.
 You are welcome to redistribute it under certain conditions.
@@ -16,22 +17,45 @@
 Type 'q()' to quit R.
 
 > library(pomp)
+Loading required package: mvtnorm
+Loading required package: subplex
+Loading required package: deSolve
 > 
 > set.seed(921625222L)
 > 
 > data(ou2)
-> pf <- pfilter(ou2,Np=1000,save.states=TRUE)
+> tm <- time(ou2)
+> y <- obs(ou2)
+> theta <- coef(ou2)
+> pf <- pfilter(ou2,Np=1000,save.states=TRUE,save.params=TRUE)
 > ll <- cumsum(pf$cond.loglik)
-> pp <- matrix(data=coef(ou2),nrow=length(coef(ou2)),ncol=1000,dimnames=list(names(coef(ou2)),NULL))
-> for (t in seq(60,90,by=10)) {
-+   pp[c("x1.0","x2.0"),] <- pf$last.states[c("x1","x2"),,t]
-+   y <- simulate(ou2,params=pp,obs=TRUE,t0=time(ou2)[t],times=time(ou2)[(t+1):(t+10)])
-+   mn <- apply(y,c(1,3),mean)
-+   sd <- apply(y,c(1,3),sd)
-+   z <- (data.array(ou2)[,(t+1):(t+10)]-mn)/sd       ## z score
-+   mse <- (data.array(ou2)[,(t+1):(t+10)]-mn)^2+sd^2 ## mean squared error
+> pp <- matrix(
++              data=theta,
++              nrow=length(theta),
++              ncol=1000,
++              dimnames=list(names(theta),NULL)
++              )
+> z <- array(dim=c(2,9,10))
+> mse <- array(dim=c(2,9,10))
+> t0 <- seq(from=10,to=90,by=10)
+> for (k in seq_along(t0)) {
++   pp[c("x1.0","x2.0"),] <- pf$saved.states[c("x1","x2"),,tm==t0[k]]
++   inds <- which(tm>t0[k]&tm<=t0[k]+10)
++   Y <- simulate(ou2,params=pp,obs=TRUE,t0=t0[k],times=tm[inds])
++   mn <- apply(Y,c(1,3),mean)
++   sd <- apply(Y,c(1,3),sd)
++   bias <- mn-y[,inds]
++   z[,k,] <- bias/sd                               ## z score
++   mse[,k,] <- bias^2+sd^2                         ## mean squared error
 + }
+> 
 > fit <- mif(ou2,Nmif=3,rw.sd=c(alpha.1=0.1,alpha.4=0.1),Np=1000,cooling.factor=0.98,var.factor=1,ic.lag=2)
-> pf <- pfilter(fit,save.states=TRUE)
+> pf <- pfilter(fit,save.states=TRUE,save.params=TRUE)
 > 
+> pdf(file="ou2-forecast.pdf")
+> matplot(t(mse[1,,]),xlab="horizon",ylab="MSE",type='l',lty=1,col=hsv(h=seq(from=0,by=0.1,length=9)))
+> legend("topleft",bty='n',lty=1,col=hsv(h=seq(from=0,by=0.1,length=9)),legend=t0)
+> dev.off()
+null device 
+          1 
 >