[Robkalman-commits] r32 - branches/robkalman_pr/pkg/robKalman/R branches/robkalman_pr/pkg/robKalman/demo branches/robkalman_pr/pkg/robKalman/man pkg/robKalman/R pkg/robKalman/man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Jun 12 19:21:08 CEST 2009
Author: ruckdeschel
Date: 2009-06-12 19:21:07 +0200 (Fri, 12 Jun 2009)
New Revision: 32
Modified:
branches/robkalman_pr/pkg/robKalman/R/ACMfilter.R
branches/robkalman_pr/pkg/robKalman/R/classKalman.R
branches/robkalman_pr/pkg/robKalman/R/mACMfilter.R
branches/robkalman_pr/pkg/robKalman/R/rLSfilter.R
branches/robkalman_pr/pkg/robKalman/R/recFilter.R
branches/robkalman_pr/pkg/robKalman/demo/rLSdemo.R
branches/robkalman_pr/pkg/robKalman/man/internalKalman.Rd
branches/robkalman_pr/pkg/robKalman/man/internalrLS.Rd
branches/robkalman_pr/pkg/robKalman/man/recFilter.Rd
pkg/robKalman/R/ACMfilter.R
pkg/robKalman/R/classKalman.R
pkg/robKalman/R/mACMfilter.R
pkg/robKalman/R/rLSfilter.R
pkg/robKalman/R/recFilter.R
pkg/robKalman/R/simulateSScont.R
pkg/robKalman/man/internalKalman.Rd
pkg/robKalman/man/internalrLS.Rd
pkg/robKalman/man/recFilter.Rd
Log:
merged time variant SSM formulation back into trunc;
modified return values (e.g. Delta y's and their cov's are returned now)
Modified: branches/robkalman_pr/pkg/robKalman/R/ACMfilter.R
===================================================================
--- branches/robkalman_pr/pkg/robKalman/R/ACMfilter.R 2009-06-10 10:45:12 UTC (rev 31)
+++ branches/robkalman_pr/pkg/robKalman/R/ACMfilter.R 2009-06-12 17:21:07 UTC (rev 32)
@@ -93,5 +93,5 @@
S0 <- .getcorrCovACM(S1, K, Z, W = w*diag(rep(1, nrow(Z))))
- return(list(x0 = x0, K = K, S0 = S0, Ind = ind, rob0 = rob1))
+ return(list(x0 = x0, K = K, S0 = S0, Delta=NULL, Ind=ind, rob0=rob1, DeltaY = rst))
}
Modified: branches/robkalman_pr/pkg/robKalman/R/classKalman.R
===================================================================
--- branches/robkalman_pr/pkg/robKalman/R/classKalman.R 2009-06-10 10:45:12 UTC (rev 31)
+++ branches/robkalman_pr/pkg/robKalman/R/classKalman.R 2009-06-12 17:21:07 UTC (rev 32)
@@ -36,10 +36,14 @@
#Kalman filter routines
#
#
-.getKG <- function(S1, Z, V)# calculates the Kalman Gain for S1 = S_{t|t-1}, Z, V as above
+.getDelta <- function(S1, Z, V)# calculates the Cov(Delta y_t) for S1 = S_{t|t-1}, Z, V as above
{ H <- S1 %*% t(Z)
- H %*% ginv( Z%*%H + V ) }
+ ginv( Z%*%H + V ) }
+.getKG <- function(S1, Z, Delta)# calculates the Kalman Gain for S1 = S_{t|t-1}, Z, V as above
+{ H <- S1 %*% t(Z)
+ H %*% Delta }
+
.getcorrCov <- function(S1, K, Z)# calculates S_{t|t} for S1 = S_{t|t-1}, K_t, Z as above
{ S1 - K %*% Z %*% S1 }
@@ -47,14 +51,16 @@
{ F %*% S0 %*% t(F) + Q }
##steps for classical Kalman filter (cK)
-.cKinitstep <- function(a, S, ...)
+.cKinitstep <- function(a, S, ...)
{list( x0 = a, S0 = S )}
-
-.cKpredstep <- function(x0, S0, F, Q, ...)
- {list( x1 = F %*% x0, S1 = .getpredCov(S0, F, Q), Ind = FALSE)}
-.cKcorrstep <- function(y, x1, S1, Z, V, ...)
- {K <- .getKG(S1, Z, V)
- x0 <- x1 + K %*% (y - Z %*% x1)
+.cKpredstep <- function(x0, S0, F, Q, ...)
+ {list( x1 = F %*% x0, S1 = .getpredCov(S0, F, Q), Ind=1)}
+
+.cKcorrstep <- function(y, x1, S1, Z, V, ...)
+ {Delta <- .getDelta(S1, Z, V)
+ K <- .getKG(S1, Z, Delta)
+ DeltaY <- y - Z %*% x1
+ x0 <- x1 + K %*% DeltaY
S0 <- .getcorrCov(S1, K, Z)
- list(x0 = x0, K = K, S0 = S0, Ind = FALSE)}
+ list(x0 = x0, K = K, S0 = S0, Delta = Delta, Ind=1, DeltaY = yDelta)}
Modified: branches/robkalman_pr/pkg/robKalman/R/mACMfilter.R
===================================================================
--- branches/robkalman_pr/pkg/robKalman/R/mACMfilter.R 2009-06-10 10:45:12 UTC (rev 31)
+++ branches/robkalman_pr/pkg/robKalman/R/mACMfilter.R 2009-06-12 17:21:07 UTC (rev 32)
@@ -62,6 +62,6 @@
S0 <- .getcorrCovmACM(S1, K, W = W)
- return(list(x0 = x0, K = K, S0 = S0, Ind=ind, rob0=st))
+ return(list(x0 = x0, K = K, S0 = S0, Ind=ind, rob0=st, Delta =NULL, DeltaY = yDelta))
}
Modified: branches/robkalman_pr/pkg/robKalman/R/rLSfilter.R
===================================================================
--- branches/robkalman_pr/pkg/robKalman/R/rLSfilter.R 2009-06-10 10:45:12 UTC (rev 31)
+++ branches/robkalman_pr/pkg/robKalman/R/rLSfilter.R 2009-06-12 17:21:07 UTC (rev 32)
@@ -7,14 +7,14 @@
## we assume a time-invariant SSM of the following form
# Hyper-parameters:
-# F (p x p), 0<= S (p x p), 0< Q (p x p)
-# Z (q x p), 0< V (q x q)
-# distributional assumptions
+# F (p x p), 0<= S (p x p), 0< Q (p x p)
+# Z (q x p), 0< V (q x q)
+# distributional assumptions
# +initial condition x_0 ~ N_p(a,S)
# +innovations v_t ~ N_p(0,Q), t>=1
# +observation errors e_t ~ N_q(0,V), t>=1
-#
-# Model equations
+#
+# Model equations
# +state equation x_t = F x_{t-1} + v_t, t>=1
# +observation equation y_t = Z x_t + e_t , t>=1
##
@@ -23,12 +23,12 @@
### Notation for the Kalman filter
-# +initial step x_{0|0} = a
+# +initial step x_{0|0} = a
# error covariance S_{0|0} = Cov(x_0-x_{0|0}) = S
# +prediction step x_{t|t-1} = F x_{t-1|t-1}, t>=1
-# error covariance S_{t|t-1} = Cov(x_t-x_{t|t-1}) = F S_{t-1|t-1} F' + Q
+# error covariance S_{t|t-1} = Cov(x_t-x_{t|t-1}) = F S_{t-1|t-1} F' + Q
# +correction step x_{t|t} = x_{t|t-1} + K_t (y_t - Z x_{t|t-1}) t>=1
-# for Kalman Gain K_t = S_{t|t-1} Z' (Z S_{t|t-1} Z' + V )^-
+# for Kalman Gain K_t = S_{t|t-1} Z' (Z S_{t|t-1} Z' + V )^-
# error covariance S_{t|t} = Cov(x_t-x_{t|t}) = S_{t|t-1} - K_t Z S_{t|t-1}
#########################################################################################
@@ -39,16 +39,18 @@
##steps for rLS filter (rLS) (xr0, xr1)
## (simultaneously for cK (x0, x1) )
-.rLScorrstep <- function(y, x1, S1, Z, V, i, rob1 = NULL, b,
- norm = EuclideanNorm, ...)
- {K <- .getKG(S1, Z, V)
- dx <- K %*% (y - Z %*% x1)
- if(length(b)>1)
+.rLScorrstep <- function(y, x1, S1, Z, V, i, rob1 = NULL, b,
+ norm = EuclideanNorm, ...)
+ {Delta <- .getDelta(S1, Z, V)
+ K <- .getKG(S1, Z, Delta)
+ DeltaY <- y - Z %*% x1
+ dx <- K %*% DeltaY
+ if(length(b)>1)
b <- b[min(i,length(b))]
x0 <- x1 + Huberize(dx, b, norm = norm)
S0 <- .getcorrCov(S1, K, Z)
- if (ncol(x1)==1)
+ if (ncol(x1)==1)
Ind <- (norm(dx)>1)
- else
+ else
Ind <- apply(dx, 2, function(xx) norm(xx)>1)
- list(x0 = x0, K = K, S0 = S0, Ind = Ind)}
+ list(x0 = x0, K = K, S0 = S0, Delta = Delta, Ind = Ind, DeltaY = DeltaY)}
Modified: branches/robkalman_pr/pkg/robKalman/R/recFilter.R
===================================================================
--- branches/robkalman_pr/pkg/robKalman/R/recFilter.R 2009-06-10 10:45:12 UTC (rev 31)
+++ branches/robkalman_pr/pkg/robKalman/R/recFilter.R 2009-06-12 17:21:07 UTC (rev 32)
@@ -41,16 +41,16 @@
-recursiveFilter <- function(Y, a, S, F, Q, Z, V,
- initSc = .cKinitstep, predSc = .cKpredstep,
- corrSc = .cKcorrstep,
- initSr = NULL, predSr = NULL, corrSr = NULL,
+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)
# a generalization of the Kalmanfilter
-#arguments:
+#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
+# :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
@@ -73,12 +73,12 @@
########################
tt <- (dim(Y))[3]
- runs <- (dim(Y))[2]
+ runs <- (dim(Y))[2]
pd <- length(a)
IndIO <- NULL
IndAO <- NULL
- St0s <- St1s <- NULL
+ St0s <- St1s <- NULL
if(is.matrix(F)) F <- array(F, dim = c(pd,pd,tt))
if(is.matrix(Z)) Z <- array(Z, dim = c(pd,qd,tt))
@@ -86,17 +86,19 @@
if(is.matrix(V)) V <- array(V, dim = c(qd,qd,tt))
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))
+ DeltaY <- array(0, dim = c(qd, runs, tt))
St0 <- array(0, dim = c(pd, pd, tt + 1))
St1 <- array(0, dim = c(pd, pd, tt))
KG <- array(0, dim = c(pd, qd, tt))
-
+ Delta <- array(0, dim = c(qd, qd, tt))
+
if (nsim)
- {if (!exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE))
+ {if (!exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE))
runif(1)
- if (is.null(seed))
+ if (is.null(seed))
RNGstate <- get(".Random.seed", envir = .GlobalEnv)
else {
R.seed <- get(".Random.seed", envir = .GlobalEnv)
@@ -114,6 +116,8 @@
Str0 <- array(0, dim = c(pd, pd, tt + 1))
Str1 <- array(0, dim = c(pd, pd, tt))
KGr <- array(0, dim = c(pd, qd, tt))
+ Deltar <- array(0, dim = c(qd, qd, tt))
+ DeltaYr <- array(0, dim = c(qd, runs, tt))
}
if(!is.null(predSr))
@@ -121,43 +125,44 @@
if(!is.null(corrSr))
IndAO <- matrix(FALSE,runs,tt)
-
+
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(rmvnorm(nsim, a, S))
St0s <- array(0, c(pd, pd, tt))
- St1s <- 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[,, 1] <- xr0
- rob0L <- list(rob0)
-
+ rob0L <- list(rob0)
+
}
else{Xrf <- NULL
Xrp <- NULL
Str0 <- NULL
Str1 <- NULL
KGr <- NULL
+ Deltar <- NULL
rob0L <- NULL
rob1L <- NULL
}
@@ -167,15 +172,15 @@
Q0 <- matrix(Q[,,i], nrow = pd, ncol = pd)
ps <- predSc(x0 = x0, S0 = S0, F = F0, Q = Q0, i = i, ...)
- x1 <- ps$x1
+ 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,
+ {psr <- predSr(x0 = xr0, S0 = Sr0, F = F0, Q = Q0, i = i,
..., rob0 = rob0)
IndIO[,i] <- as.logical(psr$Ind)
if(nsim){
@@ -183,22 +188,22 @@
Xs <- F0 %*% Xs + vs
xr1s <- predSr(x0 = xr0s, S0 = Sr0, F = F0, Q = Q0, i = i,
..., rob0 = rob0)$x1
- St1s[,,i] <- cov(t(xr1s))
+ St1s[,,i] <- cov(t(xr1s))
}
- }else{
+ }else{
psr <- predSc(x0 = xr0, S0 = Sr0, F = F0, Q = Q0, i = i, ...)
if(nsim){
vs <- t(rmvnorm(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))
+ St1s[,,i] <- cov(t(xr1s))
}
}
xr1 <- psr$x1
Sr1 <- psr$S1
- rob1 <- psr$rob1
-
+ rob1 <- psr$rob1
+
Xrp[,, i] <- xr1
Str1[,, i]<- S1
if(i==1) rob1L <- list(rob1)
@@ -214,42 +219,46 @@
cs <- corrSc(y = Y0, x1 = x1, S1 = S1, Z = Z0, V = V0, i = i,...)
x0 <- cs$x0
S0 <- cs$S0
-
+ DeltaY[,,i] <- cs$DeltaY
+
Xf[,, i + 1] <- x0
St0[,, i + 1] <- S0
KG[,, i] <- cs$K
-
+ Delta[,, i] <- cs$Delta
+
if(robust)
{if (!is.null(corrSr))
- {csr <- corrSr(y = Y0, x1 = xr1, S1 = Sr1,
+ {csr <- corrSr(y = Y0, x1 = xr1, S1 = Sr1,
Z = Z0, V = V0, i = i, ..., rob1 = rob1)
- IndAO[,i] <- as.logical(csr$Ind)
+ IndAO[,i] <- as.logical(csr$Ind)
if(nsim){
es <- t(rmvnorm(nsim, Y[,1]*0, V0))
Ys <- Z0 %*% Xs + es
- xr0s <- corrSr(y = Ys, x1 = xr1s, S1 = Sr1,
+ xr0s <- corrSr(y = Ys, x1 = xr1s, S1 = Sr1,
Z = Z0, V = V0, i = i, ..., rob1 = rob1)$x0
- St0s[,,i] <- cov(t(xr0s))
+ St0s[,,i] <- cov(t(xr0s))
}
}else{
- csr <- corrSc(y = Y0, x1 = xr1, S1 = Sr1, Z = Z0, V = V0,
+ csr <- corrSc(y = Y0, x1 = xr1, S1 = Sr1, Z = Z0, V = V0,
i = i, ...)
if(nsim){
es <- t(rmvnorm(nsim, Y[,1]*0, V0))
Ys <- Z0 %*% Xs + es
- xr0s <- corrSc(y = Ys, x1 = xr1s, S1 = Sr1,
+ xr0s <- corrSc(y = Ys, x1 = xr1s, S1 = Sr1,
Z = Z0, V = V0, i = i, ...)$x0
- St0s[,,i] <- cov(t(xr0s))
+ St0s[,,i] <- cov(t(xr0s))
}
- }
+ }
xr0 <- csr$x0
Sr0 <- csr$S0
- rob0 <- csr$rob0
+ rob0 <- csr$rob0
+ DeltaYr[,,i] <- cs$DeltaY
- Xrf[,, i + 1] <- xr0
+ Xrf[,, i + 1] <- xr0
Str0[,, i + 1] <- S0
- rob0L[[i + 1]] <- rob0
+ rob0L[[i + 1]] <- rob0
KGr[,, i] <- csr$K
+ Deltar[,, i] <- cs$Deltar
}
}
@@ -258,14 +267,16 @@
Xp <- matrix(Xp,pd,tt)
if(!is.null(Xrp)) {
Xrf <- matrix(Xrf,pd,tt+1)
- Xrp <- matrix(Xrp,pd,tt)
+ Xrp <- matrix(Xrp,pd,tt)
IndIO <- as.logical(IndIO)
- IndAO <- as.logical(IndAO)
+ IndAO <- as.logical(IndAO)
}}
-list(Xf = Xf, Xp = Xp, Xrf = Xrf, Xrp = Xrp,
- S0 = St0, S1 = St1, KG = KG,
- Sr0 = Str0, Sr1 = Str1,
- KGr = KGr, IndIO = IndIO, IndAO = 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)
Modified: branches/robkalman_pr/pkg/robKalman/demo/rLSdemo.R
===================================================================
--- branches/robkalman_pr/pkg/robKalman/demo/rLSdemo.R 2009-06-10 10:45:12 UTC (rev 31)
+++ branches/robkalman_pr/pkg/robKalman/demo/rLSdemo.R 2009-06-12 17:21:07 UTC (rev 32)
@@ -62,9 +62,9 @@
mc = m0c, Vc = V0c, r = ract, rcalib=r1, effcalib=eff1)
{
#Simulation::
-X <- simulateState(a = a0, S = Ss, F = F, Q = Q, tt = tt)
-Yid <- simulateObs(X = X, Z = Z, Vi = Vi, mc = mc, Vc = Vc, r = 0)
-Yre <- simulateObs(X = X, Z = Z, Vi = Vi, mc = mc, Vc = Vc, r = ract)
+X <- simulateState(a = a0, S = Ss, F = F, Q = Q, tt = tt, runs = 1)
+Yid <- simulateObs(X = X, Z = Z, Vi = Vi, mc = mc, Vc = Vc, r = 0, runs = 1)
+Yre <- simulateObs(X = X, Z = Z, Vi = Vi, mc = mc, Vc = Vc, r = ract, runs = 1)
pd <- dim(X)[1]
qd <- dim(Yid)[1]
Modified: branches/robkalman_pr/pkg/robKalman/man/internalKalman.Rd
===================================================================
--- branches/robkalman_pr/pkg/robKalman/man/internalKalman.Rd 2009-06-10 10:45:12 UTC (rev 31)
+++ branches/robkalman_pr/pkg/robKalman/man/internalKalman.Rd 2009-06-12 17:21:07 UTC (rev 32)
@@ -12,7 +12,8 @@
These functions are used internally by package robKalman
}
\usage{
-.getKG(S1, Z, V)
+.getDelta(S1, Z, V)
+.getKG(S1, Z, Delta)
.getcorrCov(S1, K, Z)
.getpredCov(S0, F, Q)
.cKinitstep(a, S, ...)
@@ -27,6 +28,7 @@
\item{V}{observation error covariance (see below)}
\item{F}{innovation transition matrix (see below)}
\item{Q}{innovation covariance (see below)}
+ \item{Delta}{Covariance of \eqn{\Delta y_t}{Delta y_t}}
\item{K}{Kalman gain \eqn{K_t}}
\item{S1}{prediction error covariance \eqn{S_{t|t-1}} of the classical Kalman filter}
\item{S0}{filter error covariance \eqn{S_{t-1|t-1}} of the classical Kalman filter}
@@ -37,48 +39,73 @@
}
\details{
-We work in the setup of the time-invariant, linear, Gaussian state space model (ti-l-G-SSM)
+We work in the setup of the linear, Gaussian state space model (l-G-SSM)
with \eqn{p} dimensional states \eqn{x_t} and \eqn{q} dimensional observations \eqn{y_t},
with \strong{initial condition}
\deqn{x_0 \sim {\cal N}_p(a,S),}{x_0 ~ N_p(a,S),}
\strong{state equation}
-\deqn{x_t = F x_{t-1} + v_t, \qquad v_t \sim {\cal N}_p(0,Q),\qquad t\ge 1,}{x_t = F x_{t-1} + v_t, v_t ~ N_p(0,Q), t>=1,}
+\deqn{x_t = F_t x_{t-1} + v_t, \qquad v_t \sim {\cal N}_p(0,Q_t),\qquad t\ge 1,
+ }{x_t = F_t x_{t-1} + v_t, v_t ~ N_p(0,Q_t), t>=1,}
\strong{observation equation}
-\deqn{y_t = Z x_{t} + \varepsilon_t, \qquad \varepsilon_t \sim {\cal N}_q(0,V),\qquad t\ge 1,}{y_t = Z x_t + e_t, e_t ~ N_q(0,V), t>=1,}
+\deqn{y_t = Z_t x_{t} + \varepsilon_t, \qquad \varepsilon_t \sim
+{\cal N}_q(0,V_t),\qquad t\ge 1,}{y_t = Z_t x_t + e_t, e_t ~ N_q(0,V_t), t>=1,}
and where all random variable \eqn{x_0}, \eqn{v_t}, \eqn{\varepsilon_t}{e_t} are independent.
For notation, let us formulate the classical Kalman filter in this context:
-\strong{(0) ininitial step} \deqn{x_{0|0} = a}
-\eqn{\qquad}{\code{ }} with error covariance
+\strong{(0) ininitial step} \deqn{x_{0|0} = a}
+\eqn{\qquad}{\code{ }} with error covariance
\deqn{S_{0|0} = {\rm Cov}(x_0-x_{0|0}) = \code{S}}{S_{0|0} = Cov(x_0-x_{0|0}) = S}%
-\strong{(1) prediction step}
-\deqn{x_{t|t-1} = F x_{t-1|t-1},\qquad t\ge 1}{x_{t|t-1} = F x_{t-1|t-1}, t>=1}
-\eqn{\qquad}{\code{ }} with error covariance
-\deqn{S_{t|t-1} = {\rm Cov}(x_t-x_{t|t-1}) = F S_{t-1|t-1} F' + Q}{S_{t|t-1} = Cov(x_t-x_{t|t-1}) = F S_{t-1|t-1} F' + Q}
+\strong{(1) prediction step}
+\deqn{x_{t|t-1} = F_t x_{t-1|t-1},\qquad t\ge 1}{x_{t|t-1} = F_t x_{t-1|t-1}, t>=1}
+\eqn{\qquad}{\code{ }} with error covariance
+\deqn{S_{t|t-1} = {\rm Cov}(x_t-x_{t|t-1}) =
+F_t S_{t-1|t-1} F_t' + Q_t}{S_{t|t-1} = Cov(x_t-x_{t|t-1}) = F_t S_{t-1|t-1} F_t' + Q_t}
-\strong{(2) correction step}
-\deqn{x_{t|t} = x_{t|t-1} + K_t (y_t - Z x_{t|t-1}),\qquad t\ge 1}{x_{t|t} = x_{t|t-1} + K_t (y_t - Z x_{t|t-1}), t>=1}
-\eqn{\qquad}{\code{ }} for Kalman Gain
-\deqn{K_t = S_{t|t-1} Z' (Z S_{t|t-1} Z' + V )^-}
-\eqn{\qquad}{\code{ }} with error covariance
-\deqn{S_{t|t} = {\rm Cov}(x_t-x_{t|t}) = S_{t|t-1} - K_t Z S_{t|t-1}}{S_{t|t} = Cov(x_t-x_{t|t}) = S_{t|t-1} - K_t Z S_{t|t-1}}
+\strong{(2) correction step}
+\deqn{x_{t|t} = x_{t|t-1} + K_t (y_t - Z_t x_{t|t-1}),\qquad
+t\ge 1}{x_{t|t} = x_{t|t-1} + K_t (y_t - Z_t x_{t|t-1}), t>=1}
+\eqn{\qquad}{\code{ }} for Kalman Gain
+\deqn{K_t = S_{t|t-1} Z_t' (Z_t S_{t|t-1} Z_t' + V_t )^-}
+\eqn{\qquad}{\code{ }} with error covariance
+\deqn{S_{t|t} = {\rm Cov}(x_t-x_{t|t}) = S_{t|t-1} - K_t Z_t S_{t|t-1}}{S_{t|t} =
+Cov(x_t-x_{t|t}) = S_{t|t-1} - K_t Z_t S_{t|t-1}}
}
\value{
+\code{.getDelta} calculates the covariance of \eqn{\Delta y_t}{Delta y_t}
+ for \code{S1 =}\eqn{S_{t|t-1}}, \code{Z}, \code{V} as above \cr
+
\code{.getKG} calculates the Kalman Gain for \code{S1 =}\eqn{S_{t|t-1}}, \code{Z}, \code{V} as above \cr
+
\code{.getcorrCov} calculates \eqn{S_{t|t}} for \code{S1 =} \eqn{S_{t|t-1}}, \code{K =}\eqn{K_t} and \code{Z} as above \cr
+
\code{.getpredCov} calculates \eqn{S_{t|t-1}} for \code{S0 =} \eqn{S_{t-1|t-1}}, \code{F}, and \code{Q} as above\cr
+
\code{.cKinitstep} calculates \eqn{x_{0|0}} for \code{a}, \code{S} as above\cr
+The return value is a list with components
+\code{x0} (the filtered value)
+\code{S0} (the filter error covariance)\cr
+
\code{.cKpredstep} calculates \eqn{x_{t|t-1}} for \code{x0 =}\eqn{x_{t-1|t-1}}, \code{S0 =}\eqn{S_{t-1|t-1}}, and \code{F}, \code{Q}\cr
+The return value is a list with components
+\code{x1} (the predicted values),
+\code{S1} (the prediction error covariance),
+\code{Ind}(the indicators of clipped runs)\cr
\code{.cKcorrstep} calculates \eqn{x_{t|t}} for \code{x1 =}\eqn{x_{t|t-1}}, \code{y =}\eqn{y_{t}},
\code{S1 =}\eqn{S_{t|t-1}}, and \code{Z}, \code{V} \cr
+The return value is a list with components
+\code{x0} (the filtered values), \code{K} (the Kalman gain),
+\code{S0} (the filter error covariance),
+\code{Delta} (the covariance of \eqn{\Delta y_t}{Delta y_t}),
+\code{DeltaY} (the observation residuals \eqn{\Delta y_t}{Delta y_t}),
+\code{Ind}(the indicators of clipped runs)
}
Modified: branches/robkalman_pr/pkg/robKalman/man/internalrLS.Rd
===================================================================
--- branches/robkalman_pr/pkg/robKalman/man/internalrLS.Rd 2009-06-10 10:45:12 UTC (rev 31)
+++ branches/robkalman_pr/pkg/robKalman/man/internalrLS.Rd 2009-06-12 17:21:07 UTC (rev 32)
@@ -29,42 +29,52 @@
}
\details{
-We work in the setup of the time-invariant, linear, Gaussian state space model (ti-l-G-SSM)
+We work in the setup of the linear, Gaussian state space model (l-G-SSM)
with \eqn{p} dimensional states \eqn{x_t} and \eqn{q} dimensional observations \eqn{y_t},
with \strong{initial condition}
\deqn{x_0 \sim {\cal N}_p(a,S),}{x_0 ~ N_p(a,S),}
\strong{state equation}
-\deqn{x_t = F x_{t-1} + v_t, \qquad v_t \sim {\cal N}_p(0,Q),\qquad t\ge 1,}{x_t = F x_{t-1} + v_t, v_t ~ N_p(0,Q), t>=1,}
+\deqn{x_t = F_t x_{t-1} + v_t, \qquad v_t \sim {\cal N}_p(0,Q_t),\qquad t\ge 1,
+ }{x_t = F_t x_{t-1} + v_t, v_t ~ N_p(0,Q_t), t>=1,}
\strong{observation equation}
-\deqn{y_t = Z x_{t} + \varepsilon_t, \qquad \varepsilon_t \sim {\cal N}_q(0,V),\qquad t\ge 1,}{y_t = Z x_t + e_t, e_t ~ N_q(0,V), t>=1,}
+\deqn{y_t = Z_t x_{t} + \varepsilon_t, \qquad \varepsilon_t \sim
+{\cal N}_q(0,V_t),\qquad t\ge 1,}{y_t = Z_t x_t + e_t, e_t ~ N_q(0,V_t), t>=1,}
and where all random variable \eqn{x_0}, \eqn{v_t}, \eqn{\varepsilon_t}{e_t} are independent.
For notation, let us formulate the classical Kalman filter in this context:
-\strong{(0) ininitial step} \deqn{x_{0|0} = a}
-\eqn{\qquad}{\code{ }} with error covariance
+\strong{(0) ininitial step} \deqn{x_{0|0} = a}
+\eqn{\qquad}{\code{ }} with error covariance
\deqn{S_{0|0} = {\rm Cov}(x_0-x_{0|0}) = \code{S}}{S_{0|0} = Cov(x_0-x_{0|0}) = S}%
-\strong{(1) prediction step}
-\deqn{x_{t|t-1} = F x_{t-1|t-1},\qquad t\ge 1}{x_{t|t-1} = F x_{t-1|t-1}, t>=1}
-\eqn{\qquad}{\code{ }} with error covariance
-\deqn{S_{t|t-1} = {\rm Cov}(x_t-x_{t|t-1}) = F S_{t-1|t-1} F' + Q}{S_{t|t-1} = Cov(x_t-x_{t|t-1}) = F S_{t-1|t-1} F' + Q}
+\strong{(1) prediction step}
+\deqn{x_{t|t-1} = F_t x_{t-1|t-1},\qquad t\ge 1}{x_{t|t-1} = F_t x_{t-1|t-1}, t>=1}
+\eqn{\qquad}{\code{ }} with error covariance
+\deqn{S_{t|t-1} = {\rm Cov}(x_t-x_{t|t-1}) =
+F_t S_{t-1|t-1} F_t' + Q_t}{S_{t|t-1} = Cov(x_t-x_{t|t-1}) = F_t S_{t-1|t-1} F_t' + Q_t}
-\strong{(2) correction step}
-\deqn{x_{t|t} = x_{t|t-1} + K_t (y_t - Z x_{t|t-1}),\qquad t\ge 1}{x_{t|t} = x_{t|t-1} + K_t (y_t - Z x_{t|t-1}), t>=1}
-\eqn{\qquad}{\code{ }} for Kalman Gain
-\deqn{K_t = S_{t|t-1} Z' (Z S_{t|t-1} Z' + V )^-}
-\eqn{\qquad}{\code{ }} with error covariance
-\deqn{S_{t|t} = {\rm Cov}(x_t-x_{t|t}) = S_{t|t-1} - K_t Z S_{t|t-1}}{S_{t|t} = Cov(x_t-x_{t|t}) = S_{t|t-1} - K_t Z S_{t|t-1}}
-
+\strong{(2) correction step}
+\deqn{x_{t|t} = x_{t|t-1} + K_t (y_t - Z_t x_{t|t-1}),\qquad
+t\ge 1}{x_{t|t} = x_{t|t-1} + K_t (y_t - Z_t x_{t|t-1}), t>=1}
+\eqn{\qquad}{\code{ }} for Kalman Gain
+\deqn{K_t = S_{t|t-1} Z_t' (Z_t S_{t|t-1} Z_t' + V_t )^-}
+\eqn{\qquad}{\code{ }} with error covariance
+\deqn{S_{t|t} = {\rm Cov}(x_t-x_{t|t}) = S_{t|t-1} - K_t Z_t S_{t|t-1}}{S_{t|t} =
+Cov(x_t-x_{t|t}) = S_{t|t-1} - K_t Z_t S_{t|t-1}}
}
\value{
\code{.rLScorrstep(y, x1, S1, Z, V, b)} calculates rLS-\eqn{x_{t|t}} for arguments \code{b}, \code{x1 =}\eqn{x_{t|t-1}} (from rLS-past),
\code{y =}\eqn{y_{t}}, \code{S1 =}\eqn{S_{t|t-1}}, \code{Z}, and \code{V} as above\cr
+The return value is a list with components
+\code{x0} (the filtered values), \code{K} (the Kalman gain),
+\code{S0} (the filter error covariance),
+\code{Delta} (the covariance of \eqn{\Delta y_t}{Delta y_t}),
+\code{DeltaY} (the observation residuals \eqn{\Delta y_t}{Delta y_t}),
+\code{Ind}(the indicators of clipped runs)
}
Modified: branches/robkalman_pr/pkg/robKalman/man/recFilter.Rd
===================================================================
--- branches/robkalman_pr/pkg/robKalman/man/recFilter.Rd 2009-06-10 10:45:12 UTC (rev 31)
+++ branches/robkalman_pr/pkg/robKalman/man/recFilter.Rd 2009-06-12 17:21:07 UTC (rev 32)
@@ -13,27 +13,40 @@
}
\usage{
-recursiveFilter(Y, a, S, F, Q, Z, V,
- initSc=.cKinitstep, predSc=.cKpredstep, corrSc=.cKcorrstep,
- initSr=NULL, predSr=NULL, corrSr=NULL, ...)
+recursiveFilter(initSc = .cKinitstep, predSc = .cKpredstep,
+ corrSc = .cKcorrstep,
+ initSr = NULL, predSr = NULL, corrSr = NULL,
+ nsim = 0, seed = NULL, ..., dropRuns = TRUE)
KalmanFilter(Y, a, S, F, Q, Z, V)
rLSFilter(Y, a, S, F, Q, Z, V, b, norm=Euclideannorm)
ACMfilter(Y, a, S, F, Q, Z, V, s0, psi, apsi, bpsi, cpsi, flag)
}
\arguments{
- \item{a}{mean of the initial state}
- \item{S}{initial state covariance (see below)}
- \item{Z}{observation matrix (see below)}
- \item{V}{observation error covariance (see below)}
- \item{F}{innovation transition matrix (see below)}
- \item{Q}{innovation covariance (see below)}
+ \item{a}{mean of the initial state, in matrix form \eqn{p \times runs}{p x runs}
+ for \code{p} the state dimension
+ and \code{runs} the number of paths to be followed in parallel.}
+ \item{S}{initial state covariance (see below), in matrix form \eqn{p \times p}{p x p}
+ for \code{p} the state dimension.}
+ \item{Z}{observation matrices (see below), in array form \eqn{q \times p \times t}{q x p x t}
+ for \code{q} and \code{p} the observation resp.\ state dimension
+ and \code{t} the length of the observation series.}
+ \item{V}{observation error covariances (see below), in array form \eqn{q \times q \times t}{q x q x t}
+ for \code{q} the observation dimension
+ and \code{t} the length of the observation series.}
+ \item{F}{innovation transition matrices (see below), in array form \eqn{p \times p \times t}{p x p x t}
+ for \code{p} the state dimension
+ and \code{t} the length of the observation series.}
+ \item{Q}{innovation covariances (see below), in array form \eqn{p \times p \times t}{p x p x t}
+ for \code{p} the state dimension
+ and \code{t} the length of the observation series.}
\item{b}{clipping height \code{b} for the rLS filter}
\item{norm}{a function with a numeric vector \code{x} as first argument,
returning a norm of \code{x} - not necessarily, but defaulting to, Euclidean norm;
used by rLS filter to determine "too" large corrections}
- \item{Y}{observations \eqn{y_t}, in matrix form \eqn{q \times t}{q x t} for \code{t}
- the length of the observation series and \code{q} the observation dimension}
+ \item{Y}{observations \eqn{y_t}, in array form \eqn{q \times runs \times t}{q x runs x t} for \code{t}
+ the length of the observation series, \code{q} the observation dimension
+ and \code{runs} the number of paths to be followed in parallel.}
\item{s0}{scale of nominal Gaussian component of additive noise}
\item{psi}{influence function to be used (default: Hampel's \eqn{\psi} function, which is the only one available at the moment)}
\item{apsi,bpsi,cpsi}{tuning constants for Hampel's \eqn{\psi}-function, (default: \code{a=b=2.5}, \code{c=5.0})}
@@ -50,6 +63,9 @@
and \code{V} to produce the classical (non robust) correction value \eqn{x_{t|t}}}
\item{corrSr}{either \code{NULL} or a function with first arguments \code{y}\eqn{=y_t}, \code{x1}\eqn{=x_{t|t-1}} and \code{S1}\eqn{=S_{t|t-1}}, \code{Z},
and \code{V} to produce a robust correction value \eqn{x_{t|t}}}
+ \item{nsim}{integer; if positive, we simulate a bunch of nsim paths (acc. to ideal model) to get emp. covariances}
+ \item{seed}{seed for the simulations}
+ \item{dropruns}{logical; shall run-dimension be collapsed if it is one?}
\item{...}{further arguments to the "step"-functions}
}
@@ -60,10 +76,12 @@
\deqn{x_0 \sim {\cal N}_p(a,S),}{x_0 ~ N_p(a,S),}
\strong{state equation}
-\deqn{x_t = F x_{t-1} + v_t, \qquad v_t \sim {\cal N}_p(0,Q),\qquad t\ge 1,}{x_t = F x_{t-1} + v_t, v_t ~ N_p(0,Q), t>=1,}
+\deqn{x_t = F_t x_{t-1} + v_t, \qquad v_t \sim {\cal N}_p(0,Q_t),\qquad t\ge 1,
+ }{x_t = F_t x_{t-1} + v_t, v_t ~ N_p(0,Q_t), t>=1,}
\strong{observation equation}
-\deqn{y_t = Z x_{t} + \varepsilon_t, \qquad \varepsilon_t \sim {\cal N}_q(0,V),\qquad t\ge 1,}{y_t = Z x_t + e_t, e_t ~ N_q(0,V), t>=1,}
+\deqn{y_t = Z_t x_{t} + \varepsilon_t, \qquad \varepsilon_t \sim
+{\cal N}_q(0,V_t),\qquad t\ge 1,}{y_t = Z_t x_t + e_t, e_t ~ N_q(0,V_t), t>=1,}
and where all random variable \eqn{x_0}, \eqn{v_t}, \eqn{\varepsilon_t}{e_t} are independent.
For notation, let us formulate the classical Kalman filter in this context:
@@ -73,34 +91,46 @@
\deqn{S_{0|0} = {\rm Cov}(x_0-x_{0|0}) = \code{S}}{S_{0|0} = Cov(x_0-x_{0|0}) = S}%
\strong{(1) prediction step}
-\deqn{x_{t|t-1} = F x_{t-1|t-1},\qquad t\ge 1}{x_{t|t-1} = F x_{t-1|t-1}, t>=1}
+\deqn{x_{t|t-1} = F_t x_{t-1|t-1},\qquad t\ge 1}{x_{t|t-1} = F_t x_{t-1|t-1}, t>=1}
\eqn{\qquad}{\code{ }} with error covariance
-\deqn{S_{t|t-1} = {\rm Cov}(x_t-x_{t|t-1}) = F S_{t-1|t-1} F' + Q}{S_{t|t-1} = Cov(x_t-x_{t|t-1}) = F S_{t-1|t-1} F' + Q}
+\deqn{S_{t|t-1} = {\rm Cov}(x_t-x_{t|t-1}) =
+F_t S_{t-1|t-1} F_t' + Q_t}{S_{t|t-1} = Cov(x_t-x_{t|t-1}) = F_t S_{t-1|t-1} F_t' + Q_t}
\strong{(2) correction step}
-\deqn{x_{t|t} = x_{t|t-1} + K_t (y_t - Z x_{t|t-1}),\qquad t\ge 1}{x_{t|t} = x_{t|t-1} + K_t (y_t - Z x_{t|t-1}), t>=1}
+\deqn{x_{t|t} = x_{t|t-1} + K_t (y_t - Z_t x_{t|t-1}),\qquad
+t\ge 1}{x_{t|t} = x_{t|t-1} + K_t (y_t - Z_t x_{t|t-1}), t>=1}
\eqn{\qquad}{\code{ }} for Kalman Gain
-\deqn{K_t = S_{t|t-1} Z' (Z S_{t|t-1} Z' + V )^-}
+\deqn{K_t = S_{t|t-1} Z_t' (Z_t S_{t|t-1} Z_t' + V_t )^-}
\eqn{\qquad}{\code{ }} with error covariance
-\deqn{S_{t|t} = {\rm Cov}(x_t-x_{t|t}) = S_{t|t-1} - K_t Z S_{t|t-1}}{S_{t|t} = Cov(x_t-x_{t|t}) = S_{t|t-1} - K_t Z S_{t|t-1}}
+\deqn{S_{t|t} = {\rm Cov}(x_t-x_{t|t}) = S_{t|t-1} - K_t Z_t S_{t|t-1}}{S_{t|t} =
+Cov(x_t-x_{t|t}) = S_{t|t-1} - K_t Z_t S_{t|t-1}}
}
\value{
-\code{recursiveFilter} returns a list with the following elements
-\item{Xf:}{ the series \eqn{x_{t|t}} filtered by the classical filter --- a matrix with dimensions \eqn{p \times t+1}{p x (t+1)}}
-\item{Xp:}{ the series \eqn{x_{t|t-1}} predicted by the classical filter --- a matrix with dimensions \eqn{p \times t}{p x t}}
-\item{S0:}{ the series \eqn{S_{t|t}} of filter error covariances produced by the classical filter ---
- an array with dimensions \eqn{p \times p \times t+1}{p x (t+1)}}
+\item{Xf:}{ the series \eqn{x_{t|t}} filtered by the classical filter --- a matrix with
+ dimensions \eqn{p \times t+1}{p x (t+1)}}
+\item{Xp:}{ the series \eqn{x_{t|t-1}} predicted by the classical filter ---
+ a matrix with dimensions \eqn{p \times t}{p x t}}
+\item{Xrf:}{if any of the arguments \code{initSr}, \code{predSr}, \code{corrSr} is not \code{NULL}:
+ the series \eqn{x_{t|t}} filtered by the robust filter --- a matrix with dimensions
+ \eqn{p \times t+1}{p x (t+1)}
+ else \code{NULL}}
+\item{Xrp:}{if any of the arguments \code{initSr}, \code{predSr}, \code{corrSr} is not \code{NULL}:
+ the series \eqn{x_{t|t-1}} predicted by the robust filter
+ --- a matrix with dimensions \eqn{p \times t}{p x t}
+ else \code{NULL}}
+\item{S0:}{ the series \eqn{S_{t|t}} of filter error covariances produced by the classical filter ---
+ an array with dimensions \eqn{p \times p \times t+1}{p x p x (t+1)}}
\item{S1:}{ the series \eqn{S_{t|t-1}} of prediction error covariances produced by the classical filter ---
- an array with dimensions \eqn{p \times p \times t}{p x t}}
+ an array with dimensions \eqn{p \times p \times t}{p x p x t}}
\item{KG:}{ the series \eqn{K_{t}} of Kalman gains produced by the classical filter ---
- an array with dimensions \eqn{q \times p \times t}{p x t}}
-\item{Xrf:}{if any of the arguments \code{initSr}, \code{predSr}, \code{corrSr} is not \code{NULL}:
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/robkalman -r 32
More information about the Robkalman-commits
mailing list