[Robkalman-commits] r34 - branches/robkalman_pr/pkg/robKalman branches/robkalman_pr/pkg/robKalman/R branches/robkalman_pr/pkg/robKalman/man pkg/robKalman pkg/robKalman/R pkg/robKalman/chm pkg/robKalman/man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Jun 19 14:50:03 CEST 2009
Author: ruckdeschel
Date: 2009-06-19 14:50:02 +0200 (Fri, 19 Jun 2009)
New Revision: 34
Modified:
branches/robkalman_pr/pkg/robKalman/NAMESPACE
branches/robkalman_pr/pkg/robKalman/R/calibrateRLS.R
branches/robkalman_pr/pkg/robKalman/R/rLSfilter.R
branches/robkalman_pr/pkg/robKalman/R/recFilter.R
branches/robkalman_pr/pkg/robKalman/man/0robKalman-package.Rd
branches/robkalman_pr/pkg/robKalman/man/calibrateRLS.Rd
branches/robkalman_pr/pkg/robKalman/man/recFilter.Rd
pkg/robKalman/NAMESPACE
pkg/robKalman/R/calibrateRLS.R
pkg/robKalman/R/rLSfilter.R
pkg/robKalman/R/recFilter.R
pkg/robKalman/chm/00Index.html
pkg/robKalman/chm/0robKalman-package.html
pkg/robKalman/chm/calibrateRLS.html
pkg/robKalman/chm/recFilter.html
pkg/robKalman/chm/robKalman.chm
pkg/robKalman/chm/robKalman.toc
pkg/robKalman/man/0robKalman-package.Rd
pkg/robKalman/man/calibrateRLS.Rd
pkg/robKalman/man/recFilter.Rd
Log:
Added IO-robust rLS
Modified: branches/robkalman_pr/pkg/robKalman/NAMESPACE
===================================================================
--- branches/robkalman_pr/pkg/robKalman/NAMESPACE 2009-06-15 13:33:24 UTC (rev 33)
+++ branches/robkalman_pr/pkg/robKalman/NAMESPACE 2009-06-19 12:50:02 UTC (rev 34)
@@ -3,7 +3,8 @@
import("MASS")
import("startupmsg")
-export("ACMfilt", "ACMfilter", "arGM", "Euclidnorm",
+export("ACMfilt", "ACMfilter", "arGM", "Euclideannorm",
"simulateState", "simulateObs", "rcvmvnorm", "Huberize",
- "rLScalibrateB", "limitS", "rLSFilter", "KalmanFilter",
+ "rLScalibrateB", "limitS", "rLSFilter", "rLS.IO.Filter",
+ "rLS.AO.Filter", "KalmanFilter",
"recursiveFilter")
Modified: branches/robkalman_pr/pkg/robKalman/R/calibrateRLS.R
===================================================================
--- branches/robkalman_pr/pkg/robKalman/R/calibrateRLS.R 2009-06-15 13:33:24 UTC (rev 33)
+++ branches/robkalman_pr/pkg/robKalman/R/calibrateRLS.R 2009-06-19 12:50:02 UTC (rev 34)
@@ -1,4 +1,4 @@
-rLScalibrateB <- function(Z, S, V, repl = 100000, eff, r, upto=20)#
+rLScalibrateB <- function(Z, S, V, repl = 100000, eff, r, upto=20, IO = FALSE)#
# calibrates clipping height b to given Z, V, and S_{t|t-1}
# ---
# either to given efficiency in the ideal model
@@ -15,8 +15,6 @@
qd <- ifelse(length(Z)==1, 1, (dim(Z))[1])
pd <- ifelse(length(S)==1, 1, (dim(Z))[2])
- ep <- 1+numeric(pd)
-
dx <- t(mvrnorm(repl, numeric(pd), S))
dy <- Z %*% dx + t(mvrnorm(repl, numeric(qd), V))
K <- .getKG(S, Z, .getDelta(S, Z, V))
@@ -24,28 +22,37 @@
trS <- sum(diag(.getcorrCov(S, K, Z)))
dx0 <- K %*% dy
- no <- sqrt(t(ep) %*% dx0^2)
- }
+ if(IO){
+ dx0 <- dy - Z %*% dx0
+ }
+ no <- sqrt(colSums(dx0^2))
+
+
if (missing(r)) ## calibrated to given efficiency
- {f <- function(b, dX = dx, dX0 = dx0, no0 = no, r0=r
- eff0 = eff, trS0 = trS, repl0 = repl)
+ {f <- function(b, dX = dx, dX0 = dx0, no0 = no, r0 = r,
+ eff0 = eff, trS0 = trS, repl0 = repl, dY = dy)
{w <- ifelse(no0 < b, 1, b/no0)
dxw <- as.vector(w) * t(dX0)
+ if(IO) dxw <- (t(dY) - dxw) %*% t(ginv(Z))
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)
+ {f <- function(b, dX = dx, dX0 = dx0, no0 = no, r0 = r,
+ eff0 = eff, trS0 = trS, repl0 = repl, dY = dy)
{(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)
+ erg <- uniroot(f, interval = c(10^-6, upto*sqrt(trS)), tol = 10^-7,
+ dX = dx, dX0 = dx0, no0 = no,
+ eff0 = eff1, trS0 = trS, repl0 = repl, r0 = r1,
+ dY = dy)
b <- erg$root
if (missing(r)) ### corresponding radius is calculated
@@ -54,6 +61,7 @@
else ### corresponding effciency is calculated
{ w <- ifelse(no < b, 1, b / no)
dxw <- as.vector(w) * t(dx0)
+ if(IO) dxw <- (t(dy) - dxw) %*% t(ginv(Z))
trSb <- sum( (t(dx) - dxw)^2 )/repl
eff <- trS / trSb
}
Modified: branches/robkalman_pr/pkg/robKalman/R/rLSfilter.R
===================================================================
--- branches/robkalman_pr/pkg/robKalman/R/rLSfilter.R 2009-06-15 13:33:24 UTC (rev 33)
+++ branches/robkalman_pr/pkg/robKalman/R/rLSfilter.R 2009-06-19 12:50:02 UTC (rev 34)
@@ -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}
#########################################################################################
@@ -54,3 +54,21 @@
else
Ind <- apply(dx, 2, function(xx) norm(xx)>1)
list(x0 = x0, K = K, S0 = S0, Delta = Delta, Ind = Ind, DeltaY = DeltaY)}
+
+
+.rLS.IO.corrstep <- 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
+ de <- DeltaY-Z%*%dx
+ if(length(b)>1)
+ b <- b[min(i,length(b))]
+ x0 <- x1 + ginv(Z)%*%(DeltaY-Huberize(de, b, norm = norm))
+ S0 <- .getcorrCov(S1, K, Z)
+ if (ncol(x1)==1)
+ Ind <- (norm(de)>1)
+ else
+ Ind <- apply(de, 2, function(xx) norm(xx)>1)
+ 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-15 13:33:24 UTC (rev 33)
+++ branches/robkalman_pr/pkg/robKalman/R/recFilter.R 2009-06-19 12:50:02 UTC (rev 34)
@@ -310,7 +310,24 @@
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
Modified: branches/robkalman_pr/pkg/robKalman/man/0robKalman-package.Rd
===================================================================
--- branches/robkalman_pr/pkg/robKalman/man/0robKalman-package.Rd 2009-06-15 13:33:24 UTC (rev 33)
+++ branches/robkalman_pr/pkg/robKalman/man/0robKalman-package.Rd 2009-06-19 12:50:02 UTC (rev 34)
@@ -12,7 +12,7 @@
\details{
\tabular{ll}{
Package: \tab robKalman\cr
-Version: \tab 0.3 \cr
+Version: \tab 0.2.1 \cr
Date: \tab 2009-03-18 \cr
Depends: \tab R(>= 2.3.0), methods, graphics, startupmsg, dse1, dse2, MASS\cr
Imports: \tab stats\cr
@@ -46,8 +46,11 @@
general recursive filters
+recursiveFilter
-KalmanFilter
- -rLSFilter
+ -rLSFilter:
+ *rLS.AO.Filter
+ *rLS.IO.Filter
-ACMfilter
+ -mACMfilter
ACMfilter:
+ACMfilt
Modified: branches/robkalman_pr/pkg/robKalman/man/calibrateRLS.Rd
===================================================================
--- branches/robkalman_pr/pkg/robKalman/man/calibrateRLS.Rd 2009-06-15 13:33:24 UTC (rev 33)
+++ branches/robkalman_pr/pkg/robKalman/man/calibrateRLS.Rd 2009-06-19 12:50:02 UTC (rev 34)
@@ -7,7 +7,7 @@
calibrates the clipping height \code{b} of the rLS-filter in a time-invariant, linear, Gaussian state space model
}
\usage{
-rLScalibrateB(Z, S, V, repl = 100000, eff, r, upto=20)
+rLScalibrateB(Z, S, V, repl = 100000, eff, r, upto=20, IO = FALSE)
}
\arguments{
\item{Z}{observation matrix in the (ti-l-G-SSM); see below}
@@ -17,6 +17,7 @@
\item{eff}{efficiency w.r.t. classical Kalman filter in the ideal model}
\item{repl}{number of replicates used for a LLN-approximation of the expectations needed in this calibration}
\item{upto}{an upper bound to \code{b} used in the zero-search of \code{uniroot} within \code{rLScalibrateB}}
+ \item{IO}{logical of length 1: Is it rLS.IO (\code{TRUE}) or rLS[.AO] which is to be calibrated?}
}
\details{
We work in the setup of the time-invariant, linear, Gaussian state space model (ti-l-G-SSM)
Modified: branches/robkalman_pr/pkg/robKalman/man/recFilter.Rd
===================================================================
--- branches/robkalman_pr/pkg/robKalman/man/recFilter.Rd 2009-06-15 13:33:24 UTC (rev 33)
+++ branches/robkalman_pr/pkg/robKalman/man/recFilter.Rd 2009-06-19 12:50:02 UTC (rev 34)
@@ -1,6 +1,8 @@
\name{recFilter}
\alias{recFilter}
\alias{rLSFilter}
+\alias{rLS.AO.Filter}
+\alias{rLS.IO.Filter}
\alias{KalmanFilter}
\alias{ACMfilter}
\alias{recursiveFilter}
@@ -19,6 +21,8 @@
nsim = 0, seed = NULL, ..., dropRuns = TRUE)
KalmanFilter(Y, a, S, F, Q, Z, V, dropRuns = TRUE)
rLSFilter(Y, a, S, F, Q, Z, V, b, norm=Euclideannorm, dropRuns = TRUE)
+rLS.AO.Filter(Y, a, S, F, Q, Z, V, b, norm=Euclideannorm, dropRuns = TRUE)
+rLS.IO.Filter(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)
}
@@ -160,11 +164,17 @@
e.g. in case of the ACM filter each element contains a corresponding value of \code{st}}
\code{KalmanFilter(Y, a, S, F, Q, Z, V)} is a wrapper to \code{recursiveFilter(Y, a, S, F, Q, Z, V)}.\cr
-\code{rLSFilter(Y, a, S, F, Q, Z, V, b, norm=Euclideannorm)} is a wrapper to\cr
+\code{rLS.AO.Filter(Y, a, S, F, Q, Z, V, b, norm=Euclideannorm)} and (synonymously)
+\code{rLSFilter(Y, a, S, F, Q, Z, V, b, norm=Euclideannorm)} are wrappers to\cr
\code{recursiveFilter(Y, a, S, F, Q, Z, V,
initSc=.cKinitstep, predSc=.cKpredstep, corrSc=.cKcorrstep,
- initSr=.cKinitstep, predSr=.cKpredstep, corrSr=rLScorrstep,
+ initSr=.cKinitstep, predSr=.cKpredstep, corrSr=.rLScorrstep,
b=b, norm=norm)}.
+\code{rLS.IO.Filter(Y, a, S, F, Q, Z, V, b, norm=Euclideannorm)} is a wrapper to\cr
+\code{recursiveFilter(Y, a, S, F, Q, Z, V,
+ initSc=.cKinitstep, predSc=.cKpredstep, corrSc=.cKcorrstep,
+ initSr=.cKinitstep, predSr=.cKpredstep, corrSr=.rLS.IO.corrstep,
+ b=b, norm=norm)}.
\code{ACMFilter(Y, a, S, F, Q, Z, V, s0, psi, apsi, bpsi, cpsi, flag)} is a wrapper to \cr
\code{recursiveFilter(Y, a, S, F, Q, Z, V,
initSc=.cKinitstep, predSc=.cKpredstep, corrSc=.cKcorrstep,
@@ -211,14 +221,28 @@
# r = 0.1
(B2 <- rLScalibrateB(r = 0.1, S = SS, Z = Z0, V = V0i))
+# IO-filter
+# by efficiency in the ideal model
+# efficiency = 0.9
+(B1.IO <- rLScalibrateB(eff = 0.9, S = SS, Z = Z0, V = V0i, IO = TRUE))
+# by contamination radius
+# r = 0.1
+(B2.IO <- rLScalibrateB(r = 0.1, S = SS, Z = Z0, V = V0i, IO = TRUE))
erg <- KalmanFilter(Y, a = a0, S = SS0, F = F0, Q = Q0, Z = Z0, V = V0i)
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)
+
+rerg1.IO <- rLS.IO.Filter(Y, a = a0, S = SS0, F = F0, Q = Q0, Z = Z0, V = V0i, b = B1.IO$b)
+rerg2.IO <- rLS.IO.Filter(Y, a = a0, S = SS0, F = F0, Q = Q0, Z = Z0, V = V0i, b = B2.IO$b)
+
mean((X[,1,] - rerg1$Xf)^2) ### empirical MSE of the filters considered
mean((X[,1,] - rerg1$Xrf)^2)
mean((X[,1,] - rerg2$Xrf)^2)
+### not surprisingly IO-rob filter is not a good idea for AO - situation
+mean((X[,1,] - rerg1.IO$Xrf)^2)
+mean((X[,1,] - rerg2.IO$Xrf)^2)
}
Modified: pkg/robKalman/NAMESPACE
===================================================================
--- pkg/robKalman/NAMESPACE 2009-06-15 13:33:24 UTC (rev 33)
+++ pkg/robKalman/NAMESPACE 2009-06-19 12:50:02 UTC (rev 34)
@@ -5,5 +5,6 @@
export("ACMfilt", "ACMfilter", "arGM", "Euclideannorm",
"simulateState", "simulateObs", "rcvmvnorm", "Huberize",
- "rLScalibrateB", "limitS", "rLSFilter", "KalmanFilter",
+ "rLScalibrateB", "limitS", "rLSFilter", "rLS.IO.Filter",
+ "rLS.AO.Filter", "KalmanFilter",
"recursiveFilter")
Modified: pkg/robKalman/R/calibrateRLS.R
===================================================================
--- pkg/robKalman/R/calibrateRLS.R 2009-06-15 13:33:24 UTC (rev 33)
+++ pkg/robKalman/R/calibrateRLS.R 2009-06-19 12:50:02 UTC (rev 34)
@@ -1,4 +1,4 @@
-rLScalibrateB <- function(Z, S, V, repl = 100000, eff, r, upto=20)#
+rLScalibrateB <- function(Z, S, V, repl = 100000, eff, r, upto=20, IO = FALSE)#
# calibrates clipping height b to given Z, V, and S_{t|t-1}
# ---
# either to given efficiency in the ideal model
@@ -15,8 +15,6 @@
qd <- ifelse(length(Z)==1, 1, (dim(Z))[1])
pd <- ifelse(length(S)==1, 1, (dim(Z))[2])
- ep <- 1+numeric(pd)
-
dx <- t(mvrnorm(repl, numeric(pd), S))
dy <- Z %*% dx + t(mvrnorm(repl, numeric(qd), V))
K <- .getKG(S, Z, .getDelta(S, Z, V))
@@ -24,13 +22,19 @@
trS <- sum(diag(.getcorrCov(S, K, Z)))
dx0 <- K %*% dy
- no <- sqrt(t(ep) %*% dx0^2)
+ if(IO){
+ dx0 <- dy - Z %*% dx0
+ }
+ no <- sqrt(colSums(dx0^2))
+
+
if (missing(r)) ## calibrated to given efficiency
{f <- function(b, dX = dx, dX0 = dx0, no0 = no, r0 = r,
- eff0 = eff, trS0 = trS, repl0 = repl)
+ eff0 = eff, trS0 = trS, repl0 = repl, dY = dy)
{w <- ifelse(no0 < b, 1, b/no0)
dxw <- as.vector(w) * t(dX0)
+ if(IO) dxw <- (t(dY) - dxw) %*% t(ginv(Z))
trSb <- sum( (t(dX) - dxw)^2 )/repl0
trS0 / trSb - eff0
}
@@ -38,8 +42,8 @@
eff1 <- eff
}
else ## calibrated to given radius
- {f <- function(b, dX = dx, dX0 = dx0, no0 = no, r0=r,
- eff0 = eff, trS0 = trS, repl0 = repl)
+ {f <- function(b, dX = dx, dX0 = dx0, no0 = no, r0 = r,
+ eff0 = eff, trS0 = trS, repl0 = repl, dY = dy)
{(1 - r0)/r0 * sum(pmax(no0 / b - 1, 0))/repl0 - 1}
r1 <- r
eff1 <- NULL
@@ -47,7 +51,8 @@
erg <- uniroot(f, interval = c(10^-6, upto*sqrt(trS)), tol = 10^-7,
dX = dx, dX0 = dx0, no0 = no,
- eff0 = eff1, trS0 = trS, repl0 = repl, r0 = r1)
+ eff0 = eff1, trS0 = trS, repl0 = repl, r0 = r1,
+ dY = dy)
b <- erg$root
if (missing(r)) ### corresponding radius is calculated
@@ -56,6 +61,7 @@
else ### corresponding effciency is calculated
{ w <- ifelse(no < b, 1, b / no)
dxw <- as.vector(w) * t(dx0)
+ if(IO) dxw <- (t(dy) - dxw) %*% t(ginv(Z))
trSb <- sum( (t(dx) - dxw)^2 )/repl
eff <- trS / trSb
}
Modified: pkg/robKalman/R/rLSfilter.R
===================================================================
--- pkg/robKalman/R/rLSfilter.R 2009-06-15 13:33:24 UTC (rev 33)
+++ pkg/robKalman/R/rLSfilter.R 2009-06-19 12:50:02 UTC (rev 34)
@@ -54,3 +54,21 @@
else
Ind <- apply(dx, 2, function(xx) norm(xx)>1)
list(x0 = x0, K = K, S0 = S0, Delta = Delta, Ind = Ind, DeltaY = DeltaY)}
+
+
+.rLS.IO.corrstep <- 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
+ de <- DeltaY-Z%*%dx
+ if(length(b)>1)
+ b <- b[min(i,length(b))]
+ x0 <- x1 + ginv(Z)%*%(DeltaY-Huberize(de, b, norm = norm))
+ S0 <- .getcorrCov(S1, K, Z)
+ if (ncol(x1)==1)
+ Ind <- (norm(de)>1)
+ else
+ Ind <- apply(de, 2, function(xx) norm(xx)>1)
+ list(x0 = x0, K = K, S0 = S0, Delta = Delta, Ind = Ind, DeltaY = DeltaY)}
Modified: pkg/robKalman/R/recFilter.R
===================================================================
--- pkg/robKalman/R/recFilter.R 2009-06-15 13:33:24 UTC (rev 33)
+++ pkg/robKalman/R/recFilter.R 2009-06-19 12:50:02 UTC (rev 34)
@@ -268,7 +268,24 @@
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
Modified: pkg/robKalman/chm/00Index.html
===================================================================
--- pkg/robKalman/chm/00Index.html 2009-06-15 13:33:24 UTC (rev 33)
+++ pkg/robKalman/chm/00Index.html 2009-06-19 12:50:02 UTC (rev 34)
@@ -35,6 +35,10 @@
<td>Several recursive filters: the classical Kalman filter, the rLS filter, and the ACM filter</td></tr>
<tr><td width="25%"><a href="recFilter.html">recursiveFilter</a></td>
<td>Several recursive filters: the classical Kalman filter, the rLS filter, and the ACM filter</td></tr>
+<tr><td width="25%"><a href="recFilter.html">rLS.AO.Filter</a></td>
+<td>Several recursive filters: the classical Kalman filter, the rLS filter, and the ACM filter</td></tr>
+<tr><td width="25%"><a href="recFilter.html">rLS.IO.Filter</a></td>
+<td>Several recursive filters: the classical Kalman filter, the rLS filter, and the ACM filter</td></tr>
<tr><td width="25%"><a href="calibrateRLS.html">rLScalibrateB</a></td>
<td>Calibration of clipping height b</td></tr>
<tr><td width="25%"><a href="recFilter.html">rLSFilter</a></td>
Modified: pkg/robKalman/chm/0robKalman-package.html
===================================================================
--- pkg/robKalman/chm/0robKalman-package.html 2009-06-15 13:33:24 UTC (rev 33)
+++ pkg/robKalman/chm/0robKalman-package.html 2009-06-19 12:50:02 UTC (rev 34)
@@ -92,8 +92,11 @@
general recursive filters
+recursiveFilter
-KalmanFilter
- -rLSFilter
+ -rLSFilter:
+ *rLS.AO.Filter
+ *rLS.IO.Filter
-ACMfilter
+ -mACMfilter
ACMfilter:
+ACMfilt
@@ -181,6 +184,6 @@
-<hr><div align="center">[Package <em>robKalman</em> version 0.2.1 <a href="00Index.html">Index</a>]</div>
+<hr><div align="center">[Package <em>robKalman</em> version 0.3 <a href="00Index.html">Index</a>]</div>
</body></html>
Modified: pkg/robKalman/chm/calibrateRLS.html
===================================================================
--- pkg/robKalman/chm/calibrateRLS.html 2009-06-15 13:33:24 UTC (rev 33)
+++ pkg/robKalman/chm/calibrateRLS.html 2009-06-19 12:50:02 UTC (rev 34)
@@ -24,7 +24,7 @@
<h3>Usage</h3>
<pre>
-rLScalibrateB(Z, S, V, repl = 100000, eff, r, upto=20)
+rLScalibrateB(Z, S, V, repl = 100000, eff, r, upto=20, IO = FALSE)
</pre>
@@ -52,6 +52,9 @@
<tr valign="top"><td><code>upto</code></td>
<td>
an upper bound to <code>b</code> used in the zero-search of <code>uniroot</code> within <code>rLScalibrateB</code></td></tr>
+<tr valign="top"><td><code>IO</code></td>
+<td>
+logical of length 1: Is it rLS.IO (<code>TRUE</code>) or rLS[.AO] which is to be calibrated?</td></tr>
</table>
<h3>Details</h3>
@@ -118,6 +121,6 @@
-<hr><div align="center">[Package <em>robKalman</em> version 0.2.1 <a href="00Index.html">Index</a>]</div>
+<hr><div align="center">[Package <em>robKalman</em> version 0.3 <a href="00Index.html">Index</a>]</div>
</body></html>
Modified: pkg/robKalman/chm/recFilter.html
===================================================================
--- pkg/robKalman/chm/recFilter.html 2009-06-15 13:33:24 UTC (rev 33)
+++ pkg/robKalman/chm/recFilter.html 2009-06-19 12:50:02 UTC (rev 34)
@@ -7,6 +7,8 @@
<table width="100%"><tr><td>recFilter(robKalman)</td><td align="right">R Documentation</td></tr></table><object type="application/x-oleobject" classid="clsid:1e2a7bd0-dab9-11d0-b93a-00c04fc99f9e">
<param name="keyword" value="R: recFilter">
<param name="keyword" value="R: rLSFilter">
+<param name="keyword" value="R: rLS.AO.Filter">
+<param name="keyword" value="R: rLS.IO.Filter">
<param name="keyword" value="R: KalmanFilter">
<param name="keyword" value="R: ACMfilter">
<param name="keyword" value="R: recursiveFilter">
@@ -34,6 +36,8 @@
nsim = 0, seed = NULL, ..., dropRuns = TRUE)
KalmanFilter(Y, a, S, F, Q, Z, V, dropRuns = TRUE)
rLSFilter(Y, a, S, F, Q, Z, V, b, norm=Euclideannorm, dropRuns = TRUE)
+rLS.AO.Filter(Y, a, S, F, Q, Z, V, b, norm=Euclideannorm, dropRuns = TRUE)
+rLS.IO.Filter(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)
</pre>
@@ -265,11 +269,17 @@
<br>
<code>KalmanFilter(Y, a, S, F, Q, Z, V)</code> is a wrapper to <code>recursiveFilter(Y, a, S, F, Q, Z, V)</code>.<br>
-<code>rLSFilter(Y, a, S, F, Q, Z, V, b, norm=Euclideannorm)</code> is a wrapper to<br>
+<code>rLS.AO.Filter(Y, a, S, F, Q, Z, V, b, norm=Euclideannorm)</code> and (synonymously)
+<code>rLSFilter(Y, a, S, F, Q, Z, V, b, norm=Euclideannorm)</code> are wrappers to<br>
<code>recursiveFilter(Y, a, S, F, Q, Z, V,
initSc=.cKinitstep, predSc=.cKpredstep, corrSc=.cKcorrstep,
- initSr=.cKinitstep, predSr=.cKpredstep, corrSr=rLScorrstep,
+ initSr=.cKinitstep, predSr=.cKpredstep, corrSr=.rLScorrstep,
b=b, norm=norm)</code>.
+<code>rLS.IO.Filter(Y, a, S, F, Q, Z, V, b, norm=Euclideannorm)</code> is a wrapper to<br>
+<code>recursiveFilter(Y, a, S, F, Q, Z, V,
+ initSc=.cKinitstep, predSc=.cKpredstep, corrSc=.cKcorrstep,
+ initSr=.cKinitstep, predSr=.cKpredstep, corrSr=.rLS.IO.corrstep,
+ b=b, norm=norm)</code>.
<code>ACMFilter(Y, a, S, F, Q, Z, V, s0, psi, apsi, bpsi, cpsi, flag)</code> is a wrapper to <br>
<code>recursiveFilter(Y, a, S, F, Q, Z, V,
initSc=.cKinitstep, predSc=.cKpredstep, corrSc=.cKcorrstep,
@@ -332,13 +342,28 @@
# r = 0.1
(B2 <- rLScalibrateB(r = 0.1, S = SS, Z = Z0, V = V0i))
+# IO-filter
+# by efficiency in the ideal model
+# efficiency = 0.9
+(B1.IO <- rLScalibrateB(eff = 0.9, S = SS, Z = Z0, V = V0i, IO = TRUE))
+# by contamination radius
+# r = 0.1
+(B2.IO <- rLScalibrateB(r = 0.1, S = SS, Z = Z0, V = V0i, IO = TRUE))
+
erg <- KalmanFilter(Y, a = a0, S = SS0, F = F0, Q = Q0, Z = Z0, V = V0i)
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)
+
+rerg1.IO <- rLS.IO.Filter(Y, a = a0, S = SS0, F = F0, Q = Q0, Z = Z0, V = V0i, b = B1.IO$b)
+rerg2.IO <- rLS.IO.Filter(Y, a = a0, S = SS0, F = F0, Q = Q0, Z = Z0, V = V0i, b = B2.IO$b)
+
mean((X[,1,] - rerg1$Xf)^2) ### empirical MSE of the filters considered
mean((X[,1,] - rerg1$Xrf)^2)
mean((X[,1,] - rerg2$Xrf)^2)
+### not surprisingly IO-rob filter is not a good idea for AO - situation
+mean((X[,1,] - rerg1.IO$Xrf)^2)
+mean((X[,1,] - rerg2.IO$Xrf)^2)
</pre>
Modified: pkg/robKalman/chm/robKalman.chm
===================================================================
(Binary files differ)
Modified: pkg/robKalman/chm/robKalman.toc
===================================================================
--- pkg/robKalman/chm/robKalman.toc 2009-06-15 13:33:24 UTC (rev 33)
+++ pkg/robKalman/chm/robKalman.toc 2009-06-19 12:50:02 UTC (rev 34)
@@ -54,6 +54,14 @@
<param name="Local" value="recFilter.html">
</OBJECT>
<LI> <OBJECT type="text/sitemap">
+<param name="Name" value="rLS.AO.Filter">
+<param name="Local" value="recFilter.html">
+</OBJECT>
+<LI> <OBJECT type="text/sitemap">
+<param name="Name" value="rLS.IO.Filter">
+<param name="Local" value="recFilter.html">
+</OBJECT>
+<LI> <OBJECT type="text/sitemap">
<param name="Name" value="rLScalibrateB">
<param name="Local" value="calibrateRLS.html">
</OBJECT>
Modified: pkg/robKalman/man/0robKalman-package.Rd
===================================================================
--- pkg/robKalman/man/0robKalman-package.Rd 2009-06-15 13:33:24 UTC (rev 33)
+++ pkg/robKalman/man/0robKalman-package.Rd 2009-06-19 12:50:02 UTC (rev 34)
@@ -46,8 +46,11 @@
general recursive filters
+recursiveFilter
-KalmanFilter
- -rLSFilter
+ -rLSFilter:
+ *rLS.AO.Filter
+ *rLS.IO.Filter
-ACMfilter
+ -mACMfilter
ACMfilter:
+ACMfilt
Modified: pkg/robKalman/man/calibrateRLS.Rd
===================================================================
--- pkg/robKalman/man/calibrateRLS.Rd 2009-06-15 13:33:24 UTC (rev 33)
+++ pkg/robKalman/man/calibrateRLS.Rd 2009-06-19 12:50:02 UTC (rev 34)
@@ -7,7 +7,7 @@
calibrates the clipping height \code{b} of the rLS-filter in a time-invariant, linear, Gaussian state space model
}
\usage{
-rLScalibrateB(Z, S, V, repl = 100000, eff, r, upto=20)
+rLScalibrateB(Z, S, V, repl = 100000, eff, r, upto=20, IO = FALSE)
}
\arguments{
\item{Z}{observation matrix in the (ti-l-G-SSM); see below}
@@ -17,6 +17,7 @@
\item{eff}{efficiency w.r.t. classical Kalman filter in the ideal model}
\item{repl}{number of replicates used for a LLN-approximation of the expectations needed in this calibration}
\item{upto}{an upper bound to \code{b} used in the zero-search of \code{uniroot} within \code{rLScalibrateB}}
+ \item{IO}{logical of length 1: Is it rLS.IO (\code{TRUE}) or rLS[.AO] which is to be calibrated?}
}
\details{
We work in the setup of the time-invariant, linear, Gaussian state space model (ti-l-G-SSM)
Modified: pkg/robKalman/man/recFilter.Rd
===================================================================
--- pkg/robKalman/man/recFilter.Rd 2009-06-15 13:33:24 UTC (rev 33)
+++ pkg/robKalman/man/recFilter.Rd 2009-06-19 12:50:02 UTC (rev 34)
@@ -1,6 +1,8 @@
\name{recFilter}
\alias{recFilter}
\alias{rLSFilter}
+\alias{rLS.AO.Filter}
+\alias{rLS.IO.Filter}
\alias{KalmanFilter}
\alias{ACMfilter}
\alias{recursiveFilter}
@@ -19,6 +21,8 @@
nsim = 0, seed = NULL, ..., dropRuns = TRUE)
KalmanFilter(Y, a, S, F, Q, Z, V, dropRuns = TRUE)
rLSFilter(Y, a, S, F, Q, Z, V, b, norm=Euclideannorm, dropRuns = TRUE)
+rLS.AO.Filter(Y, a, S, F, Q, Z, V, b, norm=Euclideannorm, dropRuns = TRUE)
+rLS.IO.Filter(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)
}
@@ -160,11 +164,17 @@
e.g. in case of the ACM filter each element contains a corresponding value of \code{st}}
\code{KalmanFilter(Y, a, S, F, Q, Z, V)} is a wrapper to \code{recursiveFilter(Y, a, S, F, Q, Z, V)}.\cr
-\code{rLSFilter(Y, a, S, F, Q, Z, V, b, norm=Euclideannorm)} is a wrapper to\cr
+\code{rLS.AO.Filter(Y, a, S, F, Q, Z, V, b, norm=Euclideannorm)} and (synonymously)
+\code{rLSFilter(Y, a, S, F, Q, Z, V, b, norm=Euclideannorm)} are wrappers to\cr
\code{recursiveFilter(Y, a, S, F, Q, Z, V,
initSc=.cKinitstep, predSc=.cKpredstep, corrSc=.cKcorrstep,
- initSr=.cKinitstep, predSr=.cKpredstep, corrSr=rLScorrstep,
+ initSr=.cKinitstep, predSr=.cKpredstep, corrSr=.rLScorrstep,
b=b, norm=norm)}.
+\code{rLS.IO.Filter(Y, a, S, F, Q, Z, V, b, norm=Euclideannorm)} is a wrapper to\cr
+\code{recursiveFilter(Y, a, S, F, Q, Z, V,
+ initSc=.cKinitstep, predSc=.cKpredstep, corrSc=.cKcorrstep,
+ initSr=.cKinitstep, predSr=.cKpredstep, corrSr=.rLS.IO.corrstep,
+ b=b, norm=norm)}.
\code{ACMFilter(Y, a, S, F, Q, Z, V, s0, psi, apsi, bpsi, cpsi, flag)} is a wrapper to \cr
\code{recursiveFilter(Y, a, S, F, Q, Z, V,
initSc=.cKinitstep, predSc=.cKpredstep, corrSc=.cKcorrstep,
@@ -211,14 +221,28 @@
# r = 0.1
(B2 <- rLScalibrateB(r = 0.1, S = SS, Z = Z0, V = V0i))
+# IO-filter
+# by efficiency in the ideal model
+# efficiency = 0.9
+(B1.IO <- rLScalibrateB(eff = 0.9, S = SS, Z = Z0, V = V0i, IO = TRUE))
+# by contamination radius
+# r = 0.1
+(B2.IO <- rLScalibrateB(r = 0.1, S = SS, Z = Z0, V = V0i, IO = TRUE))
erg <- KalmanFilter(Y, a = a0, S = SS0, F = F0, Q = Q0, Z = Z0, V = V0i)
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)
+
+rerg1.IO <- rLS.IO.Filter(Y, a = a0, S = SS0, F = F0, Q = Q0, Z = Z0, V = V0i, b = B1.IO$b)
+rerg2.IO <- rLS.IO.Filter(Y, a = a0, S = SS0, F = F0, Q = Q0, Z = Z0, V = V0i, b = B2.IO$b)
+
mean((X[,1,] - rerg1$Xf)^2) ### empirical MSE of the filters considered
mean((X[,1,] - rerg1$Xrf)^2)
mean((X[,1,] - rerg2$Xrf)^2)
+### not surprisingly IO-rob filter is not a good idea for AO - situation
+mean((X[,1,] - rerg1.IO$Xrf)^2)
+mean((X[,1,] - rerg2.IO$Xrf)^2)
}
More information about the Robkalman-commits
mailing list