[Robkalman-commits] r40 - branches/robKalman_bs/pkg/robKalman/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Sep 3 11:35:14 CEST 2009
Author: bspangl
Date: 2009-09-03 11:35:13 +0200 (Thu, 03 Sep 2009)
New Revision: 40
Modified:
branches/robKalman_bs/pkg/robKalman/R/ACMfilt.R
branches/robKalman_bs/pkg/robKalman/R/ACMfilter.R
branches/robKalman_bs/pkg/robKalman/R/Psi.R
branches/robKalman_bs/pkg/robKalman/R/mACMfilter.R
branches/robKalman_bs/pkg/robKalman/R/recFilter.R
Log:
Updated branch 'branches/robKalman_bs/pkg/robKalman/R/' (Bernhard Spangl, 2009-09-03)
Modified: branches/robKalman_bs/pkg/robKalman/R/ACMfilt.R
===================================================================
--- branches/robKalman_bs/pkg/robKalman/R/ACMfilt.R 2009-07-23 19:51:16 UTC (rev 39)
+++ branches/robKalman_bs/pkg/robKalman/R/ACMfilt.R 2009-09-03 09:35:13 UTC (rev 40)
@@ -6,7 +6,7 @@
##
## R-function: ACMfilt - approximate conditional-mean filtering (wrapper)
## author: Bernhard Spangl
-## version: 1.1 (2007-08-13 and 2006-08-31)
+## version: 1.2 (2009-07-30)
## References:
## [Mart79c] R.D. Martin, Approximate Conditional-mean Type Smoothers
## and Interpolators (1979)
@@ -17,7 +17,7 @@
###########################################
## Paramters:
-## x ... univariate time series (vector)
+## x ... univariate time series (vector) or runs x tt matrix
## gm ... list as produced by function 'arGM' which includes components
## 'ar' containing the AR(p) coefficient estimates, 'sinnov' containing
## innovation scale estiamtes from AR(p) fits of orders 1 through p;
@@ -34,48 +34,54 @@
## if FALSE filtering from the top of ^X_t is performed
## Variable definitions:
-
- N <- length(x)
+ if (is.null(dim(x))) {
+ runs <- 1
+ tt <- length(x)
+ } else {
+ runs <- dim(x)[1]
+ tt <- dim(x)[2]
+ }
+ Y <- array(x, dim=c(1, dim(x)))
phi <- gm$ar
- p <- length(phi)
- si <- gm$sinnov[p]
+ pd <- length(phi)
+ m0 <- rep(0, pd)
Cx <- gm$Cx
- Phi <- cbind(rbind(phi[-p], diag(rep(1, (p-1)))), c(phi[p], rep(0, (p-1))))
- Q <- matrix(0, p, p)
-## Q <- diag(rep(0, p))
- Q[1, 1] <- si^2
+ Phi <- cbind(rbind(phi[-pd], diag(rep(1, (pd-1)))),
+ c(phi[pd], rep(0, (pd-1))))
+ F <- array(Phi, dim=c(pd, pd, tt))
+ si <- gm$sinnov[pd]
+ Q <- array(0, dim=c(pd, pd, tt))
+ Q[1, 1, ] <- si^2
- m0 <- rep(0, p)
- H <- matrix(c(1, rep(0, (p-1))), 1, p)
- V <- matrix(s0^2)
+ Z <- array(0, dim=c(1, pd, tt))
+ Z[1, 1, ] <- 1
+ V <- array(s0^2, dim=c(1, 1, tt))
psi <- .psi(psi)
-
+
## Centering:
x <- x - gm$mu
- ACMres <- ACMfilter(Y=array(x,dim=c(1,1,N)), a=m0, S=Cx,
- F=Phi, Q=Q, Z=H, V=V, s0=s0,
- psi=psi, apsi=a, bpsi=b, cpsi=c,
- flag=flag)
+ ACMres <- ACMfilter(Y=Y, a=m0, S=Cx, F=F, Q=Q, Z=Z, V=V,
+ psi=psi, apsi=a, bpsi=b, cpsi=c, flag=flag,
+ dropRuns=FALSE,
+ saveOpt=FALSE, dimsCheck=c(pd, 1, runs, tt))
- X.ck <- ACMres$Xf; X.ck <- X.ck[,2:(N+1)]
- X <- ACMres$Xrf; X <- X[,2:(N+1)]
- st <- as.numeric(unlist(ACMres$rob1L))
+ X.ck <- ACMres$Xf[, , 2:(tt+1)]
+ X <- ACMres$Xrf[, , 2:(tt+1)]
+ st <- matrix(unlist(ACMres$rob0L), runs, tt)
-
if (!lagsmo) {
- x.ck <- X.ck[1, ]
- x <- X[1, ]
+ x.ck <- X.ck[1, , ]
+ x <- X[1, , ]
} else {
- x.ck <- c(X.ck[p, p:N], X.ck[(p-1):1, N])
- x <- c(X[p, p:N], X[(p-1):1, N])
+ if (runs==1) {
+ x.ck <- c(X.ck[pd, , pd:tt], t(X.ck[(pd-1):1, , tt]))
+ x <- c(X[pd, , pd:tt], t(X[(pd-1):1, , tt]))
+ } else {
+ x.ck <- cbind(X.ck[pd, , pd:tt], t(X.ck[(pd-1):1, , tt]))
+ x <- cbind(X[pd, , pd:tt], t(X[(pd-1):1, , tt]))
+ }
}
-## ARmodel <- .ARmodel(x, p)
-## y <- ARmodel$y
-## Z <- ARmodel$Z
-## r <- resid(lm.fit(Z, y))
-
- return(list(filt.ck=x.ck +gm$mu, filt=x + gm$mu, st=st)) #,
-## r=c(rep(NA, p), r)))
+ return(list(filt.ck=x.ck +gm$mu, filt=x + gm$mu, st=st))
}
Modified: branches/robKalman_bs/pkg/robKalman/R/ACMfilter.R
===================================================================
--- branches/robKalman_bs/pkg/robKalman/R/ACMfilter.R 2009-07-23 19:51:16 UTC (rev 39)
+++ branches/robKalman_bs/pkg/robKalman/R/ACMfilter.R 2009-09-03 09:35:13 UTC (rev 40)
@@ -1,31 +1,24 @@
-.getcorrCovACM <- function (S1, K, Z, W=diag(nrow(Z)))
+.getcorrCovACM <- function (S1, K, Z, W=1)
{
###########################################
##
## R-function: .corrCov - computes filtering error covarince matrix
## (internal function)
## author: Bernhard Spangl
-## version: 1.0 (2006-05-22)
+## version: 1.1 (2009-07-29)
##
###########################################
## Paramters:
## S1 ... prediction error covariance matrix
## K ... Kalman gain
-## W ... weight matrix
+## W ... weight (scalar)
## Z ... observation matrix
- S1 - K %*% W %*% Z %*% S1
+ S1 - (K * W) %*% Z %*% S1
}
-##steps for classical Kalman filter (cK)
-.ACMinitstep <- function(a, S, ...)
- {dots <- list(...)
- if(hasArg("s0"))
- s0<-dots$"s0"
- else
- s0<-NULL
- list( x0 = a, S0 = S, s0 = s0)}
+## better: .mACMinitstep <- .cKinitstep
.ACMpredstep <- function (x0, S0, F, Q, i, rob0, s0, ...) ### S=P F= Phi
{
@@ -33,59 +26,83 @@
##
## R-function: .ACMpredstep - prediction step (internal function)
## author: Bernhard Spangl
-## version: 1.0 (2006-05-22)
+## version: 1.1 (2009-07-29)
##
###########################################
## Paramters:
-## x0 ... state vector (filter estimate)
-## F=Phi ... design matrix of state equation
-## S0 ... filtering error covariance matrix
+## x0 ... state vector (filter estimate), pd x runs matrix
+## S0 ... filtering error covariance matrix (formerly: P),
+## pd x pd x runs array (!!!)
+## F ... design matrix of state equation (formerly: Phi)
## Q ... covariance matrix of state innovation process
-## rob0 ... general robust parameter --- here: scale s0 of nominal Gaussain component of additive noise
- S1 <- .getpredCov(S0, F, Q)
- return(list(x1 = F %*% x0, S1 = S1, rob1 = sqrt(S1[1, 1] + s0), Ind=1))
+## i ... time index, currently not used
+## rob0 ... list, general robust parameter
+## here: sqrt of variance 'Rt',
+## i.e., of Var(Delta y_t)
+
+ S1 <- array(apply(S0, 3, .getpredCov, F, Q), dim=dim(S0))
+ rob1 <- NULL
+ return(list(x1 = F %*% x0, S1 = S1, rob1 = rob1, Ind = 0))
}
-.ACMcorrstep <- function (y, x1, S1, Z, V, i, rob1, dum=NULL, psi, apsi, bpsi, cpsi, flag, ...)
+.ACMcorrstep <- function (y, x1, S1, Z, V, i, rob1=NULL,
+ psi, apsi, bpsi, cpsi, flag, ...)
{
###########################################
##
## R-function: .ACMcorrstep - correction step (internal function)
## author: Bernhard Spangl
-## version: 1.0 (2006-05-22)
+## version: 1.1 (2009-07-29)
##
###########################################
## Paramters:
-## y ... univariate time series
-## x1 ... state vector (one-step-ahead predictor)
-## rob1 ... general robust parameter --- here st ... time-dependent scale parameter
-## S1 ... prediction error covariance matrix
-## Z ... observation matrix
-## dum ... dummy variable for compatibility with ... argument of calling function
-## V ... covariance matrix of observation noise
+## y ... observed univariate time series, qd(=1) x runs matrix
+## x1 ... state vector (one-step-ahead predictor), pd x runs matix
+## S1 ... prediction error covariance matrix (formerly: M),
+## pd x pd x runs array (!!!)
+## Z ... observation matrix (formerly: H), qd(=1) x pd
+## V ... covariance matrix of observation noise (formerly: R), qd(=1) x qd(=1)
+## i ... time index, currently not used
## psi ... influence function to be used
## a, b, c ... tuning constants for Hampel's psi-function
-## (defaul: a=b=2.5, c=5.0)
+## (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)
- st <- rob1
- K <- .getKG(S1, Z, V)
+## Warning: Algorithm can not be completely vectorized!
- rst <- (y - x1[1])/st
+ runs <- dim(S1)[3]
+ foo <- function (x, m) x %*% t(m)
+ S0 <- array(0, dim=dim(S1))
- ps <- psi(rst, apsi, bpsi, cpsi)[1,1]
- dx <- K * st * ps
+ D <- apply(S1, 3, foo, m=Z) # \
+ Rt <- drop(Z %*% D) + V # > st^2 = Z %*% S1 %*% t(Z) + V
+ st <- sqrt(Rt) # /
+
+ K <- t(t(D)/Rt) # Kalman Gain
+
+ rst <- drop(y - Z %*% x1)/st
+
+ ps <- psi(rst, apsi, bpsi, cpsi)
+ dx <- t(t(K) * (st * ps))
x0 <- x1 + dx
- ind <- (abs(rst-ps)>10^-8)
+ ind <- abs(rst) > apsi
- w <- psi(rst, apsi, bpsi, cpsi, flag)
-
- S0 <- .getcorrCovACM(S1, K, Z, W = w*diag(rep(1, nrow(Z))))
- Delta <- Z %*% S0 %*% t(Z) + V
+ w <- psi(rst, apsi, bpsi, cpsi, flag)
- return(list(x0 = x0, K = K, S0 = S0, Delta=Delta, Ind=ind, rob0=rob1, DeltaY = rst))
+ for (j in 1:runs) {
+ S0[, , j] <- .getcorrCovACM(S1=S1[, , j], K=K[, j], Z=Z, W=w[j])
+ }
+
+ return(list(x0 = x0,
+ K = aperm(array(K, dim=c(dim(K), 1)), c(1, 3, 2)),
+ S0 = S0,
+ Ind=ind,
+ rob0=st,
+ Delta=array(Rt, dim=c(1, 1, runs)),
+ DeltaY = matrix(rst*st, 1, runs)))
}
+
Modified: branches/robKalman_bs/pkg/robKalman/R/Psi.R
===================================================================
--- branches/robKalman_bs/pkg/robKalman/R/Psi.R 2009-07-23 19:51:16 UTC (rev 39)
+++ branches/robKalman_bs/pkg/robKalman/R/Psi.R 2009-09-03 09:35:13 UTC (rev 40)
@@ -128,25 +128,27 @@
Hampel = get("psiHampel", mode="function"))
}
-mvpsiHampel <- function (x, a=2, b=4, c=8)
+mvpsiHampel <- function (x, a=2, b=4, c=8, Norm)
{
###########################################
##
## R-function: mvpsiHampel - multivariate analogue of
## Hampel's psi-function
## author: Bernhard Spangl
-## version: 0.1 (2008-02-23)
+## version: 0.3 (2009-07-24)
##
###########################################
## Paramters:
-## x ... vector
+## x ... qd x runs matrix
## a, b, c ... tuning constants
- x.norm <- EuclideanNorm(x)
+## Norm ... which norm should be used? (default: EuclideanNorm)
+ x.norm <- Norm(x)
small <- (x.norm <= a)
dummy <- pmin(a, a/(c-b)*(c-x.norm))
dummy <- pmax(dummy, 0)/(x.norm+small)*(!small) + small
- return(x*dummy)
+ ## return(x*dummy)
+ return(t(t(x)*dummy))
}
jacobian.Hampel <- function (x, ...)
Modified: branches/robKalman_bs/pkg/robKalman/R/mACMfilter.R
===================================================================
--- branches/robKalman_bs/pkg/robKalman/R/mACMfilter.R 2009-07-23 19:51:16 UTC (rev 39)
+++ branches/robKalman_bs/pkg/robKalman/R/mACMfilter.R 2009-09-03 09:35:13 UTC (rev 40)
@@ -17,23 +17,60 @@
S1 - K %*% W %*% t(K)
}
-.mACMcorrstep <- function (y, x1, S1, Z, V, rob1=NULL,
- psi, apsi, bpsi, cpsi, flag, ...)
+## better: .mACMinitstep <- .cKinitstep
+
+.mACMpredstep <- function (x0, S0, F, Q, i, ..., rob0)
{
###########################################
##
-## R-function: .mACMcorrstep - correction step (internal function)
+## R-function: .mACMpredstep - prediction step (internal function)
## author: Bernhard Spangl
-## version: 0.2 (2008-03-31)
+## version: 0.3 (2009-07-27)
##
###########################################
## Paramters:
+## x0 ... state vector (filter estimate), pd x runs matrix
+## S0 ... filtering error covariance matrix (formerly: P),
+## pd x pd x runs array (!!!)
+## F ... design matrix of state equation (formerly: Phi)
+## Q ... covariance matrix of state innovation process
+## i ... time index, currently not used
+## rob0 ... list, general robust parameter
+## here: inverse sqrt of covariance matrix 'Rt',
+## i.e., of Cov(Delta y_t)
+
+## Warning: Algorithm can not be completely vectorized!
+
+ runs <- dim(x0)[2]
+ S1 <- array(0, dim=dim(S0))
+
+ for (j in 1:runs) {
+ S1[, , j] <- .getpredCov(S0[, , j], F, Q)
+ }
+ rob1 <- NULL
+ return(list(x1 = F %*% x0, S1 = S1, rob1 = rob1, Ind=0))
+}
+
+.mACMcorrstepInt <- function (y, x1, S1, Z, V, rob1=NULL, norm,
+ psi, apsi, bpsi, cpsi, flag, ...)
+{
+###########################################
+##
+## R-function: .mACMcorrstepInt - inner correction step
+## (internal function)
+## author: Bernhard Spangl
+## version: 0.3 (2009-07-24)
+##
+###########################################
+
+## Paramters:
## y ... observed vector-valued time series
## x1 ... state vector (one-step-ahead predictor)
## S1 ... prediction error covariance matrix (formerly: M)
## Z ... observation matrix (formerly: H)
## V ... covariance matrix of observation noise (formerly: R)
+## norm ... which norm should be used? (default: EuclideanNorm)
## psi ... influence function to be used
## a, b, c ... tuning constants for Hampel's psi-function
## (default: a=b=2.5, c=5.0)
@@ -47,21 +84,85 @@
sqrtR <- rootMatrix(Rt)
st <- sqrtR$X.sqrt.inv
- K <- D %*% st
+ K <- D %*% st # not exactly the classical Kalman Gain,
+ # 'K %*% st' wolud be
- yDelta <- drop(st %*% (y - Z %*% x1))
+ DeltaY <- y - Z %*% x1
+ yDelta <- st %*% DeltaY
+ ## yDelta <- drop(st %*% (y - Z %*% x1))
- yPsi <- psi(yDelta, apsi, bpsi, cpsi)
+ yPsi <- psi(yDelta, apsi, bpsi, cpsi, norm)
xDelta <- K %*% yPsi
x0 <- x1 + xDelta
- ind <- sqrt(yDelta%*%yDelta) > apsi
+ ind <- norm(yDelta) > apsi
jacobian.psi <- .weighting(flag)
- W <- jacobian.psi(yDelta, a=apsi, b=bpsi, c=cpsi)
+ W <- jacobian.psi(drop(yDelta), a=apsi, b=bpsi, c=cpsi, norm)
S0 <- .getcorrCovmACM(S1, K, W = W)
- return(list(x0 = x0, K = K, S0 = S0, Ind=ind, rob0=st, Delta =NULL, DeltaY = yDelta))
+ return(list(x0 = x0, K = K, S0 = S0, Ind = ind, rob0 = st, Delta = Rt,
+ DeltaY = DeltaY))
}
+.mACMcorrstep <- function (y, x1, S1, Z, V, i, rob1=NULL, norm,
+ psi, apsi, bpsi, cpsi, flag, ...)
+{
+###########################################
+##
+## R-function: .mACMcorrstep - correction step (internal function)
+## author: Bernhard Spangl
+## version: 0.3 (2009-07-24)
+##
+###########################################
+
+## Paramters:
+## y ... observed vector-valued time series, qd x runs matrix
+## x1 ... state vector (one-step-ahead predictor), pd x runs matix
+## S1 ... prediction error covariance matrix (formerly: M),
+## pd x pd x runs array (!!!)
+## Z ... observation matrix (formerly: H)
+## V ... covariance matrix of observation noise (formerly: R)
+## i ... time index, currently not used
+## norm ... which norm should be used? (default: EuclideanNorm)
+## psi ... influence function to be used
+## a, b, c ... 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)
+
+## Warning: Algorithm can not be completely vectorized!
+
+ pd <- dim(x1)[1]
+ qd <- dim(y)[1]
+ runs <- dim(S1)[3]
+
+ x0 <- matrix(0, pd, runs)
+ K <- array(0, dim = c(pd, qd, runs))
+ S0 <- array(0, dim = c(pd, pd, runs))
+ Ind <- NULL
+ rob0 <- array(0, dim = c(qd, qd, runs))
+ Delta <- array(0, dim = c(qd, qd, runs))
+ DeltaY <- matrix(0, qd, runs)
+
+ for (j in 1:runs) {
+ mACMcs <- .mACMcorrstepInt(y[, j], x1[, j], S1[, , j], Z, V,
+ rob1, norm, psi, apsi, bpsi, cpsi, flag, ...)
+
+ x0[, j] <- mACMcs$x0
+ K[, , j] <- mACMcs$K
+ S0[, , j] <- mACMcs$S0
+ Ind[j] <- mACMcs$Ind
+ rob0[, , j] <- mACMcs$rob0
+ Delta[, , j] <- mACMcs$Delta
+ DeltaY[, j] <- mACMcs$DeltaY
+ }
+
+ return(list(x0 = x0, K = K, S0 = S0, Ind=Ind, rob0 = rob0,
+ Delta = Delta, DeltaY = DeltaY))
+
+}
+
Modified: branches/robKalman_bs/pkg/robKalman/R/recFilter.R
===================================================================
--- branches/robKalman_bs/pkg/robKalman/R/recFilter.R 2009-07-23 19:51:16 UTC (rev 39)
+++ branches/robKalman_bs/pkg/robKalman/R/recFilter.R 2009-09-03 09:35:13 UTC (rev 40)
@@ -338,77 +338,83 @@
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)}
+ACMfilter <- function (Y, a, S, F, Q, Z, V,
+ nsim=0, seed=NULL,
+ psi=psiHampel, apsi=2.5, bpsi=2.5, cpsi=5.0,
+ flag="weights", dropRuns = TRUE,
+ saveOpt=TRUE, dimsCheck=NULL)
+{
###########################################
##
## R-function: ACMfilter - approximate conditional-mean filtering
## author: Bernhard Spangl
-## version: 1.0 (2006-05-21)
+## version: 1.1 (2009-07-29)
##
###########################################
## 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
+## Y ... observed univariate time series
+## (array with dimensions qd(=1!) x runs x tt)
+## a ... unconditional mean vector (formerly: m0)
+## S ... covariance matrix (formerly: Cx)
+## F ... design matrix of state equation (formerly: Phi), array, pd x pd x tt
+## Q ... covariance matrix of state innovation process, array, pd x pd x tt
+## Z ... observation matrix (formerly: H), array, qd x pd x tt
+## V ... covariance matrix of observation noise (formerly: R),
+## array, qd x qd x tt
+## nsim ... if >0 we simulate a bunch of nsim paths (acc. to ideal model) to
+## get emp. covariances
+## seed ... seed for the simulations
## 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)
+## 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)
+## dropRuns ... shall run-dimension be collapsed if it is one?
+## saveOpt ... logical (default: TRUE), should the stuff really be saved?
+## dimsCheck ... either 'NULL' (default) or vector [pd, qd, runs, tt] with
+## correct dimensions
-mACMfilter <- function(Y, a, S, F, Q, Z, V,
- psi=mvpsiHampel, apsi=2.5, bpsi=2.5, cpsi=5.0,
- flag="deriv", dropRuns = TRUE)
+ warning("The algorithm is not completely vectorized!
+ Please be patient... \n")
+
+ recursiveFilter(Y, a, S, F, Q, Z, V,
+ initSc=.cKinitstep, predSc=.cKpredstep, corrSc=.cKcorrstep,
+ initSr=.cKinitstep, predSr=.ACMpredstep, corrSr=.ACMcorrstep,
+ nsim=nsim, seed=seed,
+ psi=psi, apsi=apsi, bpsi=bpsi, cpsi=cpsi, flag=flag, dropRuns=dropRuns,
+ CovRunDep=TRUE, saveOpt=saveOpt, dimsCheck=dimsCheck)
+}
+
+mACMfilter <- function (Y, a, S, F, Q, Z, V,
+ nsim=0, seed=NULL, norm=EuclideanNorm,
+ psi=mvpsiHampel, apsi=2.5, bpsi=2.5, cpsi=5.0,
+ flag="deriv", dropRuns = TRUE,
+ saveOpt=TRUE, dimsCheck=NULL)
{
###########################################
##
## R-function: mACMfilter - approximate conditional-mean filtering
## author: Bernhard Spangl
-## version: 0.2 (2008-03-31)
+## version: 0.3 (2009-07-28)
##
###########################################
## Paramters:
## Y ... observed vector-valued time series
-## (column-wise, matrix: q rows, number of columns equal to time points)
+## (array with dimensions qd x runs x tt)
## 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)
+## F ... design matrix of state equation (formerly: Phi), array, pd x pd x tt
+## Q ... covariance matrix of state innovation process, array, pd x pd x tt
+## Z ... observation matrix (formerly: H), array, qd x pd x tt
+## V ... covariance matrix of observation noise (formerly: R),
+## array, qd x qd x tt
+## nsim ... if >0 we simulate a bunch of nsim paths (acc. to ideal model) to
+## get emp. covariances
+## seed ... seed for the simulations
+## norm ... which norm should be used? (default: EuclideanNorm)
## 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
@@ -417,10 +423,19 @@
## if "deriv" (default), Jacobian matrix of multivariate analogue
## of Hampel's psi-function is used (only default is available
## at the moment)
+## dropRuns ... shall run-dimension be collapsed if it is one?
+## saveOpt ... logical (default: TRUE), should the stuff really be saved?
+## dimsCheck ... either 'NULL' (default) or vector [pd, qd, runs, tt] with
+## correct dimensions
+ warning("The algorithm is not completely vectorized!
+ Please be patient... \n")
+
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)
+ initSc=.cKinitstep, predSc=.cKpredstep, corrSc=.cKcorrstep,
+ initSr=.cKinitstep, predSr=.mACMpredstep, corrSr=.mACMcorrstep,
+ nsim=nsim, seed=seed, norm=norm,
+ psi=psi, apsi=apsi, bpsi=bpsi, cpsi=cpsi, flag=flag, dropRuns=dropRuns,
+ CovRunDep=TRUE, saveOpt=saveOpt, dimsCheck=dimsCheck)
}
More information about the Robkalman-commits
mailing list