[Pomp-commits] r181 - in pkg: R man tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Dec 9 18:20:40 CET 2009
Author: kingaa
Date: 2009-12-09 18:20:39 +0100 (Wed, 09 Dec 2009)
New Revision: 181
Added:
pkg/tests/ou2-forecast.R
pkg/tests/ou2-forecast.Rout.save
Modified:
pkg/R/mif.R
pkg/R/pfilter-mif.R
pkg/R/pfilter.R
pkg/man/pfilter.Rd
Log:
- in pfilter, add ability to save individual particles' state vectors. this feature is useful for forecasting.
Modified: pkg/R/mif.R
===================================================================
--- pkg/R/mif.R 2009-12-09 13:04:54 UTC (rev 180)
+++ pkg/R/mif.R 2009-12-09 17:20:39 UTC (rev 181)
@@ -239,6 +239,7 @@
pred.mean=(n==Nmif),
pred.var=(weighted||(n==Nmif)),
filter.mean=TRUE,
+ save.states=FALSE,
.rw.sd=sigma.n[pars],
verbose=verbose
),
Modified: pkg/R/pfilter-mif.R
===================================================================
--- pkg/R/pfilter-mif.R 2009-12-09 13:04:54 UTC (rev 180)
+++ pkg/R/pfilter-mif.R 2009-12-09 17:20:39 UTC (rev 181)
@@ -5,11 +5,18 @@
tol = 1e-17, warn = TRUE, max.fail = 0,
pred.mean = FALSE,
pred.var = FALSE,
- filter.mean = FALSE, ...) {
+ filter.mean = FALSE,
+ ...) {
if (missing(params))
params <- coef(object)
if (missing(Np))
Np <- object at alg.pars$Np
+ if ("save.states"%in%names(list(...)))
+ stop(
+ "pfilter error: you cannot set ",sQuote("save.states"),
+ " when ",sQuote("object")," is of class mif",
+ call.=FALSE
+ )
pfilter(
object=as(object,"pomp"),
params=params,
@@ -20,6 +27,7 @@
pred.mean=pred.mean,
pred.var=pred.var,
filter.mean=filter.mean,
+ save.states=FALSE,
...)
}
)
Modified: pkg/R/pfilter.R
===================================================================
--- pkg/R/pfilter.R 2009-12-09 13:04:54 UTC (rev 180)
+++ pkg/R/pfilter.R 2009-12-09 17:20:39 UTC (rev 181)
@@ -6,7 +6,8 @@
pfilter.internal <- function (object, params, Np,
tol, warn, max.fail,
pred.mean, pred.var, filter.mean,
- .rw.sd, verbose) {
+ .rw.sd, verbose,
+ save.states) {
if (missing(params)) {
params <- coef(object)
if (length(params)==0) {
@@ -35,6 +36,15 @@
x <- init.state(object,params=params)
statenames <- rownames(x)
nvars <- nrow(x)
+ if (save.states) {
+ xparticles <- array(
+ data=NA,
+ dim=c(nvars,Np,ntimes),
+ dimnames=list(statenames,NULL,NULL)
+ )
+ } else {
+ xparticles <- NULL
+ }
random.walk <- !missing(.rw.sd)
if (random.walk) {
@@ -215,6 +225,10 @@
params[rw.names,] <- params[rw.names,]+rnorm(n=Np*length(sigma),mean=0,sd=sigma)
}
+ if (save.states) {
+ xparticles[,,nt] <- x
+ }
+
if (verbose && ((ntimes-nt)%%5==0))
cat("pfilter timestep",nt,"of",ntimes,"finished\n")
@@ -226,6 +240,7 @@
filter.mean=filt.m,
eff.sample.size=eff.sample.size,
cond.loglik=loglik,
+ states=xparticles,
nfail=nfail,
loglik=sum(loglik)
)
@@ -239,7 +254,9 @@
pred.mean = FALSE,
pred.var = FALSE,
filter.mean = FALSE,
- verbose = FALSE, ...) {
+ save.states = FALSE,
+ verbose = FALSE,
+ ...) {
pfilter.internal(
object=object,
params=params,
@@ -250,6 +267,7 @@
pred.mean=pred.mean,
pred.var=pred.var,
filter.mean=filter.mean,
+ save.states=save.states,
verbose=verbose
)
}
Modified: pkg/man/pfilter.Rd
===================================================================
--- pkg/man/pfilter.Rd 2009-12-09 13:04:54 UTC (rev 180)
+++ pkg/man/pfilter.Rd 2009-12-09 17:20:39 UTC (rev 181)
@@ -13,10 +13,9 @@
pfilter(object, \dots)
\S4method{pfilter}{pomp}(object, params, Np, tol = 1e-17,
warn = TRUE, max.fail = 0, pred.mean = FALSE, pred.var = FALSE,
- filter.mean = FALSE, verbose = FALSE, \dots)
+ filter.mean = FALSE, save.states = FALSE, verbose = FALSE, \dots)
\S4method{pfilter}{mif}(object, params, Np, tol = 1e-17, warn = TRUE,
- max.fail = 0, pred.mean = FALSE, pred.var = FALSE,
- filter.mean = FALSE, \dots)
+ pred.mean = FALSE, pred.var = FALSE, filter.mean = FALSE, \dots)
}
\arguments{
\item{object}{
@@ -53,6 +52,9 @@
\item{filter.mean}{
If \code{TRUE}, the filtering means are calculated for the state variables and parameters.
}
+ \item{save.states}{
+ If \code{TRUE}, the state-vector for each particle is saved and returned.
+ }
\item{verbose}{
If \code{TRUE}, progress information is reported as \code{pfilter} works.
}
@@ -78,6 +80,11 @@
\item{cond.loglik}{
A vector containing the conditional log likelihoods at each time point.
}
+ \item{states}{
+ If \code{saves.states=TRUE}, the array of state-vectors at each time point, for each particle.
+ An array with dimensions \code{nvars}-by-\code{Np}-by-\code{ntimes}.
+ In particular, \code{states[,i,t]} can be considered a sample from \eqn{f[X|y_{1:t}]}.
+ }
\item{nfail}{
The number of filtering failures encountered.
}
Added: pkg/tests/ou2-forecast.R
===================================================================
--- pkg/tests/ou2-forecast.R (rev 0)
+++ pkg/tests/ou2-forecast.R 2009-12-09 17:20:39 UTC (rev 181)
@@ -0,0 +1,20 @@
+library(pomp)
+
+set.seed(921625222L)
+
+pf <- pfilter(ou2,Np=1000,save.states=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$states[c("x1","x2"),,t]
+ y <- simulate(ou2,params=pp,obs=TRUE,times=time(ou2)[t:(t+10)])
+ mn <- apply(y,c(1,3),mean)
+ sd <- apply(y,c(1,3),sd)
+ z <- (data.array(ou2)[,t:(t+10)]-mn)/sd ## z score
+ mse <- (data.array(ou2)[,t:(t+10)]-mn)^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)
+try(
+ pfilter(fit,save.states=TRUE)
+ )
+
Added: pkg/tests/ou2-forecast.Rout.save
===================================================================
--- pkg/tests/ou2-forecast.Rout.save (rev 0)
+++ pkg/tests/ou2-forecast.Rout.save 2009-12-09 17:20:39 UTC (rev 181)
@@ -0,0 +1,42 @@
+
+R version 2.9.1 (2009-06-26)
+Copyright (C) 2009 The R Foundation for Statistical Computing
+ISBN 3-900051-07-0
+
+R is free software and comes with ABSOLUTELY NO WARRANTY.
+You are welcome to redistribute it under certain conditions.
+Type 'license()' or 'licence()' for distribution details.
+
+R is a collaborative project with many contributors.
+Type 'contributors()' for more information and
+'citation()' on how to cite R or R packages in publications.
+
+Type 'demo()' for some demos, 'help()' for on-line help, or
+'help.start()' for an HTML browser interface to help.
+Type 'q()' to quit R.
+
+> library(pomp)
+Loading required package: deSolve
+Loading required package: subplex
+Loading required package: mvtnorm
+>
+> set.seed(921625222L)
+>
+> pf <- pfilter(ou2,Np=1000,save.states=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$states[c("x1","x2"),,t]
++ y <- simulate(ou2,params=pp,obs=TRUE,times=time(ou2)[t:(t+10)])
++ mn <- apply(y,c(1,3),mean)
++ sd <- apply(y,c(1,3),sd)
++ z <- (data.array(ou2)[,t:(t+10)]-mn)/sd ## z score
++ mse <- (data.array(ou2)[,t:(t+10)]-mn)^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)
+> try(
++ pfilter(fit,save.states=TRUE)
++ )
+Error : pfilter error: you cannot set 'save.states' when 'object' is of class mif
+>
+>
More information about the pomp-commits
mailing list