pfilter(..., save.states) + documentation edits

R/pfilter.R

@@ -28,10 +28,10 @@
 ##' @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 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"}.
+##' If \code{save.states="filter"}, the state-vector for each filtered particle \eqn{X_{n, j}^F} at each time \eqn{n} is saved.
+##' If \code{save.states="prediction"}, the state-vector for each prediction particle \eqn{X_{n, j}^P} at each time \eqn{n} is saved, along with the corresponding weight \eqn{w_{n, j} = f_{Y_n|X_n}(y^*|X_{n, j}^P;\theta)}.
+##' The old values "unweighted", "weighted", TRUE, and FALSE are deprecated and will issue a warning if used, mapping to the new values for backward compatibility.
 ##' To retrieve the saved states, apply \code{\link{saved_states}} to the result of the \code{pfilter} computation.
 ##' @param ... additional arguments are passed to \code{\link{pomp}}.
 ##' This allows one to set, unset, or modify \link[=basic_components]{basic model components} within a call to this function.
@@ -129,7 +129,7 @@ setMethod(
     pred.var = FALSE,
     filter.mean = FALSE,
     filter.traj = FALSE,
-    save.states = c("no", "weighted", "unweighted", "FALSE", "TRUE"),
+    save.states = c("no", "filter", "prediction"),
     verbose = getOption("verbose", FALSE)
   ) {
 
@@ -168,7 +168,7 @@ setMethod(
     pred.var = FALSE,
     filter.mean = FALSE,
     filter.traj = FALSE,
-    save.states = c("no", "weighted", "unweighted", "FALSE", "TRUE"),
+    save.states = c("no", "filter", "prediction"),
     verbose = getOption("verbose", FALSE)
   ) {
 
@@ -213,7 +213,7 @@ pfilter_internal <- function (
   ...,
   pred.mean = FALSE, pred.var = FALSE, filter.mean = FALSE,
   filter.traj = FALSE, cooling, cooling.m,
-  save.states = c("no", "weighted", "unweighted", "FALSE", "TRUE"),
+  save.states = c("no", "filter", "prediction"),
   .gnsi = TRUE, verbose = FALSE
 ) {
 
@@ -229,8 +229,27 @@ pfilter_internal <- function (
   pred.var <- as.logical(pred.var)
   filter.mean <- as.logical(filter.mean)
   filter.traj <- as.logical(filter.traj)
-  save.states <- as.character(save.states)
-  save.states <- match.arg(save.states)
+
+  # Necessary for backwards compatability with FALSE and TRUE inputs.
+  save.states <- as.character(save.states)   
+
+  # Handle default case where entire vector is passed; typically handled by 
+  # match.arg, but I want to allow for deprecated inputs: TRUE, FALSE, weighted, 
+  # unweighted. This condition once previous arguments are removed. 
+  if (length(save.states) > 1) {
+    save.states <- save.states[1]
+  }
+
+  state_arg_map <- c("FALSE" = "no", "TRUE" = "filter", "weighted" = "prediction", "unweighted" = "filter")
+
+  # Check if input is deprecated argument, and return warning. 
+  if (save.states %in% names(state_arg_map)) {
+    pWarn_(paste("The", sQuote("save.states"), "value", sQuote(save.states), "is deprecated and will be removed in a future version. Please use", sQuote(unname(state_arg_map[save.states])), "instead." ))
+    save.states <- unname(state_arg_map[save.states])
+  }
+
+  # Additional check to ensure final value is valid. 
+  save.states <- match.arg(save.states, c("no", "filter", "prediction"))
 
   params <- coef(object)
   times <- time(object,t0=TRUE)
@@ -247,8 +266,8 @@ pfilter_internal <- function (
   x <- init.x
 
   ## set up storage for saving samples from filtering distributions
-  stsav <- save.states %in% c("unweighted","TRUE")
-  wtsav <- save.states == "weighted"
+  stsav <- save.states %in% c("filter")
+  wtsav <- save.states == "prediction"
   if (stsav || wtsav || filter.traj) {
     xparticles <- vector(mode="list",length=ntimes)
     if (wtsav) xweights <- xparticles

R/saved_states.R

@@ -2,9 +2,10 @@
 ##'
 ##' Retrieve latent state trajectories from a particle filter calculation.
 ##'
-##' When one calls \code{\link{pfilter}} with \code{save.states=TRUE}, the latent state vector associated with each particle is saved.
+##' When one calls \code{\link{pfilter}} with \code{save.states="filter"} or \code{save.states="prediction"}, the latent state vector associated with each particle is saved.
 ##' This can be extracted by calling \code{saved_states} on the \sQuote{pfilterd.pomp} object.
-##' These are the \emph{unweighted} particles, saved \emph{after} resampling.
+##' If the filtered particles are saved, these particles are \emph{unweighted}, saved \emph{after} resampling using their normalized weights.
+##' If the argument  \code{save.states="prediction"} was used, the particles correspond to simulations from \code{rprocess}, and their corresponding unnormalized weights are included in the output. 
 ##'
 ##' @name saved_states
 ##' @aliases saved_states,ANY-method saved_states,missing-method

examples/pfilter.R

@@ -7,7 +7,7 @@ eff_sample_size(pf)             ## effective sample size
 logLik(pfilter(pf))      	## run it again with 1000 particles
 
 ## run it again with 2000 particles
-pf <- pfilter(pf,Np=2000,filter.mean=TRUE,filter.traj=TRUE,save.states="weighted")
+pf <- pfilter(pf,Np=2000,filter.mean=TRUE,filter.traj=TRUE,save.states="filter")
 fm <- filter_mean(pf) ## extract the filtering means
 ft <- filter_traj(pf) ## one draw from the smoothing distribution
 ss <- saved_states(pf,format="d") ## the latent-state portion of each particle

examples/traj_match.R

@@ -4,7 +4,7 @@
     traj_objfun(
       est=c("r","sigma","N_0"),
       partrans=parameter_trans(log=c("r","sigma","N_0")),
-      paramnames=c("r","sigma","N_0"),
+      paramnames=c("r","sigma","N_0")
       ) -> f
 
   f(log(c(20,0.3,10)))

examples/userdata.R

@@ -6,7 +6,7 @@
   ## C snippet approach:
 
   simulate(times=1:100,t0=0,
-    phi=as.double(100),
+    userdata=list(phi=as.double(100)),
     params=c(r=3.8,sigma=0.3,N.0=7),
     rprocess=discrete_time(
       step.fun=Csnippet(r"{
@@ -46,7 +46,7 @@
   ## Finally, the R function approach:
 
   simulate(times=1:100,t0=0,
-    phi=100,
+    userdata=list(phi=100),
     params=c(r=3.8,sigma=0.3,N_0=7),
     rprocess=discrete_time(
       step.fun=function (r, N, sigma, ...) {

tests/dacca.R

@@ -33,7 +33,7 @@ pfilter(
   pred.mean=TRUE,
   pred.var=TRUE,
   filter.traj=TRUE,
-  save.states=TRUE
+  save.states='filter'
 ) -> pf
 
 stopifnot(

tests/ebola.R

@@ -40,7 +40,7 @@ pfilter(
   pred.mean=TRUE,
   pred.var=TRUE,
   filter.traj=TRUE,
-  save.states=TRUE
+  save.states='filter'
 ) -> pf
 
 logLik(pf)

tests/gompertz.R

@@ -25,7 +25,7 @@ pfilter(
   pred.mean=TRUE,
   pred.var=TRUE,
   filter.traj=TRUE,
-  save.states=TRUE
+  save.states='filter'
 ) -> pf
 
 stopifnot(

tests/ou2.R

@@ -30,7 +30,7 @@ pfilter(
   pred.mean=TRUE,
   pred.var=TRUE,
   filter.traj=TRUE,
-  save.states=TRUE
+  save.states='filter'
 ) -> pf
 
 plot(pf,yax.flip=TRUE)

tests/pfilter.R

@@ -66,8 +66,8 @@ try(pfilter(pf,dmeasure=Csnippet("error(\"ouch!\");")))
 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="weighted")
+pf1 <- pfilter(pf,save.states='filter',filter.traj=TRUE)
+pf2 <- pfilter(pf,pred.mean=TRUE,pred.var=TRUE,filter.mean=TRUE,save.states="prediction")
 pf3 <- pfilter(pf,t0=1,filter.traj=TRUE)
 pf4 <- pfilter(pf,dmeasure=Csnippet("lik = (give_log) ? R_NegInf : 0;"),
   filter.traj=TRUE)
@@ -151,3 +151,18 @@ stopifnot(
   dim(y)==c(2,100),
   dim(y1)==c(1,100)
 )
+
+set.seed(1); pf_no <- pfilter(ou2, Np = 100, save.states = 'no')
+set.seed(1); pf_false <- pfilter(ou2, Np = 100, save.states = FALSE)
+set.seed(1); pf_true <- pfilter(ou2, Np = 100, save.states = TRUE)
+set.seed(1); pf_unweighted <- pfilter(ou2, Np = 100, save.states = 'unweighted')
+set.seed(1); pf_filter <- pfilter(ou2, Np = 100, save.states = 'filter')
+set.seed(1); pf_weighted <- pfilter(ou2, Np = 100, save.states = 'weighted')
+set.seed(1); pf_predict <- pfilter(ou2, Np = 100, save.states = 'prediction')
+stopifnot(
+  all.equal(unlist(saved_states(pf_unweighted)), unlist(saved_states(pf_filter))),
+  all.equal(unlist(saved_states(pf_true)), unlist(saved_states(pf_filter))),
+  all.equal(unlist(saved_states(pf_weighted)[2]), unlist(saved_states(pf_predict)[2])),
+  length(saved_states(pf_no)) == 0,
+  length(saved_states(pf_false)) == 0
+)

tests/ricker.R

@@ -25,7 +25,7 @@ pfilter(
   pred.mean=TRUE,
   pred.var=TRUE,
   filter.traj=TRUE,
-  save.states=TRUE
+  save.states="filter"
 ) -> pf
 
 stopifnot(

tests/rw2.R

@@ -30,7 +30,7 @@ pfilter(
   pred.mean=TRUE,
   pred.var=TRUE,
   filter.traj=TRUE,
-  save.states=TRUE
+  save.states="filter"
 ) -> pf
 
 plot(pf,yax.flip=TRUE)

tests/sir.R

@@ -28,7 +28,7 @@ pfilter(
   pred.mean=TRUE,
   pred.var=TRUE,
   filter.traj=TRUE,
-  save.states=TRUE
+  save.states="filter"
 ) -> pf
 
 stopifnot(

tests/sir2.R

@@ -28,7 +28,7 @@ pfilter(
   pred.mean=TRUE,
   pred.var=TRUE,
   filter.traj=TRUE,
-  save.states=TRUE
+  save.states="filter"
 ) -> pf
 
 stopifnot(

tests/steps.R

@@ -130,7 +130,7 @@ s1 |>
   )
 
 s2 |>
-  pfilter(Np=1,save.states=TRUE) |>
+  pfilter(Np=1,save.states='filter') |>
   saved_states(format="list") |>
   melt() |>
   pivot_wider() |>

tests/verhulst.R

@@ -25,7 +25,7 @@ pfilter(
   pred.mean=TRUE,
   pred.var=TRUE,
   filter.traj=TRUE,
-  save.states=TRUE
+  save.states='filter'
 ) -> pf
 
 print(logLik(pf))