Merge pull request #106 from nwfsc-timeseries/MARSSkfss

MARSSkfss
atsa-es · Sep 14, 2020 · c3e6f57 · c3e6f57
2 parents 88e5a64 + c39642d
commit c3e6f57
Show file tree

Hide file tree

Showing 6 changed files with 58 additions and 32 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,8 +1,8 @@
 Package: MARSS
 Type: Package
 Title: Multivariate Autoregressive State-Space Modeling
-Version: 3.11.2
-Date: 2020-09-12
+Version: 3.11.3
+Date: 2020-09-13
 Depends: R (>= 3.5.0)
 Imports: 
     graphics, grDevices, KFAS (>= 1.0.1), mvtnorm, nlme, stats, utils

diff --git a/R/MARSS.R b/R/MARSS.R
@@ -125,29 +125,29 @@ MARSS <- function(y,
       return(MLEobj)
     }
 
+    # Check that Kalman filter/smoother will run
     if (MLEobj$control$trace != -1) {
       MLEobj.test <- MLEobj
       MLEobj.test$par <- MLEobj$start
-      # Do this test with MARSSkfss since it has internal tests for problems
-      kftest <- try(MARSSkfss(MLEobj.test), silent = TRUE)
+      kftest <- try(MARSSkf(MLEobj.test), silent = TRUE)
       if (inherits(kftest, "try-error")) {
-        cat("Error: Stopped in MARSS() before fitting because MARSSkfss stopped.  Something in the model structure prevents the kfss Kalman filter running.\n Try using control$trace=-1 to turn off error-checking (and MARSSkfss calls) or use fit=FALSE and check the model you are trying to fit.\n")
+        cat("Error: Stopped in MARSS() before fitting because", fun.kf, "stopped.  Something in the model structure prevents the Kalman filter or smoother running.\n Try setting fun.kf to use a different KF function (MARSSkfss or MARSSkfas) or use fit=FALSE and check the model you are trying to fit. You can also try trace=1 to get more progress output.\n")
         MLEobj$convergence <- 2
         return(MLEobj.test)
       }
       if (!kftest$ok) {
         cat(kftest$msg)
-        cat(paste("Error: Stopped in MARSS() before fitting because MARSSkfss stopped.  Something in the model structure prevents the kfss Kalman filter running.\n Try using control$trace=-1 to turn off error-checking (and MARSSkfss calls) or use fit=FALSE and check the model you are trying to fit.\n", kftest$errors, "\n", sep = ""))
+        cat(paste("Error: Stopped in MARSS() before fitting because", fun.kf, "stopped.  Something in the model structure prevents the Kalman filter or smoother running.\n Try setting fun.kf to use a different KF function (MARSSkfss or MARSSkfas) or use fit=FALSE and check the model you are trying to fit. You can also try trace=1 to get more progress output.\n", kftest$errors, "\n", sep = ""))
         MLEobj$convergence <- 2
         return(MLEobj.test)
       }
       MLEobj.test$kf <- kftest
     }
     # Ey is needed for method=kem
