allow saving of weighted particles and their weights in 'pfilter' computation

Author: kingaa

DESCRIPTION

@@ -1,8 +1,8 @@
 Package: pomp
 Type: Package
 Title: Statistical Inference for Partially Observed Markov Processes
-Version: 4.4.0.1
-Date: 2022-11-29
+Version: 4.4.1.0
+Date: 2022-11-30
 Authors@R: c(person(given=c("Aaron","A."),family="King",
 		role=c("aut","cre"),email="kingaa@umich.edu"),
 	  person(given=c("Edward","L."),family="Ionides",role=c("aut")),

R/pfilter.R

@@ -32,8 +32,12 @@
 ##' @param filter.traj logical; if \code{TRUE}, a filtered trajectory is returned for the state variables and parameters.
 ##' See \code{\link{filter.traj}} for more information.
 ##'
-##' @param save.states logical.
-##' If \code{save.states=TRUE}, the state-vector for each particle at each time is saved.
+##' @param save.states character;
+##' If \code{save.states="unweighted"}, the state-vector for each unweighted particle at each time is saved.
+##' If \code{save.states="weighted"}, the state-vector for each weighted particle at each time is saved, along with the corresponding weight.
+##' If \code{save.states="no"}, information on the latent states is not saved.
+##' \code{"FALSE"} is a synonym for \code{"no"} and \code{"TRUE"} is a synonym for \code{"unweighted"}.
+##' To retrieve the saved states, applying \code{\link{saved.states}} to the result of the \code{pfilter} computation.
 ##'
 ##' @return
 ##' An object of class \sQuote{pfilterd_pomp}, which extends class \sQuote{pomp}.
@@ -51,8 +55,8 @@
 ##'   retrieve one particle trajectory.
 ##'   Useful for building up the smoothing distribution.
 ##' }
-##' \item{\code{\link{saved.states}}}{retrieve list of saved states.}
-##' \item{\code{\link{as.data.frame}}}{ coerce to a data frame }
+##' \item{\code{\link{saved.states}}}{retrieve saved states}
+##' \item{\code{\link{as.data.frame}}}{coerce to a data frame}
 ##' \item{\code{\link{plot}}}{diagnostic plots}
 ##' }
 ##'
@@ -132,7 +136,7 @@ setMethod(
     pred.var = FALSE,
     filter.mean = FALSE,
     filter.traj = FALSE,
-    save.states = FALSE,
+    save.states = c("no", "weighted", "unweighted", "FALSE", "TRUE"),
     ...,
     verbose = getOption("verbose", FALSE)) {
 
@@ -170,7 +174,7 @@ setMethod(
     pred.var = FALSE,
     filter.mean = FALSE,
     filter.traj = FALSE,
-    save.states = FALSE,
+    save.states = c("no", "weighted", "unweighted", "FALSE", "TRUE"),
     ...,
     verbose = getOption("verbose", FALSE)) {
 
@@ -209,8 +213,9 @@ setMethod(
 
 pfilter.internal <- function (object, Np,
   pred.mean = FALSE, pred.var = FALSE, filter.mean = FALSE,
-  filter.traj = FALSE, cooling, cooling.m, save.states = FALSE, ...,
-  .gnsi = TRUE, verbose = FALSE) {
+  filter.traj = FALSE, cooling, cooling.m,
+  save.states = c("no", "weighted", "unweighted", "FALSE", "TRUE"),
+  ..., .gnsi = TRUE, verbose = FALSE) {
 
   verbose <- as.logical(verbose)
 
@@ -224,7 +229,8 @@ pfilter.internal <- function (object, Np,
   pred.var <- as.logical(pred.var)
   filter.mean <- as.logical(filter.mean)
   filter.traj <- as.logical(filter.traj)
-  save.states <- as.logical(save.states)
+  save.states <- as.character(save.states)
+  save.states <- match.arg(save.states)
 
   params <- coef(object)
   times <- time(object,t0=TRUE)
@@ -241,8 +247,11 @@ pfilter.internal <- function (object, Np,
   x <- init.x
 
   ## set up storage for saving samples from filtering distributions
-  if (save.states || filter.traj) {
+  stsav <- save.states %in% c("unweighted","TRUE")
+  wtsav <- save.states == "weighted"
+  if (stsav || wtsav || filter.traj) {
     xparticles <- setNames(vector(mode="list",length=ntimes),time(object))
+    if (wtsav) xweights <- xparticles
   }
   if (filter.traj) {
     pedigree <- vector(mode="list",length=ntimes+1)
@@ -298,6 +307,12 @@ pfilter.internal <- function (object, Np,
       times=times[nt+1],params=params,log=TRUE,.gnsi=gnsi)
     gnsi <- FALSE
 
+    ## store unweighted particles and their weights
+    if (wtsav) {
+      xparticles[[nt]] <- x
+      xweights[[nt]] <- weights
+    }
+
     ## compute prediction mean, prediction variance, filtering mean,
     ## effective sample size, log-likelihood.
     ## also do resampling.
@@ -328,7 +343,7 @@ pfilter.internal <- function (object, Np,
     if (filter.mean) filt.m[,nt] <- xx$fm
     if (filter.traj) pedigree[[nt]] <- xx$ancestry
 
-    if (save.states || filter.traj) {
+    if (stsav || filter.traj) {
       xparticles[[nt]] <- x
       dimnames(xparticles[[nt]]) <- setNames(dimnames(xparticles[[nt]]),
         c("variable",".id"))
@@ -354,7 +369,13 @@ pfilter.internal <- function (object, Np,
     }
   }
 
-  if (!save.states) xparticles <- list()
+  if (stsav) {
+    stsav <- xparticles
+  } else if (wtsav) {
+    stsav <- list(states=xparticles,weights=xweights)
+  } else {
+    stsav <- list()
+  }
 
   new(
     "pfilterd_pomp",
@@ -366,7 +387,7 @@ pfilter.internal <- function (object, Np,
     paramMatrix=array(data=numeric(0),dim=c(0,0)),
     eff.sample.size=eff.sample.size,
     cond.logLik=loglik,
-    saved.states=xparticles,
+    saved.states=stsav,
     Np=as.integer(Np),
     loglik=sum(loglik)
   )

R/saved_states.R

@@ -13,10 +13,20 @@
 ##' @family particle filter methods
 ##' @family extraction methods
 ##' @inheritParams filter.mean
+##' @param format character;
+##' format of the returned object (see below).
 ##'
-##' @return The saved states are returned in the form of a list, with one element per time-point.
-##' Each element consists of a matrix, with one row for each state variable and one column for each particle.
+##' @return According to the \code{format} argument, the saved states are returned either as a list or a data frame.
 ##' 
+##' If \code{format="data.frame"}, then the returned data frame holds the state variables and (optionally) the log weight of each particle at each observation time.
+##' The \code{.id} variable distinguishes particles.
+##'
+##' If \code{format="list"} and \code{\link{pfilter}} was called with \code{save.states="unweighted"} or \code{save.states="TRUE"}, the returned list contains one element per observation time.
+##' Each element consists of a matrix, with one row for each state variable and one column for each particle.
+##' If \code{\link{pfilter}} was called with \code{save.states="weighted"}, the list itself contains two lists:
+##' the first holds the particles as above, the second holds the corresponding log weights.
+##' In particular, it has one element per observation time; each element is the vector of per-particle log weights.
+##'
 NULL
 
 setGeneric(
@@ -45,8 +55,31 @@ setMethod(
 setMethod(
   "saved.states",
   signature=signature(object="pfilterd_pomp"),
-  definition=function (object, ...) {
-    object@saved.states
+  definition=function (object, ...,
+    format = c("list","data.frame")) {
+    format <- match.arg(format)
+    if (format=="list") {
+      object@saved.states
+    } else if (length(object@saved.states)==2L) {
+      s <- melt(object@saved.states$states)
+      w <- melt(object@saved.states$weights)
+      s$time <- time(object)[as.integer(s$L1)]
+      w$time <- time(object)[as.integer(w$L1)]
+      w$variable <- "weight"
+      x <- rbind(
+        s[,c("time",".id","variable","value")],
+        w[,c("time",".id","variable","value")]
+      )
+      x <- x[order(x$time,x$.id),]
+      rownames(x) <- NULL
+      x
+    } else {
+      s <- melt(object@saved.states)
+      s$time <- time(object)[as.integer(s$L1)]
+      s <- s[,c("time",".id","variable","value")]
+      rownames(s) <- NULL
+      s
+    }
   }
 )
 

inst/NEWS

@@ -1,5 +1,14 @@
 _N_e_w_s _f_o_r _p_a_c_k_a_g_e '_p_o_m_p'
 
+_C_h_a_n_g_e_s _i_n '_p_o_m_p' _v_e_r_s_i_o_n _4._4._1:
+
+        • It is now possible to retrieve the weighted particles
+          computed in the course of a ‘pfilter’ computation.  To
+          accomplish this, set ‘save.states="weighted"’ in the call to
+          ‘pfilter’ and retrieve the particles and their weights using
+          ‘saved.states’.  Previously, one could obtain only the
+          unweighted particles.
+
 _C_h_a_n_g_e_s _i_n '_p_o_m_p' _v_e_r_s_i_o_n _4._4._0:
 
         • Some documentation improvements.

inst/NEWS.Rd

@@ -1,5 +1,12 @@
 \name{NEWS}
 \title{News for package `pomp'}
+\section{Changes in \pkg{pomp} version 4.4.1}{
+  \itemize{
+    \item It is now possible to retrieve the weighted particles computed in the course of a \code{pfilter} computation.
+    To accomplish this, set \code{save.states="weighted"} in the call to \code{pfilter} and retrieve the particles and their weights using \code{saved.states}.
+    Previously, one could obtain only the unweighted particles.
+  }
+}
 \section{Changes in \pkg{pomp} version 4.4.0}{
   \itemize{
     \item Some documentation improvements.

man/pfilter.Rd

@@ -67,12 +67,15 @@ pfilter(pf,dmeasure=function(log,...) -Inf)
 pfilter(pf,dmeasure=function(log,...) -Inf,filter.mean=TRUE)
 
 pf1 <- pfilter(pf,save.states=TRUE,filter.traj=TRUE)
-pf2 <- pfilter(pf,pred.mean=TRUE,pred.var=TRUE,filter.mean=TRUE,save.states=TRUE)
+pf2 <- pfilter(pf,pred.mean=TRUE,pred.var=TRUE,filter.mean=TRUE,save.states="unweighted")
 pf3 <- pfilter(pf,t0=1,filter.traj=TRUE)
 pf4 <- pfilter(pf,dmeasure=Csnippet("lik = (give_log) ? R_NegInf : 0;"),
   filter.traj=TRUE)
+pf5 <- pfilter(pf,save.states="weighted")
 pf1 %>% saved.states() %>% melt() %>% names()
 pf1 %>% saved.states() %>% melt() %>% dim()
+pf1 %>% saved.states(format="data") %>% names()
+pf1 %>% saved.states(format="data") %>% dim()
 c(A=pf1,B=pf2) %>% saved.states() %>% melt() %>% names()
 c(A=pf1,B=pf2) %>% saved.states() %>% melt() %>% sapply(class)
 c(A=pf1,B=pf2) %>% as.data.frame() %>% sapply(class)
@@ -87,7 +90,8 @@ names(dimnames(filter.traj(c(pf1,pf4))))
 names(melt(as(c(pf1,pf4),"data.frame")))
 pf2 %>% melt() %>% names()
 pf2 %>% melt(id="time") %>% names()
-
+pf5 %>% saved.states(format="d") %>% names()
+pf5 %>% saved.states(format="d") %>% dim()
 try(saved.states())
 try(saved.states(NULL))
 try(saved.states("bob"))

man/saved_states.Rd

@@ -148,14 +148,19 @@ Error : in 'pfilter': in 'dmeasure': ouch!
 <object of class 'pfilterd_pomp'>
 > 
 > pf1 <- pfilter(pf,save.states=TRUE,filter.traj=TRUE)
-> pf2 <- pfilter(pf,pred.mean=TRUE,pred.var=TRUE,filter.mean=TRUE,save.states=TRUE)
+> pf2 <- pfilter(pf,pred.mean=TRUE,pred.var=TRUE,filter.mean=TRUE,save.states="unweighted")
 > pf3 <- pfilter(pf,t0=1,filter.traj=TRUE)
 > pf4 <- pfilter(pf,dmeasure=Csnippet("lik = (give_log) ? R_NegInf : 0;"),
 +   filter.traj=TRUE)
+> pf5 <- pfilter(pf,save.states="weighted")
 > pf1 %>% saved.states() %>% melt() %>% names()
 [1] "variable" ".id"      "value"    "L1"      
 > pf1 %>% saved.states() %>% melt() %>% dim()
 [1] 20000     4
+> pf1 %>% saved.states(format="data") %>% names()
+[1] "time"     ".id"      "variable" "value"   
+> pf1 %>% saved.states(format="data") %>% dim()
+[1] 20000     4
 > c(A=pf1,B=pf2) %>% saved.states() %>% melt() %>% names()
 [1] "variable" ".id"      "value"    "L2"       "L1"      
 > c(A=pf1,B=pf2) %>% saved.states() %>% melt() %>% sapply(class)
@@ -211,7 +216,10 @@ No id variables; using all as measure variables
 [1] "variable" "value"   
 > pf2 %>% melt(id="time") %>% names()
 [1] "time"     "variable" "value"   
-> 
+> pf5 %>% saved.states(format="d") %>% names()
+[1] "time"     ".id"      "variable" "value"   
+> pf5 %>% saved.states(format="d") %>% dim()
+[1] 30000     4
 > try(saved.states())
 Error : in 'saved.states': 'object' is a required argument.
 > try(saved.states(NULL))