[Robkalman-commits] r90 - in pkg/robKalman: . R demo man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu May 15 22:28:59 CEST 2025
Author: ruckdeschel
Date: 2025-05-15 22:28:58 +0200 (Thu, 15 May 2025)
New Revision: 90
Added:
pkg/robKalman/R/recEFilter.R
pkg/robKalman/R/recFilterWrapper.R
Removed:
pkg/robKalman/R/recFilter.R
pkg/robKalman/chm/
Modified:
pkg/robKalman/R/generate.R
pkg/robKalman/R/rLSfilter.R
pkg/robKalman/demo/rLSdemo.R
pkg/robKalman/man/0robKalman-package.Rd
pkg/robKalman/man/calibrateRLS.Rd
pkg/robKalman/man/internalACM.Rd
pkg/robKalman/man/internalKalman.Rd
pkg/robKalman/man/internalrLS.Rd
pkg/robKalman/man/recFilter.Rd
pkg/robKalman/man/simulateSScont.Rd
pkg/robKalman/man/util.Rd
Log:
delete chm files
Modified: pkg/robKalman/R/generate.R
===================================================================
--- pkg/robKalman/R/generate.R 2025-05-15 14:52:46 UTC (rev 89)
+++ pkg/robKalman/R/generate.R 2025-05-15 20:28:58 UTC (rev 90)
@@ -370,14 +370,16 @@
timestamps.x = timestamps.x))
}
-createFctCtrl <- function(fct, control) return(list(fct = fct, control = control))
+createFctCtrlPrepPost <- function(fct, control, prep, prep.ctl, post, post.ctl)
+ return(list(fct = fct, control = control,
+ prep = prep, prep.ctl = prep.ctl, post = post, post.ctl = post.ctl))
-createFilterProc <- function(init.fct, init.ctrl = NULL,
- pred.fct, pred.ctrl = NULL,
- corr.fct, corr.ctrl = NULL){
- return(list(init = createFctCtrl(init.fct, init.ctrl),
- pred = createFctCtrl(corr.fct, corr.ctrl)
- corr = createFctCtrl(pred.fct, pred.ctrl)))
+createFilterProc <- function(init.fct, init.ctrl = NULL, init.prep.fct = NULL, init.prep.ctl = NULL, init.post.fct = NULL, init.post.ctl = NULL,
+ pred.fct, pred.ctrl = NULL, pred.prep.fct = NULL, prep.prep.ctl = NULL, prep.post.fct = NULL, prep.post.ctl = NULL,
+ corr.fct, corr.ctrl = NULL, corr.prep.fct = NULL, corr.prep.ctl = NULL, corr.post.fct = NULL, corr.post.ctl = NULL ){
+ return(list(init = createFctCtrlPrepPost(init.fct, init.ctrl, init.prep.fct, init.prep.ctl, init.post.fct, init.post.ctl),
+ pred = createFctCtrlPrepPost(pred.fct, pred.ctrl, pred.prep.fct, prep.prep.ctl, prep.post.fct, prep.post.ctl),
+ corr = createFctCtrlPrepPost(corr.fct, corr.ctrl, corr.prep.fct, corr.prep.ctl, corr.post.fct, corr.post.ctl)))
}
createRobFilterProc <- function(init.fct.cla, init.ctrl.cla = NULL,
Modified: pkg/robKalman/R/rLSfilter.R
===================================================================
--- pkg/robKalman/R/rLSfilter.R 2025-05-15 14:52:46 UTC (rev 89)
+++ pkg/robKalman/R/rLSfilter.R 2025-05-15 20:28:58 UTC (rev 90)
@@ -3,6 +3,7 @@
################################################################
# P. Ruckdeschel 10.06.06
+# revised 15.08.11
#
## we assume a time-invariant SSM of the following form
Added: pkg/robKalman/R/recEFilter.R
===================================================================
--- pkg/robKalman/R/recEFilter.R (rev 0)
+++ pkg/robKalman/R/recEFilter.R 2025-05-15 20:28:58 UTC (rev 90)
@@ -0,0 +1,234 @@
+#######################################################
+##
+## recursive filter algorithm for Kalman filter routines
+## author: Bernhard Spangl & Peter Ruckdeschel
+## version: 0.3 (changed: 2011-08-16, created: 2011-06-21)
+##
+#######################################################
+
+recEFilter <- function(obs, model, proc)
+#
+#recEFilter <- function (Y, a, S, F, Q, Z, V,
+#% u=matrix(0, nrow=length(a), ncol=ncol(Y)),
+# w=matrix(0, nrow=nrow(Y), ncol=ncol(Y)),
+# v=matrix(0, nrow=length(a), ncol=ncol(Y)),
+# eps=matrix(0, nrow=nrow(Y), ncol=ncol(Y)),
+# R=NULL, T=NULL, exQ=NULL, exV=NULL,
+# initSc=.cEKFinitstep, predSc=.cEKFpredstep, corrSc=.cEKFcorrstep,
+# initSr=NULL, predSr=NULL, corrSr=NULL,
+# controlc=NULL, controlr=NULL,
+# controlF=NULL, controlQ=NULL, controlZ=NULL, controlV=NULL,
+# ...)
+{
+## a generalization of the extended Kalman filter, no vectorization!
+## arguments:
+## + Y : observations in a matrix with dimensions qd x tt
+## + a, S, F, Q, Z, V: Hyper-parameters of the ssm
+## + u, w : exogenous variables, ?? x tt matrix
+## + mu.v, mu.eps : expectation of innivations and observation noise
+## + R, T : selection matrices of innovations and observation
+## noise
+## + exQ, exV : exogenous variables of Q and V
+## + initSc, predSc, corrSc: (classical) initialization-, prediction-, and
+## correction-step function
+## + initSr, predSr, corrSr: (robust) initialization-, prediction-, and
+## correction-step function
+## + control_ : control paramaters, list
+## + ...: additional arguments for initS, predS, corrS
+
+# if (!(is.function(F))) createF(F=F, R=R)
+# if (!(is.function(Z))) createZ(Z=Z, T=T)
+# if (!(is.function(Q))) createQ(Q=Q)
+# if (!(is.function(V))) createV(V=V)
+
+
+ pd <- length(model$Start$a0)
+ qd <- (dim(obs$y))[1]
+ tt.y <- length(obs$timestamps.y)
+ tt.x <- length(obs$timestamps.x)-1
+
+ robust <- !(is.null(proc$robust$init) && is.null(proc$robust$pred) &&
+ is.null(proc$robust$corr))
+
+ Xf <- array(0, dim = c(pd, tt.x + 1))
+ Xp <- array(0, dim = c(pd, tt.x))
+ DeltaY <- array(0, dim = c(qd, tt.x))
+ St0 <- array(0, dim = c(pd, pd, tt.x + 1))
+ St1 <- array(0, dim = c(pd, pd, tt.x))
+ KG <- array(0, dim = c(pd, qd, tt.x))
+ Delta <- array(0, dim = c(qd, qd, tt.x))
+
+ if (robust) {
+ Xrf <- array(0, dim = c(pd, tt.x + 1))
+ Xrp <- array(0, dim = c(pd, tt.x))
+ DeltaYr <- array(0, dim = c(qd, tt.x))
+ Str0 <- array(0, dim = c(pd, pd, tt.x + 1))
+ Str1 <- array(0, dim = c(pd, pd, tt.x))
+ KGr <- array(0, dim = c(pd, qd, tt.x))
+ Deltar <- array(0, dim = c(qd, qd, tt.x))
+ }
+
+ ## initialization
+
+ ini <- proc$classic$init$fct(model$Start$a0, model$Start$Sigma0,
+ proc$classic$init$control, t = tt.x[1], ...)
+ x0 <- ini$x0
+ S0 <- ini$S0
+ controlc <- ini$control
+
+ Xf[, 1] <- ini$x0
+ St0[, , 1] <- ini$S0
+
+ if (robust) {
+ if (!is.null(proc$robust$init)) {
+ inir <- proc$robust$init$fct(model$Start$a0, model$Start$Sigma0,
+ proc$robust$init$control, t = tt.x[1], ...)
+ xr0 <- inir$x0
+ Sr0 <- inir$S0
+ controlr <- inir$control
+ } else {
+ xr0 <- x0
+ Sr0 <- S0
+ controlr <- controlc
+ }
+ Xrf[, 1] <- xr0
+ Str0[, , 1] <- Sr0
+ } else {
+ Xrf <- NULL
+ Xrp <- NULL
+ DeltaYr <- NULL
+ Str0 <- NULL
+ Str1 <- NULL
+ KGr <- NULL
+ Deltar <- NULL
+ }
+ j <- 1
+ for (i in (1:tt.x)) {
+
+
+
+ ## prediction
+ ps <- proc$classic$pred$fct(x0 = x0, S0 = S0, stateEq=model$StateEq,
+ mu.v=mu.v[,i], t = tt.x[i+1], i = i,
+ additinfofromCorr = controlc, ...)
+ x1 <- ps$x1
+ S1 <- ps$S1
+ controlc <- ps$additinfofromPred
+
+ Xp[, i] <- x1
+ St1[, , i] <- S1
+
+ if (robust) {
+ if (!is.null(proc$robust$pred)) {
+ psr <- proc$robust$pred$fct(x0 = xr0, S0 = Sr0, stateEq=model$StateEq,
+ mu.v=mu.v[,i], t = tt.x[i+1], i = i,
+ additinfofromCorr = controlr, ...)
+ } else {
+ psr <- proc$classic$pred$fct(x0 = xr0, S0 = Sr0, stateEq=model$StateEq,
+ mu.v=mu.v[,i], t = tt.x[i+1], i = i,
+ additinfofromCorr = controlr, ...)
+ }
+ xr1 <- psr$x1
+ Sr1 <- psr$S1
+ controlr <- psr$controlPred
+
+ Xrp[, i] <- xr1
+ Str1[, , i] <- Sr1
+ }
+
+ if (tt.x[i+1] %in% tt.y){
+ ## correction
+ cs <- proc$classic$corr(y = Y0[,j], x1 = x1, S1 = S1, obsEq=model$ObsEq,
+ mu.eps = mu.eps[, i], t = tt.x[i+1], i = i,
+ additinfofromPred = controlc, ...)
+ x0 <- cs$x0
+ S0 <- cs$S0
+ controlc <- cs$additinfofromCorr
+
+ Xf[, i + 1] <- x0
+ St0[, , i + 1] <- S0
+ KG[, , i] <- cs$K
+ Delta[, , i] <- cs$Delta
+ DeltaY[, i] <- cs$DeltaY
+
+ if (robust) {
+ if (!is.null(proc$robust$corr)) {
+ csr <- proc$robust$corr(y = Y0[,j], x1 = xr1, S1 = Sr1, obsEq=model$ObsEq,
+ mu.eps = mu.eps[, i], t = tt.x[i+1], i = i,
+ additinfofromPred = controlr, ...)
+ } else {
+ csr <- proc$classic$corr(y = Y0[,j], x1 = xr1, S1 = Sr1, obsEq=model$ObsEq,
+ mu.eps = mu.eps[, i], t = tt.x[i+1], i = i,
+ additinfofromPred = controlr, ...)
+ }
+ xr0 <- csr$x0
+ Sr0 <- csr$S0
+ controlr <- csr$controlCorr
+
+ Xrf[, i + 1] <- xr0
+ Str0[, , i + 1] <- Sr0
+ KGr[, , i] <- csr$K
+ Deltar[, , i] <- csr$Delta
+ DeltaYr[, i] <- csr$DeltaY
+ }
+ j <- j + 1
+ }
+ }
+ return(list(Xf = Xf, Xp = Xp, Xrf = Xrf, Xrp = Xrp,
+ S0 = St0, S1 = St1, KG = KG, Delta = Delta,
+ DeltaY = DeltaY,
+ Sr0 = Str0, Sr1 = Str1, KGr = KGr, Deltar = Deltar,
+ DeltaYr = DeltaYr))
+
+}
+
+
+#######################################################
+## simple wrappers:
+#######################################################
+
+ExtendedKF <- function (Y, a, S, F, Q, Z, V,
+ u=matrix(0, nrow=length(a), ncol=ncol(Y)),
+ w=matrix(0, nrow=nrow(Y), ncol=ncol(Y)),
+ v=matrix(0, nrow=length(a), ncol=ncol(Y)),
+ eps=matrix(0, nrow=nrow(Y), ncol=ncol(Y)),
+ R=NULL, T=NULL, exQ=NULL, exV=NULL,
+ controlF=NULL, controlQ=NULL, controlZ=NULL, controlV=NULL,
+ ...)
+{
+## arguments:
+## + Y : observations in a matrix with dimensions qd x tt
+## + a, S, F, Q, Z, V: Hyper-parameters of the ssm
+
+ recEFilter(Y, a, S, F, Q, Z, V,
+ u=u, w=w, v=v, eps=eps, R=R, T=T, exQ=exQ, exV=exV,
+ controlF=controlF, controlQ=controlQ,
+ controlZ=controlZ, controlV=controlV)
+
+}
+
+UnscentedKF <- function (Y, a, S, F, Q, Z, V,
+ u=matrix(0, nrow=length(a), ncol=ncol(Y)),
+ w=matrix(0, nrow=nrow(Y), ncol=ncol(Y)),
+ v=matrix(0, nrow=length(a), ncol=ncol(Y)),
+ eps=matrix(0, nrow=nrow(Y), ncol=ncol(Y)),
+ R=NULL, T=NULL, exQ=NULL, exV=NULL,
+ controlc=list(alpha=0.0001, beta=2, kappa=length(a)-3),
+ controlF=NULL, controlQ=NULL, controlZ=NULL, controlV=NULL,
+ ...)
+{
+## arguments:
+## + Y : observations in a matrix with dimensions qd x tt
+## + a, S, F, Q, Z, V: Hyper-parameters of the ssm
+
+ recEFilter(Y, a, S, F, Q, Z, V,
+ u=u, w=w, v=v, eps=eps, R=R, T=T, exQ=exQ, exV=exV,
+ initSc=.cUKFinitstep, predSc=.cUKFpredstep,
+ corrSc=.cUKFcorrstep,
+ controlc=controlc,
+ controlF=controlF, controlQ=controlQ,
+ controlZ=controlZ, controlV=controlV)
+
+}
+
+
Deleted: pkg/robKalman/R/recFilter.R
===================================================================
--- pkg/robKalman/R/recFilter.R 2025-05-15 14:52:46 UTC (rev 89)
+++ pkg/robKalman/R/recFilter.R 2025-05-15 20:28:58 UTC (rev 90)
@@ -1,426 +0,0 @@
-recursiveFilter <- function (Y, a, S, F, Q, Z, V,
- initSc = .cKinitstep, predSc = .cKpredstep,
- corrSc = .cKcorrstep,
- initSr = NULL, predSr = NULL, corrSr = NULL,
- nsim = 0, seed = NULL, ..., dropRuns = TRUE,
- CovRunDep = FALSE, saveOpt = TRUE, dimsCheck = NULL)
-# a generalization of the Kalmanfilter
-# arguments:
-# + Y : observations in an array with dimensions qd x runs x tt
-# (for backward compatibility: or a vector /a matrix in
-# dimensions qd x tt)
-# + a, S, F, Q, Z, V: Hyper-parameters of the ssm
-# + initSc, predSc, corrSc: (classical) initialization-, prediction-, and
-# correction-step function
-# + initSr, predSr, corrSr: (robust) initialization-, prediction-, and
-# correction-step function
-# + robustIO: if TRUE indicators are recorded whether prediction step does
-# clipping
-# + robustAO: if TRUE indicators are recorded whether correction step does
-# clipping
-# + nsim: if >0 we simulate a bunch of nsim paths (acc. to ideal model) to
-# get emp. covariances
-# + seed: seed for the simulations
-# + ...: additional arguments for initS, predS, corrS
-# + dropRuns: shall run-dimension be collapsed if it is one?
-# + CovRunDep: logical (default: FALSE), are there different prediction
-# and filter error covariances for each simulated path,
-# i.e., 'runs' > 1, as e.g. in the mACMfilter (-> complete
-# vectorization may not be possible!)
-# + saveOpt: logical (default: TRUE), should the stuff really be saved?
-# + dimsCheck: either 'NULL' (default) or vector [pd, qd, runs, tt] with
-# correct dimensions
-{
- if (is.null(dimsCheck)) {
- if (is.null(dim(Z)) && length(Z) > 1)
- stop("Error: Z has to be scalar or matrix!")
- if (!is.null(dim(Z))) qd <- dim(Z)[1]
- else qd <- 1
-
- pd <- length(a)
-
- ########################
- # for backward compatibility
- l0 <- length(dim(Y))
- if (l0 < 3) {
- if (l0 == 0) Y <- array(Y, dim = c(1, 1, length(Y)))
- if (l0 == 2) {
- if (dim(Y)[2] == qd) {
- Y <- aperm(array(Y, dim = c(dim(Y), 1)), c(1, 3, 2))
- } else {
- stop("Error: Dimensions of Z and/or Y do not fit!")
- }
- }
- }
- ########################
-
- tt <- (dim(Y))[3]
- runs <- (dim(Y))[2]
-
- if(pd==1) F <- array(F, dim = c(pd, pd, tt))
- if(pd==1 && qd==1) Z <- array(Z, dim = c(qd, pd, tt))
- if(pd==1) Q <- array(Q, dim = c(pd, pd, tt))
- if(qd==1) V <- array(V, dim = c(qd, qd, tt))
- if(is.matrix(F)) F <- array(F, dim = c(pd, pd, tt))
- if(is.matrix(Z)) Z <- array(Z, dim = c(qd, pd, tt))
- if(is.matrix(Q)) Q <- array(Q, dim = c(pd, pd, tt))
- if(is.matrix(V)) V <- array(V, dim = c(qd, qd, tt))
- } else {
- pd <- dimsCheck[1]
- qd <- dimsCheck[2]
- runs <- dimsCheck[3]
- tt <- dimsCheck[4]
- }
-
- WriteRecF <- WriteRecF(CovRunDep)
-
- saveOpt <- rep(saveOpt, length.out=8)
- if (length(names(saveOpt))==0) {
- names(saveOpt) <- c("KG", "KGr", "Delta", "Deltar",
- "DeltaY", "DeltaYr", "IndIO", "IndAO")
- }
-
- IndIO <- NULL
- IndAO <- NULL
- St0s <- St1s <- NULL
-
- robust <- !(is.null(initSr) && is.null(predSr) && is.null(corrSr))
-
- Xf <- array(0, dim = c(pd, runs, tt + 1))
- Xp <- array(0, dim = c(pd, runs, tt))
- St0 <- array(0, dim = c(pd, pd, tt + 1))
- St1 <- array(0, dim = c(pd, pd, tt))
- if (saveOpt["KG"]) KG <- array(0, dim = c(pd, qd, tt))
- else KG <- NULL
- if (saveOpt["Delta"]) Delta <- array(0, dim = c(qd, qd, tt))
- else Delta <- NULL
- if (saveOpt["DeltaY"]) DeltaY <- array(0, dim = c(qd, runs, tt))
- else DeltaY <- NULL
-
- if (nsim) {
- if (!exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE))
- runif(1)
- if (is.null(seed))
- RNGstate <- get(".Random.seed", envir = .GlobalEnv)
- else {
- R.seed <- get(".Random.seed", envir = .GlobalEnv)
- set.seed(seed)
- RNGstate <- structure(seed, kind = as.list(RNGkind()))
- on.exit(assign(".Random.seed", R.seed, envir = .GlobalEnv))
- }
- } else {
- RNGstate <- NULL
- }
-
- if (robust) {
- Xrf <- array(0, dim = c(pd, runs, tt + 1))
- Xrp <- array(0, dim = c(pd, runs, tt))
- if (!CovRunDep) {
- Str0 <- array(0, dim = c(pd, pd, tt + 1))
- Str1 <- array(0, dim = c(pd, pd, tt))
- if (saveOpt["KGr"]) KGr <- array(0, dim = c(pd, qd, tt))
- else KGr <- NULL
- if (saveOpt["Deltar"]) Deltar <- array(0, dim = c(qd, qd, tt))
- else Deltar <- NULL
- } else {
- Str0 <- array(0, dim = c(pd, pd, runs, tt + 1))
- Str1 <- array(0, dim = c(pd, pd, runs, tt))
- if (saveOpt["KGr"]) KGr <- array(0, dim = c(pd, qd, runs, tt))
- else KGr <- NULL
- if (saveOpt["Deltar"]) Deltar <- array(0, dim = c(qd, qd, runs, tt))
- else Deltar <- NULL
- }
- if (saveOpt["DeltaYr"]) DeltaYr <- array(0, dim = c(qd, runs, tt))
- else DeltaYr <- NULL
- }
-
- if (!is.null(predSr) && saveOpt["IndIO"])
- IndIO <- matrix(FALSE, runs, tt)
-
- if (!is.null(corrSr) && saveOpt["IndAO"])
- IndAO <- matrix(FALSE, runs, tt)
-
- #initialization
- ini <- initSc(a, S, i = 0, ...)
- x0 <- ini$x0
- S0 <- ini$S0
-
- Xf[, , 1] <- ini$x0
- St0[, , 1] <- ini$S0
-
- if (robust) {
- if (nsim) {
- Xs <- t(mvrnorm(nsim, a, S))
- St0s <- array(0, c(pd, pd, tt))
- St1s <- array(0, c(pd, pd, tt))
- }
- if (!is.null(initSr)) {
- inir <- initSr(a, S, i = 0, ...)
- xr0 <- inir$x0
- Sr0 <- inir$S0
- rob0 <- inir$rob
- } else {
- xr0 <- x0
- Sr0 <- S0
- rob0 <- NULL
- }
- Xrf[, , 1] <- xr0
- Str0 <- WriteRecF(Str0, Sr0, 1, opt=TRUE) # Str0[, , 1] <- Sr0
- rob0L <- list(rob0)
- } else {
- Xrf <- NULL
- Xrp <- NULL
- Str0 <- NULL
- Str1 <- NULL
- KGr <- NULL
- Deltar <- NULL
- DeltaYr <- NULL
- rob0L <- NULL
- rob1L <- NULL
- }
- for (i in (1:tt)) {
-
- #prediction
- F0 <- matrix(F[, , i], nrow = pd, ncol = pd)
- Q0 <- matrix(Q[, , i], nrow = pd, ncol = pd)
-
- ps <- predSc(x0 = x0, S0 = S0, F = F0, Q = Q0, i = i, ...)
- x1 <- matrix(ps$x1, nrow = pd, ncol = runs)
- S1 <- ps$S1
-
- Xp[, , i] <- x1
- St1[, , i] <- S1
-
- if (robust) {
- if (!is.null(predSr)) {
- psr <- predSr(x0 = xr0, S0 = Sr0, F = F0, Q = Q0, i = i,
- ..., rob0 = rob0)
- if (saveOpt["IndIO"]) IndIO[, i] <- as.logical(psr$Ind)
- if (nsim) {
- vs <- t(mvrnorm(nsim, a*0, Q0))
- Xs <- F0 %*% Xs + vs
- xr1s <- predSr(x0 = xr0s, S0 = Sr0, F = F0, Q = Q0, i = i,
- ..., rob0 = rob0)$x1
- St1s[, , i] <- cov(t(xr1s))
- }
- } else {
- psr <- predSc(x0 = xr0, S0 = Sr0, F = F0, Q = Q0, i = i, ...)
- if (nsim) {
- vs <- t(mvrnorm(nsim, a*0, Q0))
- Xs <- F %*% Xs + vs
- xr1s <- predSc(x0 = xr0s, S0 = Sr0, F = F0, Q = Q0, i = i,
- ...)$x1
- St1s[, , i] <- cov(t(xr1s))
- }
- }
- xr1 <- psr$x1
- Sr1 <- psr$S1
- rob1 <- psr$rob1
-
- Xrp[, , i] <- xr1
- Str1 <- WriteRecF(Str1, Sr1, i, opt=TRUE) # Str1[, , i] <- Sr1
- if (i==1) rob1L <- list(rob1)
- else rob1L[[i]] <- rob1
- }
-
- #correction
- Z0 <- matrix(Z[, , i], nrow = qd, ncol = pd)
- V0 <- matrix(V[, , i], nrow = qd, ncol = qd)
- Y0 <- matrix(Y[, , i], nrow = qd, ncol = runs)
-
- cs <- corrSc(y = Y0, x1 = x1, S1 = S1, Z = Z0, V = V0, i = i, ...)
- x0 <- cs$x0
- S0 <- cs$S0
-
- Xf[, , i + 1] <- x0
- St0[, , i + 1] <- S0
- if (saveOpt["KG"]) KG[, , i] <- cs$K
- if (saveOpt["Delta"]) Delta[, , i] <- cs$Delta
- if (saveOpt["DeltaY"]) DeltaY[, , i] <- cs$DeltaY
-
- if (robust) {
- if (!is.null(corrSr)) {
- csr <- corrSr(y = Y0, x1 = xr1, S1 = Sr1,
- Z = Z0, V = V0, i = i, ..., rob1 = rob1)
- if (saveOpt["IndAO"]) IndAO[, i] <- as.logical(csr$Ind)
- if (nsim) {
- es <- t(mvrnorm(nsim, Y[, 1]*0, V0))
- Ys <- Z0 %*% Xs + es
- xr0s <- corrSr(y = Ys, x1 = xr1s, S1 = Sr1,
- Z = Z0, V = V0, i = i, ..., rob1 = rob1)$x0
- St0s[, , i] <- cov(t(xr0s))
- }
- } else {
- csr <- corrSc(y = Y0, x1 = xr1, S1 = Sr1, Z = Z0, V = V0,
- i = i, ...)
- if (nsim) {
- es <- t(mvrnorm(nsim, Y[, 1]*0, V0))
- Ys <- Z0 %*% Xs + es
- xr0s <- corrSc(y = Ys, x1 = xr1s, S1 = Sr1,
- Z = Z0, V = V0, i = i, ...)$x0
- St0s[, , i] <- cov(t(xr0s))
- }
- }
- xr0 <- csr$x0
- Sr0 <- csr$S0
- rob0 <- csr$rob0
-
- Xrf[, , i + 1] <- xr0
- Str0 <- WriteRecF(Str0, Sr0, i+1, opt=TRUE)
- KGr <- WriteRecF(KGr, csr$K, i, opt=saveOpt["KGr"])
- Deltar <- WriteRecF(Deltar, csr$Delta, i, opt=saveOpt["Deltar"])
- # Str0[, , i + 1] <- Sr0
- # KGr[, , i] <- csr$K
- # Deltar[, , i] <- csr$Delta
- if (saveOpt["DeltaYr"]) DeltaYr[, , i] <- csr$DeltaY
- rob0L[[i + 1]] <- rob0
- }
- }
-
- if ((runs==1)&&(dropRuns)) {
- Xf <- matrix(Xf, pd, tt+1)
- Xp <- matrix(Xp, pd, tt)
- if (!is.null(Xrp)) {
- Xrf <- matrix(Xrf, pd, tt+1)
- Xrp <- matrix(Xrp, pd, tt)
- IndIO <- as.logical(IndIO)
- IndAO <- as.logical(IndAO)
- }
- }
- list(Xf = Xf, Xp = Xp, Xrf = Xrf, Xrp = Xrp,
- S0 = St0, S1 = St1, KG = KG, Delta = Delta,
- DeltaY = DeltaY,
- Sr0 = Str0, Sr1 = Str1,
- KGr = KGr, Deltar = Deltar, DeltaYr = DeltaYr,
- IndIO = IndIO, IndAO = IndAO,
- rob0L = rob0L, rob1L = rob1L,
- nsim = nsim, RNGstate = RNGstate,
- St0s = St0s, St1s = St1s)
-}
-
-
-######################################################
-# simple wrappers:
-######################################################
-
-KalmanFilter <- function(Y, a, S, F, Q, Z, V, dropRuns = TRUE)#
-#arguments:
-# + Y :observations
-# + a, S, F, Q, Z, V:Hyper-parameters of the ssm
-{recursiveFilter(Y, a, S, F, Q, Z, V, dropRuns = dropRuns)}
-
-rLSFilter <- function(Y, a, S, F, Q, Z, V, b, norm = EuclideanNorm, dropRuns = TRUE)#
-#arguments:
-# + Y :observations
-# + a, S, F, Q, Z, V:Hyper-parameters of the ssm
-# + b :clipping height
-{recursiveFilter(Y, a, S, F, Q, Z, V,
- initSc = .cKinitstep, predSc = .cKpredstep,
- corrSc = .cKcorrstep,
- #initSr=NULL, predSr=NULL,
- initSr = .cKinitstep, predSr = .cKpredstep,
- corrSr = .rLScorrstep, b = b, norm = norm, dropRuns = dropRuns)}
-
-##just a synonym for AO/SO robust filter
-rLS.AO.Filter <- rLSFilter
-
-## IO robust filter
-rLS.IO.Filter <- function(Y, a, S, F, Q, Z, V, b, norm = EuclideanNorm,
- dropRuns = TRUE)#
-#arguments:
-# + Y :observations
-# + a, S, F, Q, Z, V:Hyper-parameters of the ssm
-# + b :clipping height
-{recursiveFilter(Y, a, S, F, Q, Z, V,
- initSc = .cKinitstep, predSc = .cKpredstep,
- corrSc = .cKcorrstep,
- #initSr=NULL, predSr=NULL,
- initSr = .cKinitstep, predSr = .cKpredstep,
- corrSr = .rLS.IO.corrstep, b = b, norm = norm, dropRuns = dropRuns)}
-
-
-ACMfilter <- function(Y, a, S, F, Q, Z, V, s0, psi, apsi, bpsi, cpsi, flag, dropRuns = TRUE)#
-#arguments:
-# + Y :observations
-# + a, S, F, Q, Z, V:Hyper-parameters of the ssm
-# + b :clippingheight
-## Y=x ... observed time series
-## a=m0 ... unconditional mean
-## S=Cx ... covariance matrix
-## F=Phi ... design matrix of state equation
-## Q ... covariance matrix of state innovation process
-## Z=H ... observation matrix
-## V ... covariance matrix of observation noise
-## s0 ... scale of nominal Gaussian component of additive noise
-## psi ... influence function to be used (default: Hampel's psi function,
-## only that is available at the moment)
-## apsi, bpsi, cpsi ... tuning constants for Hampel's psi-function
-## (default: a=b=2.5, c=5.0)
-## flag ... character, if "weights" (default), use psi(t)/t to calculate
-## the weights; if "deriv", use psi'(t)
-{ recursiveFilter(Y, a, S, F, Q, Z, V,
- initSc = .cKinitstep, predSc = .cKpredstep,
- corrSc = .cKcorrstep,
- initSr = .cKinitstep, predSr = .ACMpredstep,
- corrSr = .ACMcorrstep, s0=s0, psi=psi,
- apsi=2.5, bpsi=2.5, cpsi=5.0, flag=flag, dropRuns = dropRuns)}
-###########################################
-##
-## R-function: ACMfilter - approximate conditional-mean filtering
-## author: Bernhard Spangl
-## version: 1.0 (2006-05-21)
-##
-###########################################
-
-## Paramters:
-## Y=x ... observed time series
-## a=m0 ... unconditional mean
-## S=Cx ... covariance matrix
-## F=Phi ... design matrix of state equation
-## Q ... covariance matrix of state innovation process
-## Z=H ... observation matrix
-## V ... covariance matrix of observation noise
-## s0 ... scale of nominal Gaussian component of additive noise
-## psi ... influence function to be used (default: Hampel's psi function,
-## only that is available at the moment)
-## a, b, c ... tuning constants for Hampel's psi-function
-## (defaul: a=b=2.5, c=5.0)
-## flag ... character, if "weights" (default), use psi(t)/t to calculate
-## the weights; if "deriv", use psi'(t)
-
-mACMfilter <- function(Y, a, S, F, Q, Z, V,
- psi=mvpsiHampel, apsi=2.5, bpsi=2.5, cpsi=5.0,
- flag="deriv", dropRuns = TRUE)
-{
-###########################################
-##
-## R-function: mACMfilter - approximate conditional-mean filtering
-## author: Bernhard Spangl
-## version: 0.2 (2008-03-31)
-##
-###########################################
-
-## Paramters:
-## Y ... observed vector-valued time series
-## (column-wise, matrix: q rows, number of columns equal to time points)
-## a ... unconditional mean vector (formerly: m0)
-## S ... covariance matrix (formerly: Cx)
-## F ... design matrix of state equation (formerly: Phi)
-## Q ... covariance matrix of state innovation process
-## Z ... observation matrix (formerly: H)
-## V ... covariance matrix of observation noise (formerly: R)
-## psi ... influence function to be used (default: Hampel's psi function,
-## only that is available at the moment)
-## apsi, bpsi, cpsi ... tuning constants for Hampel's psi-function
-## (default: a=b=2.5, c=5.0)
-## flag ... character, weight matrix to be used in correction step,
-## if "deriv" (default), Jacobian matrix of multivariate analogue
-## of Hampel's psi-function is used (only default is available
-## at the moment)
-
- recursiveFilter(Y, a, S, F, Q, Z, V,
- initSc=.cKinitstep, predSc=.cKpredstep, corrSc=.cKcorrstep,
- initSr=.cKinitstep, predSr=.cKpredstep, corrSr=.mACMcorrstep,
- psi, apsi, bpsi, cpsi, flag, dropRuns = dropRuns)
-}
-
Copied: pkg/robKalman/R/recFilterWrapper.R (from rev 41, pkg/robKalman/R/recFilter.R)
===================================================================
--- pkg/robKalman/R/recFilterWrapper.R (rev 0)
+++ pkg/robKalman/R/recFilterWrapper.R 2025-05-15 20:28:58 UTC (rev 90)
@@ -0,0 +1,426 @@
+recursiveFilter <- function (Y, a, S, F, Q, Z, V,
+ initSc = .cKinitstep, predSc = .cKpredstep,
+ corrSc = .cKcorrstep,
+ initSr = NULL, predSr = NULL, corrSr = NULL,
+ nsim = 0, seed = NULL, ..., dropRuns = TRUE,
+ CovRunDep = FALSE, saveOpt = TRUE, dimsCheck = NULL)
+# a generalization of the Kalmanfilter
+# arguments:
+# + Y : observations in an array with dimensions qd x runs x tt
+# (for backward compatibility: or a vector /a matrix in
+# dimensions qd x tt)
+# + a, S, F, Q, Z, V: Hyper-parameters of the ssm
+# + initSc, predSc, corrSc: (classical) initialization-, prediction-, and
+# correction-step function
+# + initSr, predSr, corrSr: (robust) initialization-, prediction-, and
+# correction-step function
+# + robustIO: if TRUE indicators are recorded whether prediction step does
+# clipping
+# + robustAO: if TRUE indicators are recorded whether correction step does
+# clipping
+# + nsim: if >0 we simulate a bunch of nsim paths (acc. to ideal model) to
+# get emp. covariances
+# + seed: seed for the simulations
+# + ...: additional arguments for initS, predS, corrS
+# + dropRuns: shall run-dimension be collapsed if it is one?
+# + CovRunDep: logical (default: FALSE), are there different prediction
+# and filter error covariances for each simulated path,
+# i.e., 'runs' > 1, as e.g. in the mACMfilter (-> complete
+# vectorization may not be possible!)
+# + saveOpt: logical (default: TRUE), should the stuff really be saved?
+# + dimsCheck: either 'NULL' (default) or vector [pd, qd, runs, tt] with
+# correct dimensions
+{
+ if (is.null(dimsCheck)) {
+ if (is.null(dim(Z)) && length(Z) > 1)
+ stop("Error: Z has to be scalar or matrix!")
+ if (!is.null(dim(Z))) qd <- dim(Z)[1]
+ else qd <- 1
+
+ pd <- length(a)
+
+ ########################
+ # for backward compatibility
+ l0 <- length(dim(Y))
+ if (l0 < 3) {
+ if (l0 == 0) Y <- array(Y, dim = c(1, 1, length(Y)))
+ if (l0 == 2) {
+ if (dim(Y)[2] == qd) {
+ Y <- aperm(array(Y, dim = c(dim(Y), 1)), c(1, 3, 2))
+ } else {
+ stop("Error: Dimensions of Z and/or Y do not fit!")
+ }
+ }
+ }
+ ########################
+
+ tt <- (dim(Y))[3]
+ runs <- (dim(Y))[2]
+
+ if(pd==1) F <- array(F, dim = c(pd, pd, tt))
+ if(pd==1 && qd==1) Z <- array(Z, dim = c(qd, pd, tt))
+ if(pd==1) Q <- array(Q, dim = c(pd, pd, tt))
+ if(qd==1) V <- array(V, dim = c(qd, qd, tt))
+ if(is.matrix(F)) F <- array(F, dim = c(pd, pd, tt))
+ if(is.matrix(Z)) Z <- array(Z, dim = c(qd, pd, tt))
+ if(is.matrix(Q)) Q <- array(Q, dim = c(pd, pd, tt))
+ if(is.matrix(V)) V <- array(V, dim = c(qd, qd, tt))
+ } else {
+ pd <- dimsCheck[1]
+ qd <- dimsCheck[2]
+ runs <- dimsCheck[3]
+ tt <- dimsCheck[4]
+ }
+
+ WriteRecF <- WriteRecF(CovRunDep)
+
+ saveOpt <- rep(saveOpt, length.out=8)
+ if (length(names(saveOpt))==0) {
+ names(saveOpt) <- c("KG", "KGr", "Delta", "Deltar",
+ "DeltaY", "DeltaYr", "IndIO", "IndAO")
+ }
+
+ IndIO <- NULL
+ IndAO <- NULL
+ St0s <- St1s <- NULL
+
+ robust <- !(is.null(initSr) && is.null(predSr) && is.null(corrSr))
+
+ Xf <- array(0, dim = c(pd, runs, tt + 1))
+ Xp <- array(0, dim = c(pd, runs, tt))
+ St0 <- array(0, dim = c(pd, pd, tt + 1))
+ St1 <- array(0, dim = c(pd, pd, tt))
+ if (saveOpt["KG"]) KG <- array(0, dim = c(pd, qd, tt))
+ else KG <- NULL
+ if (saveOpt["Delta"]) Delta <- array(0, dim = c(qd, qd, tt))
+ else Delta <- NULL
+ if (saveOpt["DeltaY"]) DeltaY <- array(0, dim = c(qd, runs, tt))
+ else DeltaY <- NULL
+
+ if (nsim) {
+ if (!exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE))
+ runif(1)
+ if (is.null(seed))
+ RNGstate <- get(".Random.seed", envir = .GlobalEnv)
+ else {
+ R.seed <- get(".Random.seed", envir = .GlobalEnv)
+ set.seed(seed)
+ RNGstate <- structure(seed, kind = as.list(RNGkind()))
+ on.exit(assign(".Random.seed", R.seed, envir = .GlobalEnv))
+ }
+ } else {
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/robkalman -r 90
More information about the Robkalman-commits
mailing list