-    if (MLEobj$control$trace != -1 & MLEobj$method == "kem") {
+    if (MLEobj$control$trace != -1 && MLEobj$method %in% kem.methods) {
       Eytest <- try(MARSShatyt(MLEobj.test), silent = TRUE)
       if (inherits(Eytest, "try-error")) {
-        cat("Error: Stopped in MARSS() before fitting because MARSShatyt stopped.  Something is wrong with the model structure that prevents MARSShatyt running.\n\n")
+        cat("Error: Stopped in MARSS() before fitting because MARSShatyt() stopped.  Something is wrong with the model structure that prevents MARSShatyt() running.\n\n")
         MLEobj.test$convergence <- 2
         MLEobj.test$Ey <- Eytest
         return(MLEobj.test)
@@ -159,7 +159,6 @@ MARSS <- function(y,
         MLEobj.test$Ey <- Eytest
         return(MLEobj.test)
       }
-      # MLEobj$Ey=Eytest
     }
 
     if (!fit) MLEobj$convergence <- 3
@@ -185,10 +184,11 @@ MARSS <- function(y,
         MLEobj <- MARSSaic(MLEobj)
         MLEobj$coef <- coef(MLEobj, type = "vector")
       }
+
       ## Add states.se and ytT.se if no errors.  Return kf and Ey if trace>0
-      if ((MLEobj$convergence %in% c(0, 1, 3)) || (MLEobj$convergence %in% c(10, 11) && method %in% kem.methods)) {
+      if (MLEobj$convergence %in% c(0, 1, 3) || (MLEobj$convergence %in% c(10, 11) && MLEobj$method %in% kem.methods)) {
         if (silent == 2) cat("Adding states and states.se.\n")
-        kf <- MARSSkf(MLEobj) # use function requested by user
+        kf <- MARSSkf(MLEobj) # use function requested by user; default smoother=TRUE
         MLEobj$states <- kf$xtT
         MLEobj$logLik <- kf$logLik
         if (!is.null(kf[["VtT"]])) {
@@ -222,23 +222,40 @@ MARSS <- function(y,
           ytT.se <- NULL
         }
         MLEobj$ytT.se <- ytT.se
-        if (MLEobj$control$trace > 0) { # then return kf and Ey
-          if (fun.kf == "MARSSkfas") kfss <- MARSSkfss(MLEobj) else kfss <- kf
+
+        # Return kf and Ey, if trace = 1
+        if (MLEobj$control$trace > 0) {
+          if (fun.kf == "MARSSkfas") {
+            kfss <- try(MARSSkfss(MLEobj, smoother = FALSE), silent = TRUE)
+            if (inherits(kfss, "try-error") || !kfss$ok) {
+              msg <- "Not available. MARSSkfss() returned error."
+              kfss <- list(Innov = msg, Sigma = msg, xtt = msg, Vtt = msg, J = msg, Kt = msg)
+            }
+          } else {
+            kfss <- kf
+          }
           MLEobj$kf <- kf # from above will use function requested by user
-          # except these are only returned by MARSSkfss
+          MLEobj$Ey <- Ey # from above
+          # these are only returned by MARSSkfss
           MLEobj$Innov <- kfss$Innov
           MLEobj$Sigma <- kfss$Sigma
           MLEobj$Vtt <- kfss$Vtt
           MLEobj$xtt <- kfss$xtt
           MLEobj$J <- kfss$J
-          MLEobj$J0 <- kfss$J0
           MLEobj$Kt <- kfss$Kt
-          MLEobj$Ey <- MARSShatyt(MLEobj)
+          if (fun.kf == "MARSSkfss") {
+            MLEobj$J0 <- kfss$J0
+          } else {
+            # From kfss smoother so won't be available if fun.kf=MARSSkfas
+            J0 <- try(MARSSkfss(MLEobj), silent = TRUE)
+            if (!inherits(J0, "try-error") && J0$ok) MLEobj$J0 <- J0$J0 else MLEobj$J0 <- "Not available. MARSSkfss() smoother returned error."
+          }
         }
         # apply X and Y names various X and Y related elements
         MLEobj <- MARSSapplynames(MLEobj)
       }
-    } # fit the model
+    }
+    # END fit the model #################
 
     if ((!silent || silent == 2) & MLEobj$convergence %in% c(0, 1, 3, 10, 11, 12)) {
       print(MLEobj)

diff --git a/R/MARSSkf.r b/R/MARSSkf.r
@@ -2,7 +2,7 @@
 #   MARSSkf function
 #   Utility function to choose the Kalman filter and smoother
 #######################################################################################################
-MARSSkf <- function(MLEobj, only.logLik = FALSE, return.lag.one = TRUE, return.kfas.model = FALSE, newdata = NULL) {
+MARSSkf <- function(MLEobj, only.logLik = FALSE, return.lag.one = TRUE, return.kfas.model = FALSE, newdata = NULL, smoother = TRUE) {
   if (is.null(MLEobj$par)) {
     stop("Stopped in MARSSkf(): par element of marssMLE object is required.\n")
   }
@@ -18,7 +18,7 @@ MARSSkf <- function(MLEobj, only.logLik = FALSE, return.lag.one = TRUE, return.k
     MLEobj[["marss"]][["data"]] <- newdata
   }
   if (MLEobj$fun.kf == "MARSSkfss") {
-    return(MARSSkfss(MLEobj))
+    return(MARSSkfss(MLEobj, smoother = smoother))
   }
   if (MLEobj$fun.kf == "MARSSkfas") {
     return(MARSSkfas(MLEobj, only.logLik = only.logLik, return.lag.one = return.lag.one, return.kfas.model = return.kfas.model))

diff --git a/R/MARSSkfas.r b/R/MARSSkfas.r
@@ -237,7 +237,7 @@ MARSSkfas <- function(MLEobj, only.logLik = FALSE, return.lag.one = TRUE, return
 
   if (!return.kfas.model) kfas.model <- NULL
 
-  # not using ks.out$v (Innov) and ks.out$F (Sigma) since I think there might be a bug (I misunderstand KFAS) when R is not diagonal.
+  # not using ks.out$v (Innov) and ks.out$F (Sigma) since I think there might be a bug (or I misunderstand KFS output) when R is not diagonal.
   rtn.list <- list(
     xtT = ks.out$alphahat[1:m, , drop = FALSE],
     VtT = VtT,

diff --git a/R/MARSSkfss.R b/R/MARSSkfss.R
@@ -4,7 +4,7 @@
 #   ** All eqn refs are to 2nd ed of Shumway & Stoffer (2006): Time Series Analysis and Its Applications
 # 5-17-12  I removed the tests re inversion of Z when R=0 and put that in MARSSkem
 #######################################################################################################
-MARSSkfss <- function(MLEobj) {
+MARSSkfss <- function(MLEobj, smoother=TRUE) {
   condition.limit <- 1E10
   condition.limit.Ft <- 1E5 # because the Ft is used to compute LL and LL drop limit is about 2E-8
 
@@ -268,6 +268,7 @@ MARSSkfss <- function(MLEobj) {
   ######################################################
   # BACKWARD PASS (Kalman smoother) gets you E[x(t)|y(1:T)] from E[x(t)|y(1:t)]
   ######################################################
+  if (smoother){
   xtT[, TT] <- xtt[, TT, drop = FALSE]
   VtT[, , TT] <- Vtt[, , TT]
   # indexing is 0 to T for the backwards smoother recursions
@@ -406,8 +407,11 @@ MARSSkfss <- function(MLEobj) {
     Vtt1T[, , 1] <- Vtt[, , 1] %*% t(J0) + J[, , 1] %*% (Vtt1T[, , 2] - B %*% Vtt[, , 1]) %*% t(J0)
   }
   if (init.state == "x10") Vtt1T[, , 1] <- NA
-
-  ###########################################################
+  } else {
+  x0T <- V0T <- VtT <- xtT <- J0 <- Vtt1T <- NULL
+  }
+  ## END SMOOTHER #########################################################
+
   rtn.list <- list(
     xtT = xtT, VtT = VtT, Vtt1T = Vtt1T, x0T = x0T, V0T = V0T, Vtt = Vtt,
     Vtt1 = Vtt1, J = J, J0 = J0, Kt = Kt, xtt1 = xtt1, xtt = xtt, Innov = vt, Sigma = Ft
@@ -417,7 +421,11 @@ MARSSkfss <- function(MLEobj) {
   MODELobj <- MLEobj[["marss"]]
   X.names <- attr(MODELobj, "X.names")
   Y.names <- attr(MODELobj, "Y.names")
-  for (el in c("xtT", "xtt1", "xtt", "x0T")) rownames(rtn.list[[el]]) <- X.names
+  if (smoother){
+    for (el in c("xtT", "xtt1", "xtt", "x0T")) rownames(rtn.list[[el]]) <- X.names
+  } else {
+    for (el in c("xtt1", "xtt")) rownames(rtn.list[[el]]) <- X.names
+  }
   rownames(rtn.list[["Innov"]]) <- Y.names
 
   # Calculate log likelihood, see eqn 6.62
@@ -431,11 +439,11 @@ MARSSkfss <- function(MLEobj) {
           ok = FALSE, logLik = NaN,
           errors = paste("One of the diagonal elements of Sigma[,,", t, "]=0. That should never happen when t>1 or t=1 and tinitx=0.  \n Are both Q[i,i] and R[i,i] being set to 0?\n", sep = "")
         )))
-      } else { # t=1 so ok. get the det of Ft and deal with 0s that might appear on diag of Ft when t=1 and V0=0 and R=0 and tinitx=1
-        # 2-5-15 This isn't an error.  x10 can be != y[1] when R=0; x11 cannot be.
-        #          if(any(abs(vt[diag.Ft==0,1])>1E-16)){
-        #              return(c(rtn.list,list(ok=FALSE, logLik = -Inf, errors = "V0=0, tinitx=1, R=0 and y[1] does not match x0. You shouldn't estimate x0 when R=0.\n")))
-        #          }
+      } else { 
+        # t=1 so ok. get the det of Ft and deal with 0s that might appear on 
+        #   diag of Ft when t=1 and V0=0 and R=0 and tinitx=1
+        #   Note, x10 can be != y[1] when R=0; x11 cannot be.
+
         OmgF1 <- makediag(1, n)[diag.Ft != 0, , drop = FALSE] # permutation matrix
         # need to remove those y[1] associated with Ft[,,1]==0 that were non-missing
         loglike <- loglike + sum(diag.Ft == 0 & YM[, 1] != 0) / 2 * log(2 * base::pi)
@@ -451,7 +459,7 @@ MARSSkfss <- function(MLEobj) {
             detFt <- det(OmgF1 %*% tcrossprod(Ft[, , t], OmgF1))
           }
         }
-        # get the inv of Ft
+        # get the inverse of Ft
         if (n == 1) {
           Ftinv <- pcholinv(matrix(Ft[, , t], 1, 1))
         } else {

diff --git a/man/MARSSkf.Rd b/man/MARSSkf.Rd
@@ -8,12 +8,13 @@
 }
 \usage{
 MARSSkf(MLEobj, only.logLik=FALSE, return.lag.one=TRUE, return.kfas.model=FALSE, 
-           newdata=NULL)
-MARSSkfss(MLEobj)
+           newdata=NULL, smoother=TRUE)
+MARSSkfss(MLEobj, smoother=TRUE)
 MARSSkfas(MLEobj, only.logLik=FALSE, return.lag.one=TRUE, return.kfas.model=FALSE)
 }
 \arguments{
   \item{ MLEobj }{ A \code{\link{marssMLE}} object with the \code{par} element of estimated parameters, \code{marss} element with the model description (in marss form) and data, and \code{control} element for the fitting algorithm specifications.  \code{control$debugkf} specifies that detailed error reporting will be returned (only used by \code{MARSSkf}).  \code{model$diffuse=TRUE} specifies that a diffuse prior be used (only used by \code{MARSSkfas}). See \link[KFAS]{KFS} documentation. When the diffuse prior is set, V0 should be non-zero since the diffuse prior variance is V0*kappa, where kappa goes to infinity.}
+  \item{ smoother }{ Used by \code{MARSSkfss}.  If set to FALSE, only the Kalman filter is run. xtT, VtT, x0T, Vtt1T, V0T, and J0 will be NULL. }
   \item{ only.logLik }{ Used by \code{MARSSkfas}.  If set, only the log-likelihood is returned using the \link[KFAS]{KFAS} package function \code{\link[KFAS]{logLik.SSModel}}.  This is much faster if only the log-likelihood is needed. }
   \item{ return.lag.one }{ Used by \code{MARSSkfas}.  If set to FALSE, the smoothed lag-one covariance  values are not returned (Vtt1T is set to NULL).  This speeds up \code{MARSSkfas} because to return the smoothed lag-one covariance a stacked MARSS model is used with twice the number of state vectors---thus the state matrices are larger and take more time to work with. }
   \item{ return.kfas.model }{ Used by \code{MARSSkfas}.  If set to TRUE, it returns the MARSS model in \link[KFAS]{KFAS} model form (class \code{\link[KFAS]{SSModel}}).  This is useful if you want to use other KFAS functions or write your own functions to work with \code{\link{optim}} to do optimization.  This can speed things up since there is a bit of code overhead in \code{\link{MARSSoptim}} associated with the \code{\link{marssMODEL}} model specification needed for the constrained EM algorithm (but not strictly needed for \code{\link{optim}}; useful but not required.). }