[Robkalman-commits] r33 - branches/robkalman_pr/pkg/robKalman branches/robkalman_pr/pkg/robKalman/R branches/robkalman_pr/pkg/robKalman/demo branches/robkalman_pr/pkg/robKalman/man pkg/robKalman pkg/robKalman/R pkg/robKalman/chm pkg/robKalman/demo pkg/robKalman/man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Jun 15 15:33:26 CEST 2009
Author: ruckdeschel
Date: 2009-06-15 15:33:24 +0200 (Mon, 15 Jun 2009)
New Revision: 33
Modified:
branches/robkalman_pr/pkg/robKalman/DESCRIPTION
branches/robkalman_pr/pkg/robKalman/NAMESPACE
branches/robkalman_pr/pkg/robKalman/R/ACMfilt.R
branches/robkalman_pr/pkg/robKalman/R/ACMfilter.R
branches/robkalman_pr/pkg/robKalman/R/calibrateRLS.R
branches/robkalman_pr/pkg/robKalman/R/classKalman.R
branches/robkalman_pr/pkg/robKalman/R/rLSfilter.R
branches/robkalman_pr/pkg/robKalman/R/recFilter.R
branches/robkalman_pr/pkg/robKalman/R/simulateSScont.R
branches/robkalman_pr/pkg/robKalman/demo/rLSdemo.R
branches/robkalman_pr/pkg/robKalman/man/internalACM.Rd
branches/robkalman_pr/pkg/robKalman/man/internalrLS.Rd
branches/robkalman_pr/pkg/robKalman/man/recFilter.Rd
branches/robkalman_pr/pkg/robKalman/man/simulateSScont.Rd
pkg/robKalman/DESCRIPTION
pkg/robKalman/NAMESPACE
pkg/robKalman/R/ACMfilt.R
pkg/robKalman/R/ACMfilter.R
pkg/robKalman/R/calibrateRLS.R
pkg/robKalman/R/classKalman.R
pkg/robKalman/R/rLSfilter.R
pkg/robKalman/R/recFilter.R
pkg/robKalman/R/simulateSScont.R
pkg/robKalman/chm/00Index.html
pkg/robKalman/chm/internalACM.html
pkg/robKalman/chm/internalKalman.html
pkg/robKalman/chm/internalrLS.html
pkg/robKalman/chm/recFilter.html
pkg/robKalman/chm/robKalman.chm
pkg/robKalman/chm/simulateSScont.html
pkg/robKalman/demo/rLSdemo.R
pkg/robKalman/man/internalACM.Rd
pkg/robKalman/man/internalrLS.Rd
pkg/robKalman/man/recFilter.Rd
pkg/robKalman/man/simulateSScont.Rd
Log:
time-variant SSM algos (S4-free) uploaded; builds / checks error-free
Modified: branches/robkalman_pr/pkg/robKalman/DESCRIPTION
===================================================================
--- branches/robkalman_pr/pkg/robKalman/DESCRIPTION 2009-06-12 17:21:07 UTC (rev 32)
+++ branches/robkalman_pr/pkg/robKalman/DESCRIPTION 2009-06-15 13:33:24 UTC (rev 33)
@@ -1,6 +1,6 @@
Package: robKalman
-Version: 0.3
-Date: 2009-03-18
+Version: 0.4
+Date: 2009-06-15
Title: Robust Kalman Filtering
Description: Routines for Robust Kalman Filtering --- the ACM- and rLS-filter
Author: Peter Ruckdeschel, Bernhard Spangl, Irina Ursachi, Cezar Chirila
Modified: branches/robkalman_pr/pkg/robKalman/NAMESPACE
===================================================================
--- branches/robkalman_pr/pkg/robKalman/NAMESPACE 2009-06-12 17:21:07 UTC (rev 32)
+++ branches/robkalman_pr/pkg/robKalman/NAMESPACE 2009-06-15 13:33:24 UTC (rev 33)
@@ -1,5 +1,6 @@
import("methods")
import("stats")
+import("MASS")
import("startupmsg")
export("ACMfilt", "ACMfilter", "arGM", "Euclidnorm",
Modified: branches/robkalman_pr/pkg/robKalman/R/ACMfilt.R
===================================================================
--- branches/robkalman_pr/pkg/robKalman/R/ACMfilt.R 2009-06-12 17:21:07 UTC (rev 32)
+++ branches/robkalman_pr/pkg/robKalman/R/ACMfilt.R 2009-06-15 13:33:24 UTC (rev 33)
@@ -52,10 +52,10 @@
## Centering:
x <- x - gm$mu
- ACMres <- ACMfilter(Y = matrix(x,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=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)
### from version ... the return value $X[r]f is of
# dimension p x runs x (N+1)
Modified: branches/robkalman_pr/pkg/robKalman/R/ACMfilter.R
===================================================================
--- branches/robkalman_pr/pkg/robKalman/R/ACMfilter.R 2009-06-12 17:21:07 UTC (rev 32)
+++ branches/robkalman_pr/pkg/robKalman/R/ACMfilter.R 2009-06-15 13:33:24 UTC (rev 33)
@@ -28,7 +28,7 @@
list( x0 = a, S0 = S, s0 = s0)}
-.ACMpredstep <- function (x0, S0, F, Q, i, rob0, s0, ...) ### S=P F= Phi
+.ACMpredstep <- function (x0, S0, F, Q, i, rob0, s0, ...) ### S=P F= Phi
{
###########################################
##
@@ -83,7 +83,7 @@
rst <- (y - x1[1])/st
- ps <- psi(rst, apsi, bpsi, cpsi)
+ ps <- psi(rst, apsi, bpsi, cpsi)[1,1]
dx <- K * st * ps
x0 <- x1 + dx
@@ -93,5 +93,7 @@
S0 <- .getcorrCovACM(S1, K, Z, W = w*diag(rep(1, nrow(Z))))
- return(list(x0 = x0, K = K, S0 = S0, Delta=NULL, Ind=ind, rob0=rob1, DeltaY = rst))
+ Delta <- Z %*% S0 %*% t(Z) + V
+
+ return(list(x0 = x0, K = K, S0 = S0, Delta=Delta, Ind=ind, rob0=rob1, DeltaY = rst))
}
Modified: branches/robkalman_pr/pkg/robKalman/R/calibrateRLS.R
===================================================================
--- branches/robkalman_pr/pkg/robKalman/R/calibrateRLS.R 2009-06-12 17:21:07 UTC (rev 32)
+++ branches/robkalman_pr/pkg/robKalman/R/calibrateRLS.R 2009-06-15 13:33:24 UTC (rev 33)
@@ -19,25 +19,30 @@
dx <- t(mvrnorm(repl, numeric(pd), S))
dy <- Z %*% dx + t(mvrnorm(repl, numeric(qd), V))
- K <- .getKG(S, Z, V)
+ K <- .getKG(S, Z, .getDelta(S, Z, V))
trS <- sum(diag(.getcorrCov(S, K, Z)))
dx0 <- K %*% dy
no <- sqrt(t(ep) %*% dx0^2)
+ }
if (missing(r)) ## calibrated to given efficiency
- {f <- function(b, dX = dx, dX0 = dx0, no0 = no,
+ {f <- function(b, dX = dx, dX0 = dx0, no0 = no, r0=r
eff0 = eff, trS0 = trS, repl0 = repl)
{w <- ifelse(no0 < b, 1, b/no0)
dxw <- as.vector(w) * t(dX0)
trSb <- sum( (t(dX) - dxw)^2 )/repl0
trS0 / trSb - eff0
}
+ r1 <- NULL
+ eff1 <- eff
}
else ## calibrated to given radius
{f <- function(b, no0 = no, r0=r, repl0 = repl)
{(1 - r0)/r0 * sum(pmax(no0 / b - 1, 0))/repl0 - 1}
+ r1 <- r
+ eff1 <- NULL
}
erg <- uniroot(f, interval = c(10^-6, upto*sqrt(trS)), tol = 10^-7)
Modified: branches/robkalman_pr/pkg/robKalman/R/classKalman.R
===================================================================
--- branches/robkalman_pr/pkg/robKalman/R/classKalman.R 2009-06-12 17:21:07 UTC (rev 32)
+++ branches/robkalman_pr/pkg/robKalman/R/classKalman.R 2009-06-15 13:33:24 UTC (rev 33)
@@ -63,4 +63,4 @@
DeltaY <- y - Z %*% x1
x0 <- x1 + K %*% DeltaY
S0 <- .getcorrCov(S1, K, Z)
- list(x0 = x0, K = K, S0 = S0, Delta = Delta, Ind=1, DeltaY = yDelta)}
+ list(x0 = x0, K = K, S0 = S0, Delta = Delta, Ind=1, DeltaY = DeltaY)}
Modified: branches/robkalman_pr/pkg/robKalman/R/rLSfilter.R
===================================================================
--- branches/robkalman_pr/pkg/robKalman/R/rLSfilter.R 2009-06-12 17:21:07 UTC (rev 32)
+++ branches/robkalman_pr/pkg/robKalman/R/rLSfilter.R 2009-06-15 13:33:24 UTC (rev 33)
@@ -40,7 +40,7 @@
## (simultaneously for cK (x0, x1) )
.rLScorrstep <- function(y, x1, S1, Z, V, i, rob1 = NULL, b,
- norm = EuclideanNorm, ...)
+ norm = Euclideannorm, ...)
{Delta <- .getDelta(S1, Z, V)
K <- .getKG(S1, Z, Delta)
DeltaY <- y - Z %*% x1
Modified: branches/robkalman_pr/pkg/robKalman/R/recFilter.R
===================================================================
--- branches/robkalman_pr/pkg/robKalman/R/recFilter.R 2009-06-12 17:21:07 UTC (rev 32)
+++ branches/robkalman_pr/pkg/robKalman/R/recFilter.R 2009-06-15 13:33:24 UTC (rev 33)
@@ -62,7 +62,6 @@
# + dropRuns: shall run-dimension be collapsed if it is one?
{
qd <- ifelse(length(Z)==1, 1, (dim(Y))[1])
-
########################
# for backward compatibility
if (!is.array(Y)){
@@ -79,7 +78,13 @@
IndIO <- NULL
IndAO <- NULL
St0s <- St1s <- NULL
+ DeltaYr <- NULL
+ Deltar <- NULL
+ if(pd==1) F <- array(F, dim = c(pd,pd,tt))
+ if(pd==1 && qd==1) Z <- array(Z, dim = c(pd,qd,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(pd,qd,tt))
if(is.matrix(Q)) Q <- array(Q, dim = c(pd,pd,tt))
@@ -136,7 +141,7 @@
if(robust)
{
if(nsim){
- Xs <- t(rmvnorm(nsim, a, S))
+ Xs <- t(mvrnorm(nsim, a, S))
St0s <- array(0, c(pd, pd, tt))
St1s <- array(0, c(pd, pd, tt))
}
@@ -184,7 +189,7 @@
..., rob0 = rob0)
IndIO[,i] <- as.logical(psr$Ind)
if(nsim){
- vs <- t(rmvnorm(nsim, a*0, Q0))
+ 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
@@ -193,7 +198,7 @@
}else{
psr <- predSc(x0 = xr0, S0 = Sr0, F = F0, Q = Q0, i = i, ...)
if(nsim){
- vs <- t(rmvnorm(nsim, a*0, Q0))
+ vs <- t(mvrnorm(nsim, a*0, Q0))
Xs <- F %*% Xs + vs
xr1s <- predSc(x0 = xr0s, S0 = Sr0, F = F0, Q = Q0, i = i,
...)$x1
@@ -232,7 +237,7 @@
Z = Z0, V = V0, i = i, ..., rob1 = rob1)
IndAO[,i] <- as.logical(csr$Ind)
if(nsim){
- es <- t(rmvnorm(nsim, Y[,1]*0, V0))
+ 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
@@ -242,7 +247,7 @@
csr <- corrSc(y = Y0, x1 = xr1, S1 = Sr1, Z = Z0, V = V0,
i = i, ...)
if(nsim){
- es <- t(rmvnorm(nsim, Y[,1]*0, V0))
+ 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
@@ -252,13 +257,13 @@
xr0 <- csr$x0
Sr0 <- csr$S0
rob0 <- csr$rob0
- DeltaYr[,,i] <- cs$DeltaY
+ DeltaYr[,,i] <- csr$DeltaY
Xrf[,, i + 1] <- xr0
Str0[,, i + 1] <- S0
rob0L[[i + 1]] <- rob0
KGr[,, i] <- csr$K
- Deltar[,, i] <- cs$Deltar
+ Deltar[,, i] <- csr$Delta
}
}
@@ -287,50 +292,50 @@
# simple wrappers:
######################################################
-KalmanFilter <- function(Y, a, S, F, Q, Z, V)#
-#arguments:
+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)}
+{recursiveFilter(Y, a, S, F, Q, Z, V, dropRuns = dropRuns)}
-rLSFilter <- function(Y, a, S, F, Q, Z, V, b, norm = EuclideanNorm)#
-#arguments:
+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,
+{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)}
+ initSr = .cKinitstep, predSr = .cKpredstep,
+ corrSr = .rLScorrstep, b = b, norm = norm, dropRuns = dropRuns)}
-ACMfilter <- function(Y, a, S, F, Q, Z, V, s0, psi, apsi, bpsi, cpsi, flag)#
-#arguments:
+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
+# + b :clippingheight
+## Y=x ... observed time series
## a=m0 ... unconditional mean
## S=Cx ... covariance matrix
-## F=Phi ... design matrix of state equation
+## 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,
+## 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
+## 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, psi,
- apsi=2.5, bpsi=2.5, cpsi=5.0, flag)}
+{ 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
@@ -340,24 +345,24 @@
###########################################
## Paramters:
-## Y=x ... observed time series
+## Y=x ... observed time series
## a=m0 ... unconditional mean
## S=Cx ... covariance matrix
-## F=Phi ... design matrix of state equation
+## 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,
+## 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
+## 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")
+mACMfilter <- function(Y, a, S, F, Q, Z, V,
+ psi=mvpsiHampel, apsi=2.5, bpsi=2.5, cpsi=5.0,
+ flag="deriv", dropRuns = TRUE)
{
###########################################
##
@@ -368,26 +373,26 @@
###########################################
## Paramters:
-## Y ... observed vector-valued time series
-## (column-wise, matrix: q rows, number of columns equal to time points)
+## 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
+## 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,
+## 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
+## 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)
+ 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)
}
Modified: branches/robkalman_pr/pkg/robKalman/R/simulateSScont.R
===================================================================
--- branches/robkalman_pr/pkg/robKalman/R/simulateSScont.R 2009-06-12 17:21:07 UTC (rev 32)
+++ branches/robkalman_pr/pkg/robKalman/R/simulateSScont.R 2009-06-15 13:33:24 UTC (rev 33)
@@ -201,3 +201,40 @@
Ind.AO = erg$Ind.AO, Ind.IO = erg$Ind.IO))
})
+
+
+rcvmvnorm <- function(runs, mi, Si, mc, Sc, r)
+ {U<-rbinom(runs, size = 1, prob = r);
+ (1-U) * mvrnorm(runs, mi, Si) + U * mvrnorm(runs, mc, Sc)}
+
+simulateState <- function(a, S, F, Qi, mc=0, Qc=Qi, runs = 1, tt, r=0){
+ pd <- length(a)
+ if(length(dim(F))<3) F <- array(F, dim=c(pd,pd,tt))
+ if(!is.matrix(mc)) mc <- matrix(mc, pd,tt)
+ if(length(dim(Qi))<3) Qi <- array(Qi, dim=c(pd,pd,tt))
+ if(length(dim(Qc))<3) Qc <- array(Qc, dim=c(pd,pd,tt))
+ states <- array(0, dim=c(pd,runs, tt+1))
+ states[ ,,1] <- t(matrix(mvrnorm(runs, m = a, S = S),runs,pd))
+ for (i in (1:tt))
+ states[,, i+1] <- F[,,i] %*% states[,, i] +
+ t(matrix(rcvmvnorm(runs, mi = a*0, Si = Qi[,,i], mc = mc[,i],
+ Sc = Qc[,,i], r=r),runs,pd))
+ states
+}
+
+simulateObs <- function(X, Z, Vi, mc=0, Vc=Vi, runs = 1, r=0){
+ tt <- (dim(X))[3]-1
+ qd <- if(!is.null(dim(Vi))) (dim(Vi))[1] else 1
+ pd <- (dim(X))[1]
+ if(length(dim(Z))<3) Z <- array(Z, dim=c(qd,pd,tt))
+ if(!is.matrix(mc)) mc <- matrix(mc, qd,tt)
+ if(length(dim(Vi))<3) Vi <- array(Vi, dim=c(qd,qd,tt))
+ if(length(dim(Vc))<3) Vc <- array(Vc, dim=c(qd,qd,tt))
+ obs<-array(0, dim=c(qd, runs, tt))
+ for (i in 1:tt){
+ obs[,, i] <- Z[,,i] %*% X[,, (i+1)] +
+ t(matrix(rcvmvnorm(runs, 0*mc[,1], Vi[,,i], mc[,i],
+ Vc[,,i], r=r),runs,qd))
+ }
+ obs
+}
Modified: branches/robkalman_pr/pkg/robKalman/demo/rLSdemo.R
===================================================================
--- branches/robkalman_pr/pkg/robKalman/demo/rLSdemo.R 2009-06-12 17:21:07 UTC (rev 32)
+++ branches/robkalman_pr/pkg/robKalman/demo/rLSdemo.R 2009-06-15 13:33:24 UTC (rev 33)
@@ -10,10 +10,11 @@
{
DDa<-list(NULL)
if(ncoord==1)
- {DD <- as.data.frame(cbind("act. state"=t(X), "classK."=t(erg1$Xf),
+ {pd<- if(!is.null(dim(X))) dim(X)[1] else 1
+ DD <- as.data.frame(cbind("act. state"=t(matrix(X[,1,],nrow=pd)), "classK."=t(erg1$Xf),
t(erg1$Xrf), t(erg2$Xrf)))
names(DD)[3:4] <- c(nam1,nam2)
- if (withY) {yy<-c(NA,t(Y[1,])); DD <- as.data.frame(cbind(DD,"y"=yy))}
+ if (withY) {yy<-c(NA,t(Y[1,1,])); DD <- as.data.frame(cbind(DD,"y"=yy))}
if (withInd)
{
xclip<-c(FALSE,as.logical(erg1$IndAO))
@@ -23,7 +24,7 @@
}
else
for(coord in 1:ncoord)
- { DD <- as.data.frame(cbind("act. state"=X[coord,], "classK."=erg1$Xf[coord,],
+ { DD <- as.data.frame(cbind("act. state"=X[coord,1,], "classK."=erg1$Xf[coord,],
erg1$Xrf[coord,], erg2$Xrf[coord,]))
names(DD)[3:4] <- c(nam1,nam2)
if (withInd)
@@ -62,9 +63,9 @@
mc = m0c, Vc = V0c, r = ract, rcalib=r1, effcalib=eff1)
{
#Simulation::
-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)
+X <- simulateState(a = a0, S = Ss, F = F, Qi = 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)
pd <- dim(X)[1]
qd <- dim(Yid)[1]
@@ -93,21 +94,20 @@
###evaluation of MSE
###ideal situation
-MSEid <- c("class.Kalman"=mean((X - rerg1.id$Xf)^2), ### MSE averaged over time
- "rLS.eff"=mean((X - rerg1.id$Xrf)^2),
- "rLS.r"=mean((X - rerg2.id$Xrf)^2))
+MSEid <- c("class.Kalman"=mean((X[,1,] - rerg1.id$Xf)^2), ### MSE averaged over time
+ "rLS.eff"=mean((X[,1,] - rerg1.id$Xrf)^2),
+ "rLS.r"=mean((X[,1,] - rerg2.id$Xrf)^2))
###real situation
-MSEre <- c("class.Kalman"=mean((X - rerg1.re$Xf)^2), ### MSE averaged over time
- "rLS.eff"=mean((X - rerg1.re$Xrf)^2),
- "rLS.r"=mean((X - rerg2.re$Xrf)^2))
+MSEre <- c("class.Kalman"=mean((X[,1,] - rerg1.re$Xf)^2), ### MSE averaged over time
+ "rLS.eff"=mean((X[,1,] - rerg1.re$Xrf)^2),
+ "rLS.r"=mean((X[,1,] - rerg2.re$Xrf)^2))
print(list("Ideal situation"=MSEid,"Real situation"=MSEre))
op <- par(ask = interactive())
}
-
##############################################################
###Example1
##############################################################
Modified: branches/robkalman_pr/pkg/robKalman/man/internalACM.Rd
===================================================================
--- branches/robkalman_pr/pkg/robKalman/man/internalACM.Rd 2009-06-12 17:21:07 UTC (rev 32)
+++ branches/robkalman_pr/pkg/robKalman/man/internalACM.Rd 2009-06-15 13:33:24 UTC (rev 33)
@@ -14,8 +14,8 @@
\usage{
.getcorrCovACM(S1, K, Z, W=diag(nrow(Z)))
.ACMinitstep(a, S, ...)
-.ACMpredstep(x0, S0, F, Q, rob0, s0, ...)
-.ACMcorrstep(y, x1, S1, Z, V, rob1, dum, psi, apsi, bpsi, cpsi, flag, ...)
+.ACMpredstep(x0, S0, F, Q, i, rob0, s0, ...)
+.ACMcorrstep(y, x1, S1, Z, V, i, rob1, dum, psi, apsi, bpsi, cpsi, flag, ...)
}
@@ -26,6 +26,7 @@
\item{V}{observation error covariance (see below)}
\item{F}{innovation transition matrix (see below)}
\item{Q}{innovation covariance (see below)}
+ \item{i}{the time instance}
\item{K}{Kalman gain \eqn{K_t}}
\item{W}{weight matrix}
\item{dum}{dummy variable for compatibility with \dots argument of calling function}
Modified: branches/robkalman_pr/pkg/robKalman/man/internalrLS.Rd
===================================================================
--- branches/robkalman_pr/pkg/robKalman/man/internalrLS.Rd 2009-06-12 17:21:07 UTC (rev 32)
+++ branches/robkalman_pr/pkg/robKalman/man/internalrLS.Rd 2009-06-15 13:33:24 UTC (rev 33)
@@ -11,14 +11,15 @@
}
\usage{
-.rLScorrstep(y, x1, S1, Z, V, rob1=NULL, b, norm=Euclideannorm, ...)
+.rLScorrstep(y, x1, S1, Z, V, i, rob1=NULL, b, norm=Euclideannorm, ...)
}
\arguments{
\item{Z}{observation matrix (see below)}
\item{V}{observation error covariance (see below)}
\item{b}{clipping height \code{b} for the rLS filter}
- \item{norm}{a function with a numeric vector \code{x} as first argument,
+ \item{i}{the time instance}
+ \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{S1}{prediction error covariance \eqn{S_{t|t-1}} of the classical Kalman filter}
Modified: branches/robkalman_pr/pkg/robKalman/man/recFilter.Rd
===================================================================
--- branches/robkalman_pr/pkg/robKalman/man/recFilter.Rd 2009-06-12 17:21:07 UTC (rev 32)
+++ branches/robkalman_pr/pkg/robKalman/man/recFilter.Rd 2009-06-15 13:33:24 UTC (rev 33)
@@ -13,13 +13,13 @@
}
\usage{
-recursiveFilter(initSc = .cKinitstep, predSc = .cKpredstep,
+recursiveFilter(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)
-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)
+KalmanFilter(Y, a, S, F, Q, Z, V, dropRuns = TRUE)
+rLSFilter(Y, a, S, F, Q, Z, V, b, norm=Euclideannorm, dropRuns = TRUE)
+ACMfilter(Y, a, S, F, Q, Z, V, s0, psi, apsi, bpsi, cpsi, flag, dropRuns = TRUE)
}
\arguments{
@@ -65,7 +65,7 @@
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{dropRuns}{logical; shall run-dimension be collapsed if it is one?}
\item{...}{further arguments to the "step"-functions}
}
@@ -198,7 +198,7 @@
V0c <- 0.1
ract <- 0.1
-X <- simulateState( a = a0, S = SS0, F = F0, Q = Q0, tt = TT)
+X <- simulateState( a = a0, S = SS0, F = F0, Qi = Q0, tt = TT)
Y <- simulateObs(X = X, Z = Z0, Vi = V0i, mc = m0c, Vc = V0c, r = ract)
SS <- limitS(S = SS0, F = F0, Q = Q0, Z = Z0, V = V0i)
@@ -216,9 +216,9 @@
rerg1 <- rLSFilter(Y, a = a0, S = SS0, F = F0, Q = Q0, Z = Z0, V = V0i, b = B1$b)
rerg2 <- rLSFilter(Y, a = a0, S = SS0, F = F0, Q = Q0, Z = Z0, V = V0i, b = B2$b)
-mean((X - rerg1$Xf)^2) ### empirical MSE of the filters considered
-mean((X - rerg1$Xrf)^2)
-mean((X - rerg2$Xrf)^2)
+mean((X[,1,] - rerg1$Xf)^2) ### empirical MSE of the filters considered
+mean((X[,1,] - rerg1$Xrf)^2)
+mean((X[,1,] - rerg2$Xrf)^2)
}
Modified: branches/robkalman_pr/pkg/robKalman/man/simulateSScont.Rd
===================================================================
--- branches/robkalman_pr/pkg/robKalman/man/simulateSScont.Rd 2009-06-12 17:21:07 UTC (rev 32)
+++ branches/robkalman_pr/pkg/robKalman/man/simulateSScont.Rd 2009-06-15 13:33:24 UTC (rev 33)
@@ -11,12 +11,13 @@
multivariate time-invariant, linear, Gaussian state space model may be generated
}
\usage{
-rcvmvnorm(mi, Si, mc, Sc, r)
-simulateState(a, S, F, Q, tt)
-simulateObs(X, Z, Vi, mc, Vc, r)
+rcvmvnorm(runs, mi, Si, mc, Sc, r)
+simulateState(a, S, F, Qi, mc=0, Qc=Qi, runs = 1, tt, r=0)
+simulateObs(X, Z, Vi, mc=0, Vc=Vi, runs = 1, r=0)
}
\arguments{
+ \item{runs}{number of runs to be generated}
\item{mi}{mean of the ideal multivariate normal distribution}
\item{mc}{mean of the contaminating multivariate normal distribution}
\item{Si}{covariance of the ideal multivariate normal distribution}
@@ -25,7 +26,8 @@
\item{a}{mean of the initial state}
\item{S}{initial state covariance (see below)}
\item{F}{innovation transition matrix (see below)}
- \item{Q}{innovation covariance (see below)}
+ \item{Qi}{ideal innovation covariance (see below)}
+ \item{Qc}{contaminating innovation covariance (see below)}
\item{tt}{length of the simulated series of states/observations}
\item{Z}{observation matrix (see below)}
\item{Vi}{ideal observation error covariance (see below)}
@@ -109,7 +111,7 @@
V0c <- 0.1
ract <- 0.1
-X <- simulateState( a = a0, S = SS0, F = F0, Q = Q0, tt = TT)
+X <- simulateState( a = a0, S = SS0, F = F0, Qi = Q0, tt = TT)
Y <- simulateObs(X = X, Z = Z0, Vi = V0i, mc = m0c, Vc = V0c, r = ract)
}
Modified: pkg/robKalman/DESCRIPTION
===================================================================
--- pkg/robKalman/DESCRIPTION 2009-06-12 17:21:07 UTC (rev 32)
+++ pkg/robKalman/DESCRIPTION 2009-06-15 13:33:24 UTC (rev 33)
@@ -1,6 +1,6 @@
Package: robKalman
-Version: 0.2.1
-Date: 2009-03-18
+Version: 0.3
+Date: 2009-06-14
Title: Robust Kalman Filtering
Description: Routines for Robust Kalman Filtering --- the ACM- and rLS-filter
Author: Peter Ruckdeschel, Bernhard Spangl
Modified: pkg/robKalman/NAMESPACE
===================================================================
--- pkg/robKalman/NAMESPACE 2009-06-12 17:21:07 UTC (rev 32)
+++ pkg/robKalman/NAMESPACE 2009-06-15 13:33:24 UTC (rev 33)
@@ -1,5 +1,6 @@
import("methods")
import("stats")
+import("MASS")
import("startupmsg")
export("ACMfilt", "ACMfilter", "arGM", "Euclideannorm",
Modified: pkg/robKalman/R/ACMfilt.R
===================================================================
--- pkg/robKalman/R/ACMfilt.R 2009-06-12 17:21:07 UTC (rev 32)
+++ pkg/robKalman/R/ACMfilt.R 2009-06-15 13:33:24 UTC (rev 33)
@@ -52,7 +52,10 @@
## Centering:
x <- x - gm$mu
- ACMres <- ACMfilter(Y=matrix(x,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=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)
X.ck <- ACMres$Xf; X.ck <- X.ck[,2:(N+1)]
X <- ACMres$Xrf; X <- X[,2:(N+1)]
Modified: pkg/robKalman/R/ACMfilter.R
===================================================================
--- pkg/robKalman/R/ACMfilter.R 2009-06-12 17:21:07 UTC (rev 32)
+++ pkg/robKalman/R/ACMfilter.R 2009-06-15 13:33:24 UTC (rev 33)
@@ -27,8 +27,7 @@
s0<-NULL
list( x0 = a, S0 = S, s0 = s0)}
-
-.ACMpredstep <- function (x0, S0, F, Q, rob0, s0, ...) ### S=P F= Phi
+.ACMpredstep <- function (x0, S0, F, Q, i, rob0, s0, ...) ### S=P F= Phi
{
###########################################
##
@@ -48,7 +47,7 @@
return(list(x1 = F %*% x0, S1 = S1, rob1 = sqrt(S1[1, 1] + s0), Ind=1))
}
-.ACMcorrstep <- function (y, x1, S1, Z, V, rob1, dum=NULL, psi, apsi, bpsi, cpsi, flag, ...)
+.ACMcorrstep <- function (y, x1, S1, Z, V, i, rob1, dum=NULL, psi, apsi, bpsi, cpsi, flag, ...)
{
###########################################
##
@@ -71,14 +70,13 @@
## (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)
-
st <- rob1
K <- .getKG(S1, Z, V)
rst <- (y - x1[1])/st
- ps <- psi(rst, apsi, bpsi, cpsi)
+ ps <- psi(rst, apsi, bpsi, cpsi)[1,1]
dx <- K * st * ps
x0 <- x1 + dx
@@ -87,6 +85,7 @@
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
- return(list(x0 = x0, K = K, S0 = S0, Delta=NULL, Ind=ind, rob0=rob1, DeltaY = rst))
+ return(list(x0 = x0, K = K, S0 = S0, Delta=Delta, Ind=ind, rob0=rob1, DeltaY = rst))
}
Modified: pkg/robKalman/R/calibrateRLS.R
===================================================================
--- pkg/robKalman/R/calibrateRLS.R 2009-06-12 17:21:07 UTC (rev 32)
+++ pkg/robKalman/R/calibrateRLS.R 2009-06-15 13:33:24 UTC (rev 33)
@@ -19,7 +19,7 @@
dx <- t(mvrnorm(repl, numeric(pd), S))
dy <- Z %*% dx + t(mvrnorm(repl, numeric(qd), V))
- K <- .getKG(S, Z, V)
+ K <- .getKG(S, Z, .getDelta(S, Z, V))
trS <- sum(diag(.getcorrCov(S, K, Z)))
Modified: pkg/robKalman/R/classKalman.R
===================================================================
--- pkg/robKalman/R/classKalman.R 2009-06-12 17:21:07 UTC (rev 32)
+++ pkg/robKalman/R/classKalman.R 2009-06-15 13:33:24 UTC (rev 33)
@@ -62,4 +62,4 @@
DeltaY <- y - Z %*% x1
x0 <- x1 + K %*% DeltaY
S0 <- .getcorrCov(S1, K, Z)
- list(x0 = x0, K = K, S0 = S0, Delta = Delta, Ind=1, DeltaY = yDelta)}
+ list(x0 = x0, K = K, S0 = S0, Delta = Delta, Ind=1, DeltaY = DeltaY)}
Modified: pkg/robKalman/R/rLSfilter.R
===================================================================
--- pkg/robKalman/R/rLSfilter.R 2009-06-12 17:21:07 UTC (rev 32)
+++ pkg/robKalman/R/rLSfilter.R 2009-06-15 13:33:24 UTC (rev 33)
@@ -40,7 +40,7 @@
## (simultaneously for cK (x0, x1) )
.rLScorrstep <- function(y, x1, S1, Z, V, i, rob1 = NULL, b,
- norm = EuclideanNorm, ...)
+ norm = Euclideannorm, ...)
{Delta <- .getDelta(S1, Z, V)
K <- .getKG(S1, Z, Delta)
DeltaY <- y - Z %*% x1
Modified: pkg/robKalman/R/recFilter.R
===================================================================
--- pkg/robKalman/R/recFilter.R 2009-06-12 17:21:07 UTC (rev 32)
+++ pkg/robKalman/R/recFilter.R 2009-06-15 13:33:24 UTC (rev 33)
@@ -20,7 +20,6 @@
# + dropRuns: shall run-dimension be collapsed if it is one?
{
qd <- ifelse(length(Z)==1, 1, (dim(Y))[1])
-
########################
# for backward compatibility
if (!is.array(Y)){
@@ -37,7 +36,13 @@
IndIO <- NULL
IndAO <- NULL
St0s <- St1s <- NULL
+ DeltaYr <- NULL
+ Deltar <- NULL
+ if(pd==1) F <- array(F, dim = c(pd,pd,tt))
+ if(pd==1 && qd==1) Z <- array(Z, dim = c(pd,qd,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(pd,qd,tt))
if(is.matrix(Q)) Q <- array(Q, dim = c(pd,pd,tt))
@@ -94,7 +99,7 @@
if(robust)
{
if(nsim){
- Xs <- t(rmvnorm(nsim, a, S))
+ Xs <- t(mvrnorm(nsim, a, S))
St0s <- array(0, c(pd, pd, tt))
St1s <- array(0, c(pd, pd, tt))
}
@@ -142,7 +147,7 @@
..., rob0 = rob0)
IndIO[,i] <- as.logical(psr$Ind)
if(nsim){
- vs <- t(rmvnorm(nsim, a*0, Q0))
+ 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
@@ -151,7 +156,7 @@
}else{
psr <- predSc(x0 = xr0, S0 = Sr0, F = F0, Q = Q0, i = i, ...)
if(nsim){
- vs <- t(rmvnorm(nsim, a*0, Q0))
+ vs <- t(mvrnorm(nsim, a*0, Q0))
Xs <- F %*% Xs + vs
xr1s <- predSc(x0 = xr0s, S0 = Sr0, F = F0, Q = Q0, i = i,
...)$x1
@@ -190,7 +195,7 @@
Z = Z0, V = V0, i = i, ..., rob1 = rob1)
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/robkalman -r 33
More information about the Robkalman-commits
mailing list