[Robkalman-commits] r26 - branches/robKalman_itwm/pkg/robKalman branches/robKalman_itwm/pkg/robKalman/R branches/robKalman_itwm/pkg/robKalman/man branches/robkalman_pr/pkg/robKalman 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
Wed Mar 18 18:40:44 CET 2009
Author: ruckdeschel
Date: 2009-03-18 18:40:44 +0100 (Wed, 18 Mar 2009)
New Revision: 26
Removed:
pkg/robKalman/man/0robKalman-package.Rd
pkg/robKalman/man/ACMfilt.Rd
pkg/robKalman/man/arGM.Rd
pkg/robKalman/man/calibrateRLS.Rd
pkg/robKalman/man/internalACM.Rd
pkg/robKalman/man/internalKalman.Rd
pkg/robKalman/man/internalarGM.Rd
pkg/robKalman/man/internalpsi.Rd
pkg/robKalman/man/internalrLS.Rd
pkg/robKalman/man/recFilter.Rd
pkg/robKalman/man/simulateSScont.Rd
pkg/robKalman/man/util.Rd
Modified:
branches/robKalman_itwm/pkg/robKalman/DESCRIPTION
branches/robKalman_itwm/pkg/robKalman/R/0AllClass.R
branches/robKalman_itwm/pkg/robKalman/R/ACMfilter.R
branches/robKalman_itwm/pkg/robKalman/R/Psi.R
branches/robKalman_itwm/pkg/robKalman/R/Util.R
branches/robKalman_itwm/pkg/robKalman/R/calibrateRLS.R
branches/robKalman_itwm/pkg/robKalman/R/classKalman.R
branches/robKalman_itwm/pkg/robKalman/R/recFilter.R
branches/robKalman_itwm/pkg/robKalman/R/simulateSScont.R
branches/robKalman_itwm/pkg/robKalman/man/0robKalman-package.Rd
branches/robKalman_itwm/pkg/robKalman/man/internalrLS.Rd
branches/robKalman_itwm/pkg/robKalman/man/recFilter.Rd
branches/robKalman_itwm/pkg/robKalman/man/util.Rd
branches/robkalman_pr/pkg/robKalman/DESCRIPTION
branches/robkalman_pr/pkg/robKalman/man/0robKalman-package.Rd
branches/robkalman_pr/pkg/robKalman/man/internalrLS.Rd
branches/robkalman_pr/pkg/robKalman/man/recFilter.Rd
branches/robkalman_pr/pkg/robKalman/man/util.Rd
pkg/robKalman/DESCRIPTION
pkg/robKalman/NAMESPACE
pkg/robKalman/R/0AllClass.R
pkg/robKalman/R/ACMfilt.R
pkg/robKalman/R/ACMfilter.R
pkg/robKalman/R/Psi.R
pkg/robKalman/R/Util.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/0robKalman-package.html
pkg/robKalman/chm/ACMfilt.html
pkg/robKalman/chm/arGM.html
pkg/robKalman/chm/calibrateRLS.html
pkg/robKalman/chm/internalACM.html
pkg/robKalman/chm/internalKalman.html
pkg/robKalman/chm/internalarGM.html
pkg/robKalman/chm/internalpsi.html
pkg/robKalman/chm/internalrLS.html
pkg/robKalman/chm/recFilter.html
pkg/robKalman/chm/robKalman.chm
pkg/robKalman/chm/robKalman.hhp
pkg/robKalman/chm/robKalman.toc
pkg/robKalman/chm/simulateSScont.html
pkg/robKalman/chm/util.html
Log:
reorganized trunc and branches; package pkg/robKalman now checks without error.
branch branches/robKalman_itwm starts with this package (and an updated DESCRIPTION file)
branches branches/robKalman_pr, branches/robKalman_bs receive updated DESCRIPTION files
but else only receive the package files of the former revisions sorted under
branches/robKalman_<...>/pkg/robKalman...
Modified: branches/robKalman_itwm/pkg/robKalman/DESCRIPTION
===================================================================
--- branches/robKalman_itwm/pkg/robKalman/DESCRIPTION 2009-03-18 15:41:02 UTC (rev 25)
+++ branches/robKalman_itwm/pkg/robKalman/DESCRIPTION 2009-03-18 17:40:44 UTC (rev 26)
@@ -1,13 +1,12 @@
Package: robKalman
-Version: 0.2
-Date: 2008-07-28
+Version: 0.3
+Date: 2009-03-18
Title: Robust Kalman Filtering
Description: Routines for Robust Kalman Filtering --- the ACM- and rLS-filter
-Author: Peter Ruckdeschel, Bernhard Spangl
+Author: Peter Ruckdeschel, Bernhard Spangl, Irina Ursachi, Cezar Chirila
Maintainer: Peter Ruckdeschel <Peter.Ruckdeschel at itwm.fraunhofer.de>
Depends: R(>= 2.3.0), methods, graphics, startupmsg, dse1, dse2, MASS, limma, robustbase, numDeriv
Imports: stats, MASS
-SaveImage: no
LazyLoad: yes
-License: GPL (version 2 or later)
-URL: https://r-forge.r-project.org/projects/robkalman
+License: LGPL-3
+URL: http://robkalman.r-forge.r-project.org/
Modified: branches/robKalman_itwm/pkg/robKalman/R/0AllClass.R
===================================================================
--- branches/robKalman_itwm/pkg/robKalman/R/0AllClass.R 2009-03-18 15:41:02 UTC (rev 25)
+++ branches/robKalman_itwm/pkg/robKalman/R/0AllClass.R 2009-03-18 17:40:44 UTC (rev 26)
@@ -1,5 +1,3 @@
-################ general package preparation code, nothing to do with classes:
-
.onLoad <- function(lib, pkg){
require("methods", character = TRUE, quietly = TRUE)
@@ -14,322 +12,3 @@
invisible()
}
-#.onUnload <- function(libpath)
-#{
-# library.dynam.unload("distrEx", libpath)
-#}
-#
-
-
-################ Matrix class
-
-#### Code borrowed from package distrMod
-
-if (!(isClass("PosSemDefSymmMatrix")))
-setClass("PosSemDefSymmMatrix", contains = "matrix",
- prototype = prototype(matrix(1)),
- validity = function(object){
- if(nrow(object) != ncol(object))
- stop("no square matrix")
- if(any(!is.finite(object)))
- stop("inifinite or missing values in matrix")
- if(!isTRUE(all.equal(object, t(object), .Machine$double.eps^0.5)))
- stop("matrix is not symmetric")
- if(!all(eigen(object)$values > -100*.Machine$double.eps))
- stop("matrix is (numerically) not positive semi - definite")
- return(TRUE)
- })
-
-## positive definite, symmetric matrices with finite entries
-if (!(isClass("PosDefSymmMatrix")))
-setClass("PosDefSymmMatrix", contains = "PosSemDefSymmMatrix",
- validity = function(object){
- if(!all(eigen(object)$values > 100*.Machine$double.eps))
- stop("matrix is (numerically) not positive definite")
- valid <- getValidity(getClass("PosSemDefSymmMatrix"))
- valid(as(object, "PosSemDefSymmMatrix"))
- return(TRUE)
- })
-
-
-
-#
-## register zoo as "S4"-class
-setOldClass("zoo")
-#
-
-### infra-structure classes / class unions
-
-setClassUnion("IntegerOrNULL", c("integer", "NULL"))
-setClassUnion("ArrayOrNULL", c("array", "NULL"))
-setClassUnion("ArrayOrMatrix", c("array", "matrix"))
-setClassUnion("Hyperparamtype",
- c("NULL","ArrayOrMatrix", "OptionalFunction"))
-setClassUnion("sHyperparamtype",
- c("Hyperparamtype", "numeric"))
-setClassUnion("MatrixOrLogical", c("logical", "matrix"))
-
-
-###############################################################################
-#
-# State Space Model classes (SSMs)
-#
-###############################################################################
-
-
-# class SSM --- State space model
-setClass("SSM",
- representation = representation(
- name = "character", ## name of the ssm
- F = "Hyperparamtype", ## transition matrix/ces or NULL
- Z = "Hyperparamtype", ## observation matrix/ces or NULL
- Q = "Hyperparamtype", ## innovation covariance or NULL
- V = "Hyperparamtype", ## observation error covariance or NULL
- p = "numeric", ## state dimension
- q = "numeric", ## observation dimension
- a = "numeric", ## mean value of starting state
- S = "Hyperparamtype", ## variance of starting state
- time = "zoo"), ## time index
- prototype = prototype(name = gettext("a state space"),
- F = NULL,
- Z = NULL,
- Q = NULL,
- V = NULL,
- p = 1,
- q = 1,
- a = 0,
- S = NULL,
- time = zoo(1)),
- )
-
-# class TimeInvariantSSM
-setClass("TimeInvariantSSM",
- prototype = prototype(name = gettext("a time-invariant state space"),
- F = matrix(1),
- Z = matrix(1),
- Q = matrix(1),
- V = matrix(1),
- p = 1,
- q = 1,
- a = 0,
- S = matrix(1)),
- validity = function(object){
- if(!is.matrix(object at F)|!is.matrix(object at Z)|
- !is.matrix(object at Q)|!is.matrix(object at V)|
- !is.matrix(object at S))
- stop("Hyperparameters have to be matrices")
- return(TRUE)
- },
- contains = "SSM")
-
-###############################################################################
-#
-# Filter classes
-#
-###############################################################################
-
-setClass("recFilter", representation(name = "character",
- SSM = "SSM",
- Y = "array",
- X.filtered = "array",
- X.predicted = "array",
- Cov.filtered = "array",
- Cov.predicted = "array",
- Kalman.Gain = "array",
- time = "zoo"),
- prototype = prototype(name="classical Kalman Filter",
- SSM = new("TimeInvariantSSM"),
- Y = array(1,dim=c(1,1,1)),
- X.filtered = array(1,dim=c(1,1,1)),
- X.predicted = array(1,dim=c(1,1,1)),
- Cov.filtered = array(1,dim = c(1,1,1)),
- Cov.predicted = array(1,dim = c(1,1,1)),
- Kalman.Gain = array(1,dim = c(1,1,1)),
- time = zoo(1))
-
- )
-
-
-setClass("robrecFilter", representation(
- name.rob = "character",
- X.rob.filtered = "array",
- X.rob.predicted = "array",
- Cov.rob.filtered = "array",
- Cov.rob.predicted = "array",
- Kalman.rob.Gain = "array",
- IndIO = "MatrixOrLogical",
- IndAO = "MatrixOrLogical",
- rob.correction.ctrl = "list",
- rob.prediction.ctrl = "list",
- nsim = "numeric",
- RNGstate = "IntegerOrNULL",
- Cov.rob.filtered.sim = "ArrayOrNULL",
- Cov.rob.predicted.sim = "ArrayOrNULL"
- ),
- prototype = prototype(
- name="rLS Filter",
- X.rob.filtered = array(1,dim=c(1,1,1)),
- X.rob.predicted = array(1,dim=c(1,1,1)),
- Cov.rob.filtered = array(1,dim = c(1,1,1)),
- Cov.rob.predicted = array(1,dim = c(1,1,1)),
- Kalman.rob.Gain = array(1,dim = c(1,1,1)),
- IndIO = FALSE,
- IndAO = FALSE,
- nsim = 0,
- RNGstate = as.integer(0),
- rob.correction.ctrl = list(NULL),
- rob.prediction.ctrl = list(NULL),
- Cov.rob.filtered.sim = array(1,dim = c(1,1,1)),
- Cov.rob.predicted.sim = array(1,dim = c(1,1,1))),
- contains = "recFilter")
-
-###############################################################################
-#
-# Control classes
-#
-###############################################################################
-#setClass("ACMcontrol",representation())
-#setClass("rLScontrol",representation())
-
-
-setClass("RecFiltControl", representation(
- name = "character",
- init = "function",
- predict = "function",
- correct = "function"),
- contains = "VIRTUAL")
-
-setClass("KalmanControl",
- prototype = prototype(
- name = "classical Kalman Filter",
- init = .cKinitstep,
- predict = .cKpredstep,
- correct = .cKcorrstep),
- contains="RecFiltControl"
- )
-
-setClass("robrecControl", representation(
- controls = "list",
- name.rob = "character",
- init.rob = "function",
- predict.rob = "function",
- correct.rob = "function"),
- prototype = prototype(
- name.rob = "rLS Filter",
- init.rob = .cKinitstep,
- predict.rob = .cKpredstep,
- correct.rob = .rLScorrstep,
- controls = list(b = 2, norm = EuclideanNorm)
- ),
- contains="RecFiltControl"
- )
-# = paste(gettext(
-# "Control set and init, prediction and"
-# ),gettext(
-# "correction step for the classical Kalman Filter\n"
-# ),gettext(
-# "and the rLS Filter"
-# )),
-
-###############################################################################
-#
-# multivariate Distribution classes
-#
-###############################################################################
-
-setClass("SSMDistribution.f", representation(
- r.init = "function",
- r.innov = "function",
- r.obs = "function"),
- prototype = prototype(r.init = mvrnorm,
- r.innov = mvrnorm,
- r.obs = mvrnorm))
-
-setClass("SSMellDistribution.f", representation(
- m.init = "sHyperparamtype",
- S.init = "Hyperparamtype",
- m.innov = "sHyperparamtype",
- S.innov = "Hyperparamtype",
- m.obs = "sHyperparamtype",
- S.obs = "Hyperparamtype"),
- prototype = prototype(m.init = 0, S.init = matrix(1),
- m.innov = 0, S.innov = matrix(1),
- m.obs = 0, S.obs = matrix(1)),
- contains = "SSMDistribution.f")
-
-setClass("SSMConvDistribution.f", representation(
- ideal = "SSMDistribution.f",
- cont = "SSMellDistribution.f",
- r.IO = "numeric",
- r.AO = "numeric"),
- prototype = prototype(ideal = new("SSMellDistribution.f"),
- cont = new("SSMellDistribution.f"),
- r.IO = 0, r.AO = 0),
- validity = function(object){
- if (object at r.IO<0||object at r.AO<0||object at r.IO>1||object at r.AO>1)
- stop("Radii must be between 0 and 1")
- return(TRUE)})
-
-
-###############################################################################
-#
-# SSM + Distribution classes
-#
-###############################################################################
-setClass("SSMwithDistribution", representation(
- SSM = "SSM",
- Distribution = "SSMDistribution.f"),
- prototype = prototype(SSM = new("SSM"),
- Distribution = SSMDistribution.f(new("SSM"))))
-setClass("SSMwithConvDistribution", representation(
- SSM = "SSM",
- Distribution = "SSMConvDistribution.f"),
- prototype = prototype(SSM = new("SSM"),
- Distribution = SSMContDistribution.f(new("SSM"))))
-
-setClassUnion("SSMDistr", c("SSMDistribution.f",
- "SSMConvDistribution.f"))
-
-###############################################################################
-#
-# SSM - Simulation classes
-#
-###############################################################################
-setClass("SSMsimulation", representation(
- SSM = "SSM",
- Distr = "SSMDistr",
- RNGstate = "numeric",
- states = "ArrayOrMatrix",
- obs = "ArrayOrMatrix"),
- prototype = prototype(SSM = new("SSM"),
- Distr = new("SSMDistribution.f"),
- RNGstate = structure(1, kind = as.list(RNGkind())),
- states = matrix(1), obs = matrix(1)))
-
-setClass("SSMcontSimulation", representation(
- states.id = "ArrayOrMatrix",
- obs.id = "ArrayOrMatrix",
- Ind.IO = "MatrixOrLogical",
- Ind.AO = "MatrixOrLogical"
- ),
- prototype = prototype(
- Distr = new("SSMConvDistribution.f"),
- states.id = matrix(1), obs.id = matrix(1),
- Ind.IO = FALSE, Ind.AO = FALSE
- ),
- validity = function(object){
- fct <- function(m){
- mt <- paste(deparse(substitute(m)),sep="",collapse="")
- if(is.matrix(m)){
- if(!all(is.logical(m)))
- stop(gettextf("Matrix %s has to have logical entries.", mt))
- return(TRUE)
- }
- else return(TRUE)
- }
- return(fct(object at Ind.IO)&&fct(object at Ind.AO))
- },
- contains = "SSMsimulation")
-
-
-
\ No newline at end of file
Modified: branches/robKalman_itwm/pkg/robKalman/R/ACMfilter.R
===================================================================
--- branches/robKalman_itwm/pkg/robKalman/R/ACMfilter.R 2009-03-18 15:41:02 UTC (rev 25)
+++ branches/robKalman_itwm/pkg/robKalman/R/ACMfilter.R 2009-03-18 17:40:44 UTC (rev 26)
@@ -19,16 +19,16 @@
}
##steps for classical Kalman filter (cK)
-.ACMinitstep <- function(a, S, i, ...)
+.ACMinitstep <- function(a, S, ...)
{dots <- list(...)
if(hasArg("s0"))
- s0 <- dots$"s0"
+ s0<-dots$"s0"
else
- s0 <- NULL
+ s0<-NULL
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, rob0, s0, ...) ### S=P F= Phi
{
###########################################
##
@@ -45,11 +45,10 @@
## 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 = FALSE))
+ return(list(x1 = F %*% x0, S1 = S1, rob1 = sqrt(S1[1, 1] + s0), Ind=1))
}
-.ACMcorrstep <- function (y, x1, S1, Z, V, i, rob1, dum = NULL,
- psi, apsi, bpsi, cpsi, flag, ...)
+.ACMcorrstep <- function (y, x1, S1, Z, V, rob1, dum=NULL, psi, apsi, bpsi, cpsi, flag, ...)
{
###########################################
##
@@ -68,15 +67,11 @@
## dum ... dummy variable for compatibility with ... argument of calling function
## V ... covariance matrix of observation noise
## psi ... influence function to be used
-## apsi, bpsi, cpsi ... tuning constants for Hampel's psi-function
-## (default: apsi=bpsi=2.5, cpsi=5.0)
+## a, b, c ... tuning constants for Hampel's psi-function
+## (defaul: a=b=2.5, c=5.0)
## flag ... character, if "weights" (default), use psi(t)/t to calculate
## the weights; if "deriv", use psi'(t)
-
- # to be compatible with parallel computing of a bunch of time series
- y <- y[,1]
- x1 <- x1[,1]
-
+
st <- rob1
K <- .getKG(S1, Z, V)
@@ -93,5 +88,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, Ind=ind, rob0=rob1))
}
Modified: branches/robKalman_itwm/pkg/robKalman/R/Psi.R
===================================================================
--- branches/robKalman_itwm/pkg/robKalman/R/Psi.R 2009-03-18 15:41:02 UTC (rev 25)
+++ branches/robKalman_itwm/pkg/robKalman/R/Psi.R 2009-03-18 17:40:44 UTC (rev 26)
@@ -142,7 +142,7 @@
## Paramters:
## x ... vector
## a, b, c ... tuning constants
- x.norm <- EuclidianNorm(x)
+ x.norm <- Euclideannorm(x)
small <- (x.norm <= a)
dummy <- pmin(a, a/(c-b)*(c-x.norm))
dummy <- pmax(dummy, 0)/(x.norm+small)*(!small) + small
Modified: branches/robKalman_itwm/pkg/robKalman/R/Util.R
===================================================================
--- branches/robKalman_itwm/pkg/robKalman/R/Util.R 2009-03-18 15:41:02 UTC (rev 25)
+++ branches/robKalman_itwm/pkg/robKalman/R/Util.R 2009-03-18 17:40:44 UTC (rev 26)
@@ -1,16 +1,13 @@
### the Euclidean norm
-EuclideanNorm <- function(x) {sqrt(colSums(x^2))}
+Euclideannorm <- function(x) {sqrt(sum(x^2))}
### huberizing a vector to length b
-Huberize <- function(x, b, norm = EuclideanNorm, ...)
-{ nx <- norm(x, ...)
- Ind0 <- (nx < b )
- Ind1 <- (nx < b/2)
- x*(Ind0 + (1-Ind0)* b / (nx + Ind1) )
-}
+Huberize <- function(x, b, norm=Euclideannorm, ...)
+ x*ifelse(norm(x) < b, 1, b/norm(x, ...))
+
limitS <- function(S, F, Z, Q, V, tol = 10^-4, itmax = 1000)#
## determines lim_{t->infty} S_{t|t-1}
Modified: branches/robKalman_itwm/pkg/robKalman/R/calibrateRLS.R
===================================================================
--- branches/robKalman_itwm/pkg/robKalman/R/calibrateRLS.R 2009-03-18 15:41:02 UTC (rev 25)
+++ branches/robKalman_itwm/pkg/robKalman/R/calibrateRLS.R 2009-03-18 17:40:44 UTC (rev 26)
@@ -27,20 +27,27 @@
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)
+ {f <- function(b, dX = dx, dX0 = dx0, no0 = no, r0=r,
+ eff0 = eff, trS0 = trS, 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)
+ 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)
b <- erg$root
if (missing(r)) ### corresponding radius is calculated
Modified: branches/robKalman_itwm/pkg/robKalman/R/classKalman.R
===================================================================
--- branches/robKalman_itwm/pkg/robKalman/R/classKalman.R 2009-03-18 15:41:02 UTC (rev 25)
+++ branches/robKalman_itwm/pkg/robKalman/R/classKalman.R 2009-03-18 17:40:44 UTC (rev 26)
@@ -51,10 +51,10 @@
{list( x0 = a, S0 = S )}
.cKpredstep <- function(x0, S0, F, Q, ...)
- {list( x1 = F %*% x0, S1 = .getpredCov(S0, F, Q), Ind = FALSE)}
+ {list( x1 = F %*% x0, S1 = .getpredCov(S0, F, Q), Ind=1)}
.cKcorrstep <- function(y, x1, S1, Z, V, ...)
{K <- .getKG(S1, Z, V)
x0 <- x1 + K %*% (y - Z %*% x1)
S0 <- .getcorrCov(S1, K, Z)
- list(x0 = x0, K = K, S0 = S0, Ind = FALSE)}
+ list(x0 = x0, K = K, S0 = S0, Ind=1)}
Modified: branches/robKalman_itwm/pkg/robKalman/R/recFilter.R
===================================================================
--- branches/robKalman_itwm/pkg/robKalman/R/recFilter.R 2009-03-18 15:41:02 UTC (rev 25)
+++ branches/robKalman_itwm/pkg/robKalman/R/recFilter.R 2009-03-18 17:40:44 UTC (rev 26)
@@ -1,143 +1,61 @@
-setMethod("kalmanRob", signature(method = "robrecControl", smooth = "missing"),
- function(method = rLSControl(), Y, SSM, nsim = 0, seed = NULL){
- SSMa <- makeArrayRepresentation(SSM)
- erg <- do.call(recursiveFilter,
- args = c(list(Y = Y, a = SSM at a, S = SSMa at S,
- F = SSMa at F, Q = SSMa at Q, Z = SSMa at Z, V = SSMa at V,
- nsim = nsim, seed = seed,
- initSc = init(method),
- predSc = predict(method),
- corrSc = correct(method),
- initSr = init.rob(method),
- predSr = predict.rob(method),
- corrSr = correct.rob(method)),
- controls(method), dropRuns = FALSE)
- )
-
- return(generateRobRecFilter(
- name = name(method), name.rob = name.rob(method),
- SSM = SSM, Y = Y, time = SSM at time,
- Xf = erg$Xf, Xp = erg$Xp, S0 = erg$S0, S1 = erg$S1, KG = erg$KG,
- Xrf = erg$Xrf, Xrp = erg$Xrp, Sr0 = erg$Sr0, Sr1 = erg$Sr1,
- KGr = erg$KGr, IndIO = erg$IndIO, IndAO = erg$IndAO,
- rob0L = erg$rob0L, rob1L = erg$rob1L, St0s = erg$St0s,
- St1s = erg$St1s, nsim = erg$nsim, RNGstate = erg$RNGstate
- ))
- })
+#kalmanRob <- function(Y, SSM, method=c(myRLS,myACM) functionDIEAngabegeneriert, control())
+# {methodRes <- method() # liste ()
-setMethod("kalman", signature(smooth = "missing"),
- function(Y, SSM){
- method <- KalmanControl()
- SSMa <- makeArrayRepresentation(SSM)
- erg <- recursiveFilter(Y = Y, a = SSM at a, S = SSMa at S,
- F = SSMa at F, Q = SSMa at Q, Z = SSMa at Z, V = SSMa at V,
- initSc = init(method),
- predSc = predict(method),
- corrSc = correct(method), dropRuns = FALSE)
- return(generateRecFilter(
- name = name(method), SSM = SSM, Y = Y, time = SSM at time,
- Xf = erg$Xf, Xp = erg$Xp, S0 = erg$S0, S1 = erg$S1, KG = erg$KG
- ))})
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
+ initSc=.cKinitstep, predSc=.cKpredstep, corrSc=.cKcorrstep,
+ initSr=NULL, predSr=NULL, corrSr=NULL, ...)# a generalization of the Kalmanfilter
#arguments:
-# + Y :observations in an array with dimensions qd x runs x tt
-# :for backward compatibility: or a vector /a matrix in
-# dimensions qd x tt
+# + Y :observations
# + a, S, F, Q, Z, V: Hyper-parameters of the ssm
# + initSc, predSc, corrSc: (classical) initialization-, prediction-, and correction-step function
# + initSr, predSr, corrSr: (robust) initialization-, prediction-, and correction-step function
# + robustIO: if TRUE indicators are recorded whether prediction step does clipping
# + robustAO: if TRUE indicators are recorded whether correction step does clipping
-# + nsim: if >0 we simulate a bunch of nsim paths (acc. to ideal model) to get emp. covariances
-# + seed: seed for the simulations
+# if (robustIO|)
# + ... additional arguments for initS, predS, corrS
-# + 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)){
- Y0 <- Y
- tt <- ifelse(length(Z)==1, length(Y), (dim(Y))[2])
- Y <- aperm(array(Y, dim = c(qd,tt,1)),c(1,3,2))
- }
- ########################
-
- tt <- (dim(Y))[3]
- runs <- (dim(Y))[2]
-
+{qd <- ifelse(length(Z)==1, 1, (dim(Y))[1])
pd <- length(a)
+ tt <- ifelse(length(Z)==1, length(Y), (dim(Y))[2])
+# browser()
IndIO <- NULL
IndAO <- 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))
- if(is.matrix(Q)) Q <- array(Q, dim = c(pd,pd,tt))
- 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))
- 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))
+ Xf <- matrix(0, length(a), tt + 1)
+ Xp <- matrix(0, length(a), tt)
+ St0 <- array(0, c(pd, pd, tt + 1))
+ St1 <- array(0, c(pd, pd, tt))
+ KG <- array(0, c(pd, qd, tt))
- if (nsim)
- {if (!exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE))
- runif(1)
- if (is.null(seed))
- RNGstate <- get(".Random.seed", envir = .GlobalEnv)
- else {
- R.seed <- get(".Random.seed", envir = .GlobalEnv)
- set.seed(seed)
- RNGstate <- structure(seed, kind = as.list(RNGkind()))
- on.exit(assign(".Random.seed", R.seed, envir = .GlobalEnv))
- }
- }else{
- RNGstate <- NULL
- }
-
+
if(robust)
- {Xrf <- array(0, dim = c(pd, runs, tt + 1))
- Xrp <- array(0, dim = c(pd, runs, tt))
- 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))
+ {Xrf <- matrix(0, length(a), tt + 1)
+ Xrp <- matrix(0, length(a), tt)
+ Str0 <- array(0, c(pd, pd, tt + 1))
+ Str1 <- array(0, c(pd, pd, tt))
+ KGr <- array(0, c(pd, qd, tt))
}
if(!is.null(predSr))
- IndIO <- matrix(FALSE,runs,tt)
+ IndIO <- numeric(tt)
if(!is.null(corrSr))
- IndAO <- matrix(FALSE,runs,tt)
+ IndAO <- numeric(tt)
- ini <- initSc(a, S, i = 0, ...)
+ ini <- initSc(a, S, ...)
x0 <- ini$x0
S0 <- ini$S0
- Xf[, , 1] <- ini$x0
+ 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))
- }
- if(!is.null(initSr))
- {inir <- initSr(a, S, i = 0, ...)
+ {if(!is.null(initSr))
+ {inir <- initSr(a, S, ...)
xr0 <- inir$x0
Sr0 <- inir$S0
rob0 <- inir$rob
@@ -148,7 +66,7 @@
Sr0 <- S0
rob0 <- NULL
}
- Xrf[,, 1] <- xr0
+ Xrf[, 1] <- xr0
Str0[,, 1] <- xr0
rob0L <- list(rob0)
@@ -163,43 +81,27 @@
}
for (i in (1:tt))
{#prediction
- F0 <- matrix(F[,,i], nrow = pd, ncol = pd)
- Q0 <- matrix(Q[,,i], nrow = pd, ncol = pd)
-
- ps <- predSc(x0 = x0, S0 = S0, F = F0, Q = Q0, i = i, ...)
+ ps <- predSc(x0=x0, S0=S0, F=F, Q=Q, ...)
x1 <- ps$x1
S1 <- ps$S1
- Xp[,, i] <- x1
+ Xp[, i] <- x1
St1[,, i] <- S1
+
if(robust)
{if(!is.null(predSr))
- {psr <- predSr(x0 = xr0, S0 = Sr0, F = F0, Q = Q0, i = i,
- ..., rob0 = rob0)
- IndIO[,i] <- as.logical(psr$Ind)
- if(nsim){
- vs <- t(rmvnorm(nsim, a*0, Q0))
- Xs <- F0 %*% Xs + vs
- xr1s <- predSr(x0 = xr0s, S0 = Sr0, F = F0, Q = Q0, i = i,
- ..., rob0 = rob0)$x1
- St1s[,,i] <- cov(t(xr1s))
- }
- }else{
- psr <- predSc(x0 = xr0, S0 = Sr0, F = F0, Q = Q0, i = i, ...)
- if(nsim){
- vs <- t(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))
- }
- }
+ {psr <- predSr(x0 = xr0, S0 = Sr0, F = F, Q = Q, ..., rob0 = rob0)
+ IndIO[i] <- psr$Ind}
+ else
+ psr <- predSc(x0 = xr0, S0 = Sr0, F = F, Q = Q, ...)
+
+
xr1 <- psr$x1
Sr1 <- psr$S1
rob1 <- psr$rob1
- Xrp[,, i] <- xr1
+ Xrp[, i] <- xr1
Str1[,, i]<- S1
if(i==1) rob1L <- list(rob1)
else rob1L[[i]] <- rob1
@@ -207,68 +109,39 @@
#correction
- Z0 <- matrix(Z[,,i], nrow = qd, ncol = pd)
- V0 <- matrix(V[,,i], nrow = qd, ncol = qd)
-
- Y0 <- matrix(Y[,,i], nrow = qd, ncol = runs)
- cs <- corrSc(y = Y0, x1 = x1, S1 = S1, Z = Z0, V = V0, i = i,...)
+ cs <- corrSc(y = Y[, i], x1 = x1, S1 = S1, Z = Z, V = V, ...)
x0 <- cs$x0
S0 <- cs$S0
- Xf[,, i + 1] <- x0
+ Xf[, i + 1] <- x0
St0[,, i + 1] <- S0
KG[,, i] <- cs$K
if(robust)
{if (!is.null(corrSr))
- {csr <- corrSr(y = Y0, x1 = xr1, S1 = Sr1,
- Z = Z0, V = V0, i = i, ..., rob1 = rob1)
- 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,
- Z = Z0, V = V0, i = i, ..., rob1 = rob1)$x0
- St0s[,,i] <- cov(t(xr0s))
- }
- }else{
- csr <- corrSc(y = Y0, x1 = xr1, S1 = Sr1, Z = Z0, V = V0,
- i = i, ...)
- if(nsim){
- es <- t(rmvnorm(nsim, Y[,1]*0, V0))
- Ys <- Z0 %*% Xs + es
- xr0s <- corrSc(y = Ys, x1 = xr1s, S1 = Sr1,
- Z = Z0, V = V0, i = i, ...)$x0
- St0s[,,i] <- cov(t(xr0s))
- }
- }
+ {csr <- corrSr(y = Y[, i], x1 = xr1, S1 = Sr1, Z = Z, V = V, ..., rob1 = rob1)
+ IndAO[i] <- csr$Ind }
+ else
+ csr <- corrSc(y = Y[, i], x1 = xr1, S1 = Sr1, Z = Z, V = V, ...)
+
xr0 <- csr$x0
Sr0 <- csr$S0
rob0 <- csr$rob0
- Xrf[,, i + 1] <- xr0
+ Xrf[, i + 1] <- xr0
Str0[,, i + 1] <- S0
rob0L[[i + 1]] <- rob0
KGr[,, i] <- csr$K
}
}
-if((runs==1)&&(dropRuns))
- {Xf <- matrix(Xf,pd,tt+1)
- Xp <- matrix(Xp,pd,tt)
- if(!is.null(Xrp)) {
- Xrf <- matrix(Xrf,pd,tt+1)
- Xrp <- matrix(Xrp,pd,tt)
- IndIO <- as.logical(IndIO)
- IndAO <- as.logical(IndAO)
- }}
+
list(Xf = Xf, Xp = Xp, Xrf = Xrf, Xrp = Xrp,
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/robkalman -r 26
More information about the Robkalman-commits
mailing list