From e25c5464f60cdb8941eb59777510bd042e74c382 Mon Sep 17 00:00:00 2001 From: EEH Date: Sun, 13 Sep 2020 18:27:22 -0700 Subject: [PATCH 1/3] make the error-checks a bit less picky if fun.kf=MARSSkfas, then MARSSkfss doesn't need to run. If trace=1, then run MARSSkfss filter only to get xtt etc. --- DESCRIPTION | 4 ++-- R/MARSS.R | 39 ++++++++++++++++++++++++--------------- R/MARSSkf.r | 4 ++-- R/MARSSkfas.r | 2 +- R/MARSSkfss.R | 28 ++++++++++++++++++---------- man/MARSSkf.Rd | 5 +++-- 6 files changed, 50 insertions(+), 32 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 3d075dba..a42544a6 100644 --- 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 index 0491bde3..18bd4bca 100644 --- 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,32 @@ 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 <- MARSSkfss(MLEobj, smoother=FALSE) 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"){ + # From kfss smoother so won't necessarily be available if fun.kf=MARSSkfas + MLEobj$J0 <- kfss$J0 + } else { + 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 index 93b6fdd3..ecb80630 100644 --- 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 index a5884cd0..5ffebdde 100644 --- 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 index 6de93211..b2c592b7 100644 --- 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 index 231cec71..237686b9 100644 --- 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.). } From 1cf6a2dc8cd02bdb8d157cec2bd0609357b2a6d9 Mon Sep 17 00:00:00 2001 From: EEH Date: Sun, 13 Sep 2020 18:35:58 -0700 Subject: [PATCH 2/3] make a little more robust to MARSSkfss Don't want trace=1 to fail just because MARSSkfss() won't run. --- R/MARSS.R | 28 ++++++++++++++++++---------- 1 file changed, 18 insertions(+), 10 deletions(-) diff --git a/R/MARSS.R b/R/MARSS.R index 18bd4bca..d70efe13 100644 --- a/R/MARSS.R +++ b/R/MARSS.R @@ -131,7 +131,7 @@ MARSS <- function(y, MLEobj.test$par <- MLEobj$start kftest <- try(MARSSkf(MLEobj.test), silent = TRUE) if (inherits(kftest, "try-error")) { - 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") + 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) } @@ -144,7 +144,7 @@ MARSS <- function(y, MLEobj.test$kf <- kftest } # Ey is needed for method=kem - if (MLEobj$control$trace != -1 && MLEobj$method %in% kem.methods) { + 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") @@ -184,9 +184,9 @@ 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) && MLEobj$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; default smoother=TRUE MLEobj$states <- kf$xtT @@ -222,10 +222,18 @@ MARSS <- function(y, ytT.se <- NULL } MLEobj$ytT.se <- ytT.se - + # Return kf and Ey, if trace = 1 - if (MLEobj$control$trace > 0) { - if (fun.kf == "MARSSkfas") kfss <- MARSSkfss(MLEobj, smoother=FALSE) else kfss <- kf + 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 MLEobj$Ey <- Ey # from above # these are only returned by MARSSkfss @@ -235,10 +243,10 @@ MARSS <- function(y, MLEobj$xtt <- kfss$xtt MLEobj$J <- kfss$J MLEobj$Kt <- kfss$Kt - if (fun.kf == "MARSSkfss"){ - # From kfss smoother so won't necessarily be available if fun.kf=MARSSkfas + 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." } @@ -246,7 +254,7 @@ MARSS <- function(y, # apply X and Y names various X and Y related elements MLEobj <- MARSSapplynames(MLEobj) } - } + } # END fit the model ################# if ((!silent || silent == 2) & MLEobj$convergence %in% c(0, 1, 3, 10, 11, 12)) { From c39642dde6d75cec9ea88cd438e269312c8b9fe9 Mon Sep 17 00:00:00 2001 From: EEH Date: Sun, 13 Sep 2020 18:44:36 -0700 Subject: [PATCH 3/3] typo --- R/MARSS.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/MARSS.R b/R/MARSS.R index d70efe13..898e9513 100644 --- a/R/MARSS.R +++ b/R/MARSS.R @@ -227,7 +227,7 @@ MARSS <- function(y, if (MLEobj$control$trace > 0) { if (fun.kf == "MARSSkfas") { kfss <- try(MARSSkfss(MLEobj, smoother = FALSE), silent = TRUE) - if (inherits(kfss, "try-error" || !kfss$ok)) { + 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) }