[Robkalman-commits] r43 - in branches: . robKalman_2012 robKalman_2012/pkg robKalman_2012/pkg/robKalman robKalman_2012/pkg/robKalman/R robKalman_2012/pkg/robKalman/R/toDoB robKalman_2012/pkg/robKalman/R/toDoP robKalman_2012/pkg/robKalman/demo robKalman_2012/pkg/robKalman/inst robKalman_2012/pkg/robKalman/man robKalman_2012/www
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Jun 15 17:19:26 CEST 2012
Author: ruckdeschel
Date: 2012-06-15 17:19:26 +0200 (Fri, 15 Jun 2012)
New Revision: 43
Added:
branches/robKalman_2012/
branches/robKalman_2012/README
branches/robKalman_2012/pkg/
branches/robKalman_2012/pkg/robKalman/
branches/robKalman_2012/pkg/robKalman/DESCRIPTION
branches/robKalman_2012/pkg/robKalman/NAMESPACE
branches/robKalman_2012/pkg/robKalman/R/
branches/robKalman_2012/pkg/robKalman/R/Util.R
branches/robKalman_2012/pkg/robKalman/R/UtileKF.R
branches/robKalman_2012/pkg/robKalman/R/allClass.R
branches/robKalman_2012/pkg/robKalman/R/calibrateRLS.R
branches/robKalman_2012/pkg/robKalman/R/classEKF.R
branches/robKalman_2012/pkg/robKalman/R/classUKF.R
branches/robKalman_2012/pkg/robKalman/R/generate.R
branches/robKalman_2012/pkg/robKalman/R/pkgAlt/
branches/robKalman_2012/pkg/robKalman/R/rLSeKF.R
branches/robKalman_2012/pkg/robKalman/R/rLSuKF.R
branches/robKalman_2012/pkg/robKalman/R/recESmoother.R
branches/robKalman_2012/pkg/robKalman/R/recFilter.R
branches/robKalman_2012/pkg/robKalman/R/recSmoother.R
branches/robKalman_2012/pkg/robKalman/R/simulateSScont-SO.R
branches/robKalman_2012/pkg/robKalman/R/simulateSScontEKF.R
branches/robKalman_2012/pkg/robKalman/R/toDoB/
branches/robKalman_2012/pkg/robKalman/R/toDoB/ACMfilt.R
branches/robKalman_2012/pkg/robKalman/R/toDoB/ACMfilter.R
branches/robKalman_2012/pkg/robKalman/R/toDoB/Psi.R
branches/robKalman_2012/pkg/robKalman/R/toDoB/mACMfilter.R
branches/robKalman_2012/pkg/robKalman/R/toDoB/mACMinternal.R
branches/robKalman_2012/pkg/robKalman/R/toDoD/
branches/robKalman_2012/pkg/robKalman/R/toDoP/
branches/robKalman_2012/pkg/robKalman/R/toDoP/0AllClass.R
branches/robKalman_2012/pkg/robKalman/R/toDoP/AllGenerics.R
branches/robKalman_2012/pkg/robKalman/R/toDoP/AllInitialize.R
branches/robKalman_2012/pkg/robKalman/R/toDoP/AllPlot.R
branches/robKalman_2012/pkg/robKalman/R/toDoP/AllShow.R
branches/robKalman_2012/pkg/robKalman/R/toDoP/allClass.R
branches/robKalman_2012/pkg/robKalman/R/toDoP/allClasses.R
branches/robKalman_2012/pkg/robKalman/R/toDoP/allGenerics (2).R
branches/robKalman_2012/pkg/robKalman/demo/
branches/robKalman_2012/pkg/robKalman/demo/00Index
branches/robKalman_2012/pkg/robKalman/demo/ACMdemo.R
branches/robKalman_2012/pkg/robKalman/demo/rLSdemo.R
branches/robKalman_2012/pkg/robKalman/inst/
branches/robKalman_2012/pkg/robKalman/inst/NEWS
branches/robKalman_2012/pkg/robKalman/man/
branches/robKalman_2012/pkg/robKalman/man/0robKalman-package.Rd
branches/robKalman_2012/pkg/robKalman/man/ACMfilt.Rd
branches/robKalman_2012/pkg/robKalman/man/arGM.Rd
branches/robKalman_2012/pkg/robKalman/man/calibrateRLS.Rd
branches/robKalman_2012/pkg/robKalman/man/internalACM.Rd
branches/robKalman_2012/pkg/robKalman/man/internalKalman.Rd
branches/robKalman_2012/pkg/robKalman/man/internalarGM.Rd
branches/robKalman_2012/pkg/robKalman/man/internalpsi.Rd
branches/robKalman_2012/pkg/robKalman/man/internalrLS.Rd
branches/robKalman_2012/pkg/robKalman/man/recFilter.Rd
branches/robKalman_2012/pkg/robKalman/man/simulateSScont.Rd
branches/robKalman_2012/pkg/robKalman/man/util.Rd
branches/robKalman_2012/www/
branches/robKalman_2012/www/index.php
Log:
robKalman_2012 committed: the root for the next version...
Added: branches/robKalman_2012/README
===================================================================
--- branches/robKalman_2012/README (rev 0)
+++ branches/robKalman_2012/README 2012-06-15 15:19:26 UTC (rev 43)
@@ -0,0 +1,49 @@
+ R-Forge SVN README
+
+
+(See "http://download.r-forge.r-project.org/manuals/R-Forge_Manual.pdf"
+ for detailed information on registering a new project.
+
+1. Introduction
+-----------------------------------------------------------------------
+R is free software distributed under a GNU-style copyleft. R-Forge is
+a service for R users and package developers providing certain tools
+for collaborative source code management.
+
+2. The directory you're in
+-----------------------------------------------------------------------
+This is the repository of your project. It contains two important
+pre-defined directories namely 'www' and 'pkg'. They must not be
+deleted otherwise R-Forge's core functionality will not be available
+(daily check and build of your package or project websites).
+These two directories are standardized and therefore are going to be
+described in this README. The rest of your repository can be used as
+you like.
+
+3. 'pkg' directory
+-----------------------------------------------------------------------
+Typically this directory contains the R package with the usual
+DESCRIPTION and R/, man/, data/ directories etc (see 'Writing R
+Extension' for more details). In the future it will also be possible to
+have multiple packages managed by a control file, however currently
+this feature is still under development).
+
+Furthermore, this directory will be checked out daily, the package is
+checked and if it passes this procedure it is build and made available at
+http://R-Forge.R-project.org/src/contrib/ (as source tar.gz and win32
+.zip). It should be possible to install the package via
+install.packages("foo",url="R-Forge.R-project.org") within R
+then.
+
+4. 'www' directory
+-----------------------------------------------------------------------
+This directory contains your project homepage which is available at
+http://<projectname>.R-Forge.R-project.org.
+Note that it will be checked out daily, so please take
+into consideration that it will not be available right after you
+commit your changes or updates.
+
+5. Help
+-----------------------------------------------------------------------
+If you need help don't hesitate to contact us
+(R-Forge at R-project.org)
Added: branches/robKalman_2012/pkg/robKalman/DESCRIPTION
===================================================================
--- branches/robKalman_2012/pkg/robKalman/DESCRIPTION (rev 0)
+++ branches/robKalman_2012/pkg/robKalman/DESCRIPTION 2012-06-15 15:19:26 UTC (rev 43)
@@ -0,0 +1,14 @@
+Package: robKalman
+Version: 1.0
+Date: 2011-08-14
+Title: Robust Kalman Filtering
+Description: Routines for Robust Kalman Filtering --- the ACM- and rLS-filter
+Author: Peter Ruckdeschel, Bernhard Spangl
+Maintainer: Peter Ruckdeschel <Peter.Ruckdeschel at itwm.fraunhofer.de>
+Depends: R(>= 2.3.0), methods, graphics, startupmsg, MASS, robustbase, numDeriv, robust-ts
+Imports: stats, MASS
+LazyLoad: yes
+License: LGPL-3
+URL: http://robkalman.r-forge.r-project.org/
+LastChangedDate: {$LastChangedDate: 2009-06-15 15:33:24 +0200 (Mo, 15 Jun 2009) $}
+LastChangedRevision: {$LastChangedRevision: 33 $}
Added: branches/robKalman_2012/pkg/robKalman/NAMESPACE
===================================================================
--- branches/robKalman_2012/pkg/robKalman/NAMESPACE (rev 0)
+++ branches/robKalman_2012/pkg/robKalman/NAMESPACE 2012-06-15 15:19:26 UTC (rev 43)
@@ -0,0 +1,13 @@
+import("methods")
+import("stats")
+import("MASS")
+import("startupmsg")
+import("robustbase")
+import("numDeriv")
+import("robust-ts")
+
+export("ACMfilt", "ACMfilter", "arGM", "EuclideanNorm",
+ "simulateState", "simulateObs", "rcvmvnorm", "Huberize",
+ "rLScalibrateB", "limitS", "rLSFilter", "rLS.IO.Filter",
+ "rLS.AO.Filter", "KalmanFilter",
+ "recursiveFilter")
Added: branches/robKalman_2012/pkg/robKalman/R/Util.R
===================================================================
--- branches/robKalman_2012/pkg/robKalman/R/Util.R (rev 0)
+++ branches/robKalman_2012/pkg/robKalman/R/Util.R 2012-06-15 15:19:26 UTC (rev 43)
@@ -0,0 +1,165 @@
+#######################################################
+##
+## rLS Kalman filter routines
+## [cf. Ruckdeschel, 2009]
+## author: Peter Ruckdeschel,
+## version: limitS taken from r-forge 0.3
+## rest: (P.R., 2012-04-10)
+##
+#######################################################
+
+limitS <- function(S, F, Z, Q, V, B = NULL, D = NULL, tol = 10^-4, itmax = 1000)#
+## determines lim_{t->infty} S_{t|t-1}
+#------------------------------------
+# inputs:
+#------------------------------------
+# S = Cov(x_0) (a matrix!)
+# F, Z, Q, V Hyperparameters as matrices
+# B, D notation from EKF by default set to unit matrices in code
+# tol: accuracy : when change from S_{t|t-1} to S_{t+1|t}
+# smaller than tol in abs. value, we stop
+# itmax: maximal number of iterations
+#------------------------------------
+# outputs:
+#------------------------------------
+# lim_{t->infty} S_{t|t-1} (as matrix)
+#------------------------------------
+ {SO0 <- S + 1
+ S0 <- S
+ if(is.null(B)) B <- diag(length(Q)^.5)
+ if(is.null(D)) D <- diag(length(V)^.5)
+
+
+ i <- 0
+# print(S0)
+ while( (i<3 ||(sum( (SO0 - S0)^2 ) > tol^2)) && (i < itmax) )
+ {i <- i + 1
+ S1 <- .getpredCov(S0, F, B, Q)
+ SO0 <- S0
+ K <- .getKG(S1, Z, .getDelta(S1, Z, D, V))
+ S0 <- .getcorrCov(S1, K, Z)
+
+# print(list(S0,i))
+ }
+ S1
+ }
+
+.timeInvModel <- function(S, F, Z, Q, V, B = NULL, D = NULL, tol = 10^-4,
+ itmax = 1000, repl = 100000, eff=NULL, r=0.1,
+ rlow = 0, rup = NULL, upto = 20, blow = 1e-6,
+ IO = TRUE, AO = TRUE, seed = NULL,
+ verbose = FALSE)#
+{#
+## determines lim_{t->infty} b_{t} in time-invariant model
+#------------------------------------
+# inputs:
+#------------------------------------
+# S = Cov(x_0) (a matrix!)
+# F, Z, Q, V Hyperparameters as matrices
+# B, D notation from EKF by default set to unit matrices in code
+# tol, itmax: as in limitS
+# blow, rlow, rup, upto, seed: as in rLScalibrateB
+# IO: do we compute b_IO (TRUE) or not (FALSE)
+# AO: do we compute b_AO (TRUE) or not (FALSE)
+# verbose: shall we output intermediate results?
+#------------------------------------
+# outputs:
+#------------------------------------
+# list with items
+# Slim = lim_t Sigma_{t|t-1}
+# b.AO, b.IO (respective clipping heights),
+# r.AO, r.IO (respective radii)
+#------------------------------------
+
+ #
+ Slim <- limitS(S=S, F=F, Z=Z, Q=Q, V=V, B=B, D=D, tol=tol, itmax=itmax)
+ if(verbose) print(Slim)
+ if(AO){
+ br.AO <- rLScalibrateB(Z = Z, S = Slim, V = V, repl = repl,
+ b = NULL, eff = eff, r = r,
+ rlow = rlow, rup = rup, upto = upto,
+ blow = blow, IO = FALSE,
+ seed = seed, verbose = verbose)
+ b.AO <- br.AO$b
+ r.AO <- br.AO$r
+ } else {
+ b.AO <- NULL
+ r.AO <- NULL
+ }
+ if(IO){
+ br.IO <- rLScalibrateB(Z = Z, S = Slim, V = V, repl = repl,
+ b = NULL, eff = eff, r = r,
+ rlow = rlow, rup = rup, upto = upto,
+ blow = blow, IO = TRUE,
+ seed = seed, verbose = verbose)
+ b.IO <- br.IO$b
+ r.IO <- br.IO$r
+ } else {
+ b.IO <- NULL
+ r.IO <- NULL
+ }
+ list(Slim = Slim, b.AO = b.AO, b.IO = b.IO, r.AO = r.AO, r.IO = r.IO)
+}
+
+.getb.rls <- function(b = NULL, controlCorr, Z, S, V, D=NULL,IO=FALSE){
+## determines b_{t} in time-variant model / in particular in EKF, UKF
+#------------------------------------
+# inputs:
+#------------------------------------
+# b clipping height: if given, computes r and eff else is NULL
+# controlCorr: control list for correction step rLS
+# Z, V Hyperparameters as matrices
+# S = Sigma_{t|t-1} (a matrix!)
+# D notation from EKF by default set to unit matrices in code
+# IO: do we compute b_IO (TRUE) or b_AO (FALSE)
+#------------------------------------
+# outputs:
+#------------------------------------
+# modified control list (with new entries b, eff, r in
+# list items b.v, eff.v, r.v
+#------------------------------------
+ ## vorlaeufig explizite Setzung /
+ ## mittelfristig auslesen in einer Paketoption vgl options()
+ bcal <- r <- eff <- NULL
+ if((is.null(controlCorr))||(is.null(names(controlCorr)))){
+ repl <- 10000; blow <- 1e-6; rlow <- 0; rup <- NULL
+ upto <- 20; seed <- NULL; verbose <- FALSE; r <- 0.1
+ controlCorr <- vector("list",0); nam <- NULL
+ }else{
+ nam <- names(controlCorr)
+ if(! "repl" %in% nam) repl <- 10000
+ if(! "blow" %in% nam) blow <- 1e-6
+ if(! "rlow" %in% nam) rlow <- 0
+ if(! "rup" %in% nam) rup <- NULL
+ if(! "upto" %in% nam) upto <- 20
+ if(! "seed" %in% nam) seed <- NULL
+ if(! "verbose" %in% nam) verbose <- FALSE
+ if(! (("eff" %in% nam) ||("r" %in% nam))) r <- 0.1
+ }
+ if(is.null(b)){
+ with(controlCorr,
+ { bcal <<- rLScalibrateB(Z = Z, S = S, V = V, D = D, repl = repl,
+ b = NULL, eff = eff, r = r, blow = blow,
+ rlow = rlow, rup = rup, upto = upto, IO = IO,
+ seed = seed, verbose = verbose)
+ }
+ )
+ }else{
+ with(controlCorr,
+ bcal <<- rLScalibrateB(Z = Z, S = S, V = V, D = D, repl = repl,
+ b = b, eff = NULL, r = NULL, blow = blow,
+ rlow = rlow, rup = rup, upto = upto, IO = IO,
+ seed = seed, verbose = verbose)
+ )
+ }
+ if(!is.null(controlCorr) && "b.v" %in% nam)
+ controlCorr$b.v <- c(controlCorr[["b.v"]], bcal$b)
+ else{controlCorr$b.v <- bcal$b}
+ if(!is.null(controlCorr) && "eff.v" %in% nam)
+ controlCorr$eff.v <- c(controlCorr[["eff.v"]], bcal$eff)
+ else{controlCorr$eff.v <- bcal$eff}
+ if(!is.null(controlCorr) && "r.v" %in% nam)
+ controlCorr$r.v <- c(controlCorr[["r.v"]], bcal$r)
+ else{controlCorr$r.v <- bcal$r}
+ return(controlCorr)
+}
\ No newline at end of file
Added: branches/robKalman_2012/pkg/robKalman/R/UtileKF.R
===================================================================
--- branches/robKalman_2012/pkg/robKalman/R/UtileKF.R (rev 0)
+++ branches/robKalman_2012/pkg/robKalman/R/UtileKF.R 2012-06-15 15:19:26 UTC (rev 43)
@@ -0,0 +1,24 @@
+#######################################################
+##
+## extended Kalman filter utility routines
+## [original code by Peter Ruckdeschel]
+## author: Bernhard Spangl
+## version: 0.2 (changed: 2011-12-16, created: 2010-09-07)
+##
+#######################################################
+
+## the Euclidean norm
+
+Euclideannorm <- function(x) {if(is.null(dim(x)))
+ sqrt(abs(sum(x^2))) else sqrt(colSums(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 + (nx+Ind1==0)) )
+}
+
+
Added: branches/robKalman_2012/pkg/robKalman/R/allClass.R
===================================================================
--- branches/robKalman_2012/pkg/robKalman/R/allClass.R (rev 0)
+++ branches/robKalman_2012/pkg/robKalman/R/allClass.R 2012-06-15 15:19:26 UTC (rev 43)
@@ -0,0 +1,94 @@
+########################################################
+##
+## S4 classes for robust filtering
+## author: Bernhard Spangl & Peter Ruckdeschel
+## version: 0.1 (last changed: 2012-04-27, created: 2011-08-19)
+##
+#######################################################
+
+
+.onLoad <- function(lib, pkg){
+}
+
+.onAttach <- function(library, pkg)
+{
+buildStartupMessage(pkg = "robKalman", msga, msgb, library = library, packageHelp = TRUE,
+# MANUAL="http://www.uni-bayreuth.de/departments/math/org/mathe7/DISTR/distr.pdf",
+VIGNETTE = gettext("Package \"robKalman\" provides a vignette to this package; try vignette(\"robKalman\")."))
+ invisible()
+}
+
+## ClassUnion: OptionalList
+setClassUnion("OptionalList",
+ c("list","NULL")
+ )
+
+
+## Class: FunctionWithControl
+setClass("FunctionWithControl",
+ representation = representation(fct = "function",
+ dots = "OptionalList",
+ control = "OptionalList",
+ name = "character"),
+ prototype = prototype(fct = function(x)x, dots = NULL, control = NULL,
+ name = gettext("a function with control"))
+ )
+setClass("SSVar", contains="FunctionWithControl")
+setClass("SSTransform", contains="FunctionWithControl")
+setClass("SSExo", contains="FunctionWithControl")
+setClass("SSstateEq",
+ representation = representation(F = "SSTransform",
+ Q = "SSVar",
+ Exo = "SSVar",
+ mu = "function",
+ distr = "Distribution"),
+ prototype = prototype(F = new("SSTransform",
+ fct=function(...)1, control = NULL,
+ name="state transition"),
+ Q = new("SSVar",
+ fct=function(...)1, control = NULL,
+ name="state variance"),
+ Exo = new("SSExo",
+ fct=function(...)1, control = NULL,
+ name="state Exogenous variable"),
+ mu = function(...)0,
+ distr = Norm())
+ )
+setClass("SSobsEq",
+ representation = representation(Z = "SSTransform",
+ V = "SSVar",
+ Exo = "SSVar",
+ mu = "function",
+ distr = "Distribution"),
+ prototype = prototype(Z = new("SSTransform",
+ fct=function(...)1, control = NULL,
+ name="state transition"),
+ V = new("SSVar",
+ fct=function(...)1, control = NULL,
+ name="state variance"),
+ Exo = new("SSExo",
+ fct=function(...)1, control = NULL,
+ name="state Exogenous variable"),
+ mu = function(...)0,
+ distr = Norm())
+ )
+setClass("SSstartEq",
+ representation = representation(a0 = "numeric",
+ Sigma0 = "matrix",
+ Exo = "SSVar",
+ mu = "function",
+ distr = "Distribution"),
+ prototype = prototype(a0 = 1,
+ Sigma0 = matrix(1,1,1),
+ Exo = new("SSExo",
+ fct=function(...)0, control = NULL,
+ name="state Exogenous variable"),
+ mu = function(...)0,
+ distr = Norm())
+
+ )
+setClass("SSmod",
+ representation = representation(StartEq = "SSstartEq",
+ StatesEq = "SSstateEq",
+ ObsEq = "SSobsEq")
+)
Added: branches/robKalman_2012/pkg/robKalman/R/calibrateRLS.R
===================================================================
--- branches/robKalman_2012/pkg/robKalman/R/calibrateRLS.R (rev 0)
+++ branches/robKalman_2012/pkg/robKalman/R/calibrateRLS.R 2012-06-15 15:19:26 UTC (rev 43)
@@ -0,0 +1,169 @@
+#######################################################
+##
+## rLS Kalman filter routines
+## [cf. Ruckdeschel, 2009]
+## author: Peter Ruckdeschel,
+## version: taken from r-forge 0.3
+## adopted for EKF (P.R., 2012-04-10)
+##
+#######################################################
+
+rLScalibrateB <- function(Z, S, V, D= NULL, repl = 100000, b = NULL,
+ eff = NULL, r = NULL, rlow = 0, rup = NULL, upto = 20,
+ blow = 0, IO = FALSE, seed = NULL,
+ verbose = FALSE)#
+# calibrates clipping height b to given Z, V, and S_{t|t-1}
+# ---
+# either to given efficiency in the ideal model
+# or to given (SO)-radius about the ideal model
+# expectations are calculated by a LLN approximation
+#------------------------------------
+# inputs:
+#------------------------------------
+# r: SO[AO/IO]-radius to calibrate for; may be NULL
+# eff: efficiency (< 1!) to calibrate for; may be NULL
+# b: clipping height to use; may be NULL
+# S = Sigma_{t|t-1} (a matrix!)
+# Z, V Hyperparameters as matrices
+# D notation from EKF by default set to unit matrices in code
+# repl: number of MC replicates; 100000 as default works well for me,
+# but might need a change for slower machines
+# blow: lower bound for b
+# rlow, rup: search interval for r (scalar valued each, resp. rup may be NULL)
+# upto: upper bound where to search for zeroes (scalar valued)
+# IO: do we compute b_IO (TRUE) or b_AO (FALSE)
+# seed: seed for MC integration
+# verbose: shall we output intermediate results?
+#------------------------------------
+# outputs:
+#------------------------------------
+# list with items b (clipping height), eff (efficiency), r (radius)
+{
+
+ if ( is.null(b) && is.null(eff) && is.null(r) && is.null(rup))
+ stop("You must either specify argument 'b' or 'r' or 'eff' or 'rup'")
+
+ qd <- ifelse(length(Z)==1, 1, (dim(Z))[1])
+ pd <- ifelse(length(S)==1, 1, (dim(Z))[2])
+
+ if(!is.null(seed)) set.seed(seed)
+
+ if(is.null(D)) D <- diag(qd)
+
+ dx <- t(mvrnorm(repl, numeric(pd), S))
+ eps <- t(mvrnorm(repl, numeric(qd), V))
+ dy <- Z %*% dx + eps
+ gD <- .getDelta(S, Z, D ,V)
+ K <- .getKG(S, Z, gD)
+
+# trS <- sum(diag(.getcorrCov(S, K, Z)))
+
+ dx0 <- K %*% dy
+ trS <- sum( (dx-dx0)^2 )/repl
+
+ if(IO){
+ dx0 <- dy - Z %*% dx0
+ Zs0 <- S %*% t(Z)
+ Zs1 <- Z%*%Zs0
+ Zsm <- ginv(Zs1)
+ Zs <- Zs0 %*% Zsm
+# piq <- diag(qd)-Zs1%*%Zsm
+# Kq <- Zs0 %*% piq %*% gD
+# dxq <- t(Kq %*% dy)
+ }#else dxq <- 0 * t(dx0)
+
+ no <- sqrt(colSums(dx0^2))
+
+ eff.b <- function(b){
+ w <- ifelse(no < b, 1, b / no)
+ dxw <- as.vector(w) * t(dx0)
+ if(IO) dxw <- (t(dy) - dxw) %*% t(Zs)
+ trSb <- sum( (t(dx) - dxw)^2 )/repl
+ #if(verbose) print(c(trSb,trS / trSb))
+ trS / trSb
+ }
+ r.b <- function(b){
+ ex <- mean( pmax(no/b - 1, 0) )
+ if(verbose) print(ex/(ex+1))
+ ex/(ex + 1)
+ }
+
+ todo.eff <- TRUE
+ todo.r <- TRUE
+
+ if(is.null(b)){
+ todo.r.search <- FALSE
+ if( (is.null(r)&& is.null(eff))) {
+ todo.r.search <- TRUE
+ r <- rup
+ }
+ if (is.null(r)&&!is.null(eff)){ ## calibrated to given efficiency
+ 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(Zs)
+ trSb <- sum( (t(dX) - dxw)^2 )/repl0
+ #if(verbose) print(c(b, trSb,trS / trSb, eff0))
+ trS0 / trSb - eff0
+ }
+ r1 <- NULL
+ eff1 <- eff
+ todo.eff <- FALSE
+ }else{ ## calibrated to given radius
+ f <- function(b, dX = dx, dX0 = dx0, no0 = no, r0 = r,
+ eff0 = eff, trS0 = trS, repl0 = repl, dY = dy){
+ er <- (1 - r0)/r0 * sum(pmax(no0 - b, 0))/repl0 - b
+ if(verbose) print(c(er=er))
+ er
+ }
+ r1 <- r
+ eff1 <- NULL
+ todo.r <- FALSE
+ }
+
+ if(abs(f(0))<1e-16){
+ if(is.null(eff)) eff <- 1
+ if(is.null(r)) r <- 1
+ return(list( b = 0, eff = eff, r = r ))
+ }
+ b <- uniroot(f, interval = c(blow, upto*max(sqrt(trS),1)), tol = 10^-7,
+ dX = dx, dX0 = dx0, no0 = no,
+ eff0 = eff1, trS0 = trS, repl0 = repl, r0 = r1,
+ dY = dy)$root
+
+ if(! todo.r.search){
+ if (is.null(r)) ### corresponding radius is calculated
+ r <- r.b(b)
+ else ### corresponding effciency is calculated
+ eff <- eff.b(b)
+ }else{
+# if(IO) stop("not yet implemented")
+ todo.r <- TRUE
+ b.u <- b
+
+ A.r <- function(b1) trS+ mean((pmax(no-b1,0))^2)
+ B.r <- function(b1) mean(no^2)- mean((pmax(no-b1,0)^2)) + b1^2
+
+ B.u <- B.r(b.u)
+
+
+ if(rlow <1e-6){
+ b.l <- 1e8
+ A.l <- trS
+ }else{
+ b.l <- uniroot(f, interval = c(blow, upto*max(1,sqrt(trS))), tol = 10^-7,
+ dX = dx, dX0 = dx0, no0 = no,
+ eff0 = eff1, trS0 = trS, repl0 = repl, r0 = rlow,
+ dY = dy)$root
+ A.l <- A.r(b.l)
+ }
+ AB <- function(b2) A.r(b2)/A.l - B.r(b2)/B.u
+ b <- uniroot(AB,interval = c(b.l,b.u), tol = 10^-7)$root
+ }
+ }
+ if(todo.r) r <- r.b(b)
+ if(todo.eff) eff <- eff.b(b)
+
+ list( b = b, eff = eff, r = r )
+}
Added: branches/robKalman_2012/pkg/robKalman/R/classEKF.R
===================================================================
--- branches/robKalman_2012/pkg/robKalman/R/classEKF.R (rev 0)
+++ branches/robKalman_2012/pkg/robKalman/R/classEKF.R 2012-06-15 15:19:26 UTC (rev 43)
@@ -0,0 +1,113 @@
+#######################################################
+##
+## classical extended Kalman filter routines
+## author: Bernhard Spangl
+## version: 0.3 (changed: 2011-12-13, created: 2011-06-10)
+##
+#######################################################
+
+.getDelta <- function (S1, C, D, V)
+{
+## calculates the Cov(Delta y_t)
+## for S1 = S_{t|t-1}, C (=Z), D (=Id), V as above
+ H <- S1 %*% t(C)
+ ginv( C %*% H + D %*% V %*% t(D) )
+}
+
+.getKG <- function (S1, Z, Delta)
+{
+## calculates the Kalman Gain for S1 = S_{t|t-1}, Z, V as above
+ S1 %*% t(Z) %*% 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
+}
+
+.getpredCov <- function (S0, A, B, Q)
+{
+## calculates S_{t|t-1} for S0 = S_{t-1|t-1}, A (=F), B (=Id), Q as above
+ A %*% S0 %*% t(A) + B %*% Q %*% t(B)
+}
+
+.cinitstep <- function (a, S, controlInit, ...)
+{
+## a ... E(x_0)
+## S ... S_{0|0} = Cov(x_0 - x_{0|0})
+## = E((x_0 - x_{0|0})(x_0 - x_{0|0})^\top), error covariance
+## controlInit ... control parameters of used filter
+
+ return(list(x0=a, S0=S, controlInit=controlInit))
+
+}
+
+.cpredstep <- function (x0, S0, F, Q, i,
+ v, u, controlF, # arguments of F
+ exQ, controlQ, # arguments of Q
+ controlPred, ...) # arguments of used filter
+{
+## x0 ... x_{t-1|t-1}, filter estimate
+## S0 ... S_{t-1|t-1}, conditional filtering error covariance matrix
+## F ... F(t, x_{t-1}, u_t, v_t, control), function of state equation
+## Q ... Q(t, x_{t-1}, exQ_{t-1}, control), function of innovations cov-matrix
+## i ... time index
+## v ... innovations v_t
+## u ... u_{t-1}, exogenous variables of F
+## controlF ... control parameters of F
+## exQ ... exQ_{t-1}, exogenous variables of Q
+## controlQ ... control parameters of Q
+## controlPred ... control parameters of used filter
+ Freturn <- F(t=i, x0=x0, v=v, u=u, control=controlF)
+ x1 <- Freturn$x1
+ A <- Freturn$A
+ B <- Freturn$B
+
+ Qreturn <- Q(t=i, x0=x0, exQ=exQ, control=controlQ)
+ Q <- Qreturn$Q
+
+ return(list(x1=x1, S1=.getpredCov(S0=S0, A=A, B=B, Q=Q),
+ controlPred=controlPred))
+
+}
+
+.ccorrstep <- function (y, x1, S1, Z, V, i,
+ eps, w, controlZ, # arguments of Z
+ exV, controlV, # arguments of V
+ controlCorr, ...) # arguments of used filter
+{
+## y ... observations
+## x1 ... x_{t|t-1}, one-step-ahead prediction
+## S1 ... S_{t|t-1}, conditional prediction error covariance matrix
+## Z ... Z(t, x_t, eps_t, w_t, control), function of observation equation
+## V ... V(t, x_t, exV_t, control), function of cov-matrix of observation error
+## i ... time index
+## eps ... observation error \eps_t
+## w ... exogenous variable w_t
+## controlZ ... control parameters of Z
+## exV ... exV_t, exogenous variables of V
+## controlV ... control parameters of V
+## controlCorr ... control parameters of used filter
+
+ Zreturn <- Z(t=i, x1=x1, eps=eps, w=w, control=controlZ)
+ yhat <- Zreturn$y
+ C <- Zreturn$C
+ D <- Zreturn$D
+
+ Vreturn <- V(t=i, x1=x1, exV=exV, control=controlV)
+ V <- Vreturn$V
+
+ Delta <- .getDelta(S1=S1, C=C, D=D, V=V)
+ K <- .getKG(S1=S1, Z=C, Delta=Delta)
+ DeltaY <- y - yhat
+
+ x0 <- x1 + K %*% DeltaY
+ S0 <- .getcorrCov(S1=S1, K=K, Z=C)
+
+ return(list(x0=x0, K=K, S0=S0, Delta=Delta, DeltaY=DeltaY,
+ controlCorr=controlCorr))
+
+}
+
+
Added: branches/robKalman_2012/pkg/robKalman/R/classUKF.R
===================================================================
--- branches/robKalman_2012/pkg/robKalman/R/classUKF.R (rev 0)
+++ branches/robKalman_2012/pkg/robKalman/R/classUKF.R 2012-06-15 15:19:26 UTC (rev 43)
@@ -0,0 +1,218 @@
+#######################################################
+##
+## classical unscented Kalman filter routines
+## author: Bernhard Spangl
+## version: 0.7 (changed: 2012-01-07, created: 2011-06-17)
+##
+#######################################################
+
+.SigmaPoints <- function (x, gamma, Sigma)
+{
+## x ... vector
+## gamma ... scaling parameter
+## Sigma ... covariance matrix
+ sqrtSigma <- t(chol(Sigma)) # FIXME! (doesn't work for
+ # positiv semi-definite matricies)
+## sqrtSigma <- diag(sqrt(diag(Sigma))) # Mist!
+## Q <- chol(Sigma, pivot=TRUE) # FIXED!
+## pivot <- attr(Q, "pivot")
+## sqrtSigma <- t(Q[, order(pivot)])
+## ch <- gchol(Sigma) # better!!!
+## L <- as.matrix(ch)
+## sqrtD <- diag(sqrt(diag(ch)))
+## sqrtSigma <- L %*% sqrtD
+ SP1 <- as.vector(x) %o% rep(1, ncol(Sigma)) + gamma * sqrtSigma
+ SP2 <- as.vector(x) %o% rep(1, ncol(Sigma)) - gamma * sqrtSigma
+ cbind(x, SP1, SP2)
+}
+
+## .SigmaPointsCalc <- function (x0, S0, Q, V, i,
+## v, eps,
+## exQ, controlQ, # arguments of Q
+## exV, controlV, # arguments of V
+## controlSigma, ...) # arguments of used filter
+## {
+## ## x0 ... x_{t-1|t-1}, filter estimate
+## ## S0 ... S_{t-1|t-1}, conditional filtering error covariance matrix
+## ## Q ... Q(t, x_{t-1}, exQ_{t-1}, control), function of innovations cov-matrix
+## ## V ... V(t, x_{t-1}, exV_t, control), function of cov-matrix of
+## ## observation error
+## ## i ... time index
+## ## v ... innovations v_t
+## ## eps ... observation error \eps_t
+## ## exQ ... exQ_{t-1}, exogenous variables of Q
+## ## controlQ ... control parameters of Q
+## ## exV ... exV_t, exogenous variables of V
+## ## controlV ... control parameters of V
+## ## controlSigma ... control parameters of used filter
+##
+## gamma <- controlSigma$gamma
+## pd <- controlSigma$pd
+## rd <- controlSigma$rd
+## td <- controlSigma$td
+##
+## Qreturn <- Q(t=i, x0=x0, exQ=exQ, control=controlQ)
+## Q <- Qreturn$Q
+## Vreturn <- V(t=i, x1=x0, exV=exV, control=controlV)
+## V <- Vreturn$V
+##
+## x <- c(x0, v, eps)
+## ## Sigma <- blockDiag(S0, Q, V)
+## Sigma <- bdsmatrix(c(pd, rd, td),
+## c(as.vector(S0), as.vector(Q), as.vector(V)))
+## SigmaPoints <- .SigmaPoints(x=x, gamma=gamma, Sigma=Sigma)
+## controlSigma$X0x <- SigmaPoints[1:pd, ]
+## controlSigma$Xv <- SigmaPoints[(pd+1):(pd+rd), ]
+## controlSigma$Xe <- SigmaPoints[(pd+rd+1):(pd+rd+td), ]
+##
+## return(controlSigma)
+##
+## }
+
+.CovEst <- function (X, muX, Y, muY, w)
+{
+## X ... matrix of sigma points
+## muX ... mean vector of matrix X
+## Y ... matrix of sigma points
+## muY ... mean vector of matrix Y
+## w ... vector of weights
+ X <- X - as.vector(muX) %o% rep(1, ncol(X))
+ X <- X * (rep(1, nrow(X)) %o% w)
+ Y <- Y - as.vector(muY) %o% rep(1, ncol(Y))
+ X %*% t(Y)
+}
+
+.Eval <- function (vec, G, t, err, ex, control)
+{
+ G <- G(t, vec, err, ex, control)
+ return(G[[1]])
+}
+
+.cUKFinitstep <- function (a, S, controlInit, ...)
+{
+## a ... E(x_0)
+## S ... S_{0|0} = Cov(x_0 - x_{0|0})
+## = E((x_0 - x_{0|0})(x_0 - x_{0|0})^\top), error covariance
+## controlInit ... list containing alpha, beta and kappa,
+## and dimensions of ssm
+
+ alpha <- controlInit$alpha
+ beta <- controlInit$beta
+ kappa <- controlInit$kappa
+## pd <- controlInit$pd
+
+ pd <- length(a)
+ L <- pd
+
+ lambda <- alpha^2*(L + kappa) - L
+ Wm <- Wc <- rep(1/(2*(L + lambda)), 2*L)
+ Wm <- c(lambda/(L + lambda), Wm)
+ Wc <- c(lambda/(L + lambda) + (1 + alpha^2 + beta), Wc)
+ gamma <- sqrt(L + lambda)
+
+ controlInit$gamma <-gamma
+ controlInit$Wm <- Wm
+ controlInit$Wc <- Wc
+
+ return(list(x0=a, S0=S, controlInit=controlInit))
+
+}
+
+.cUKFpredstep <- function (x0, S0, F, Q, i,
+ v, u, controlF, # arguments of F
+ exQ, controlQ, # arguments of Q
+ controlPred, ...) # arguments of used filter
+{
+## x0 ... x_{t-1|t-1}, filter estimate
+## S0 ... S_{t-1|t-1}, conditional filtering error covariance matrix
+## F ... F(t, x_{t-1}, v_t, u_t, control), function of state equation
+## Q ... Q(t, x_{t-1}, exQ_{t-1}, control), function of innovations cov-matrix
+## i ... time index
+## v ... innovations v_t
+## u ... u_{t-1}, exogenous variables of F
+## controlF ... control parameters of F
+## exQ ... exQ_{t-1}, exogenous variables of Q
+## controlQ ... control parameters of Q
+## controlPred ... control parameters of used filter
+
+ gamma <- controlPred$gamma
+ Wm <- controlPred$Wm
+ Wc <- controlPred$Wc
+## rd <- controlPred$rd # FIXME: Dimensionen bereits als Slot ins
+## # Funktionsobjekt (bzw. zugehörige
+## # S4-Klasse)!
+
+ Freturn <- F(t=i, x0=x0, v=v, u=u, control=controlF)
+ B <- Freturn$B
+
+ Qreturn <- Q(t=i, x0=x0, exQ=exQ, control=controlQ)
+ Q <- Qreturn$Q
+
+ rd <- ncol(Q)
+
+## S0 <- bdsmatrix(pd, as.vector(S0))
+ X0x <- .SigmaPoints(x=x0, gamma=gamma, Sigma=S0)
+
+ X1x <- apply(X0x, 2, .Eval,
+ G=F, t=i, err=rep(0, rd), ex=u, control=controlF)
+ x1 <- X1x %*% Wm
+ S1 <- .CovEst(X=X1x, muX=x1, Y=X1x, muY=x1, w=Wc) + B %*% Q %*% t(B)
+
+## controlPred$X1x <- X1x
+
+ return(list(x1=x1, S1=S1, controlPred=controlPred))
+
+}
+
+.cUKFcorrstep <- function (y, x1, S1, Z, V, i,
+ eps, w, controlZ, # arguments of Z
+ exV, controlV, # arguments of V
+ controlCorr, ...) # arguments of used filter
+{
+## y ... observations
+## x1 ... x_{t|t-1}, one-step-ahead prediction
+## S1 ... S_{t|t-1}, conditional prediction error covariance matrix
+## Z ... Z(t, x_t, eps_t, w_t, control), function of observation equation
+## V ... V(t, x_t, exV_t, control), function of cov-matrix of observation error
+## i ... time index
+## eps ... observation error \eps_t
+## w ... exogenous variable w_t
+## controlZ ... control parameters of Z
+## exV ... exV_t, exogenous variables of V
+## controlV ... control parameters of V
+## controlCorr ... control parameters of used filter
+
+ gamma <- controlCorr$gamma
+ Wm <- controlCorr$Wm
+ Wc <- controlCorr$Wc
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/robkalman -r 43
More information about the Robkalman-commits
mailing list