[Robkalman-commits] r42 - in branches: . robKalman-Coventry robKalman-Coventry/R robKalman-Coventry/chm robKalman-Coventry/demo robKalman-Coventry/inst robKalman-Coventry/man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Aug 16 02:05:28 CEST 2011
Author: ruckdeschel
Date: 2011-08-16 02:05:27 +0200 (Tue, 16 Aug 2011)
New Revision: 42
Added:
branches/robKalman-Coventry/
branches/robKalman-Coventry/DESCRIPTION
branches/robKalman-Coventry/NAMESPACE
branches/robKalman-Coventry/R/
branches/robKalman-Coventry/R/ACMfilt.R
branches/robKalman-Coventry/R/ACMfilter.R
branches/robKalman-Coventry/R/Psi.R
branches/robKalman-Coventry/R/Util.R
branches/robKalman-Coventry/R/calibrateRLS.R
branches/robKalman-Coventry/R/classEKF.R
branches/robKalman-Coventry/R/classKalman.R
branches/robKalman-Coventry/R/classUKF.R
branches/robKalman-Coventry/R/generate.R
branches/robKalman-Coventry/R/mACMfilter.R
branches/robKalman-Coventry/R/mACMinternal.R
branches/robKalman-Coventry/R/rLSfilter.R
branches/robKalman-Coventry/R/recEFilter.R
branches/robKalman-Coventry/R/recFilterInternal.R
branches/robKalman-Coventry/R/recFilterWrapper.R
branches/robKalman-Coventry/R/simulateSScont.R
branches/robKalman-Coventry/chm/
branches/robKalman-Coventry/chm/00Index.html
branches/robKalman-Coventry/chm/0robKalman-package.html
branches/robKalman-Coventry/chm/ACMfilt.html
branches/robKalman-Coventry/chm/Logo.jpg
branches/robKalman-Coventry/chm/Rchm.css
branches/robKalman-Coventry/chm/arGM.html
branches/robKalman-Coventry/chm/calibrateRLS.html
branches/robKalman-Coventry/chm/internalACM.html
branches/robKalman-Coventry/chm/internalKalman.html
branches/robKalman-Coventry/chm/internalarGM.html
branches/robKalman-Coventry/chm/internalpsi.html
branches/robKalman-Coventry/chm/internalrLS.html
branches/robKalman-Coventry/chm/recFilter.html
branches/robKalman-Coventry/chm/robKalman.chm
branches/robKalman-Coventry/chm/robKalman.hhp
branches/robKalman-Coventry/chm/robKalman.toc
branches/robKalman-Coventry/chm/simulateSScont.html
branches/robKalman-Coventry/chm/util.html
branches/robKalman-Coventry/demo/
branches/robKalman-Coventry/demo/00Index
branches/robKalman-Coventry/demo/ACMdemo.R
branches/robKalman-Coventry/demo/rLSdemo.R
branches/robKalman-Coventry/inst/
branches/robKalman-Coventry/inst/NEWS
branches/robKalman-Coventry/man/
branches/robKalman-Coventry/man/0robKalman-package.Rd
branches/robKalman-Coventry/man/ACMfilt.Rd
branches/robKalman-Coventry/man/calibrateRLS.Rd
branches/robKalman-Coventry/man/internalACM.Rd
branches/robKalman-Coventry/man/internalKalman.Rd
branches/robKalman-Coventry/man/internalpsi.Rd
branches/robKalman-Coventry/man/internalrLS.Rd
branches/robKalman-Coventry/man/recFilter.Rd
branches/robKalman-Coventry/man/simulateSScont.Rd
branches/robKalman-Coventry/man/util.Rd
Log:
Coventry-Version commited
Property changes on: branches/robKalman-Coventry
___________________________________________________________________
Added: bugtraq:number
+ true
Added: branches/robKalman-Coventry/DESCRIPTION
===================================================================
--- branches/robKalman-Coventry/DESCRIPTION (rev 0)
+++ branches/robKalman-Coventry/DESCRIPTION 2011-08-16 00:05:27 UTC (rev 42)
@@ -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-Coventry/NAMESPACE
===================================================================
--- branches/robKalman-Coventry/NAMESPACE (rev 0)
+++ branches/robKalman-Coventry/NAMESPACE 2011-08-16 00:05:27 UTC (rev 42)
@@ -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")
Property changes on: branches/robKalman-Coventry/R
___________________________________________________________________
Added: bugtraq:number
+ true
Added: branches/robKalman-Coventry/R/ACMfilt.R
===================================================================
--- branches/robKalman-Coventry/R/ACMfilt.R (rev 0)
+++ branches/robKalman-Coventry/R/ACMfilt.R 2011-08-16 00:05:27 UTC (rev 42)
@@ -0,0 +1,81 @@
+ACMfilt <- function (x, gm, s0=0,
+ psi="Hampel", a=2.5, b=a, c=5.0,
+ flag="weights", lagsmo=TRUE)
+{
+###########################################
+##
+## R-function: ACMfilt - approximate conditional-mean filtering (wrapper)
+## author: Bernhard Spangl
+## version: 1.1 (2007-08-13 and 2006-08-31)
+## References:
+## [Mart79c] R.D. Martin, Approximate Conditional-mean Type Smoothers
+## and Interpolators (1979)
+## [Mart81b] R.D. Martin, Robust Methods for Time Series (1981)
+## [MarT82b] R.D. Martin & D.J. Thomson, Robust-resistent Spectrum
+## Estimation (1982)
+##
+###########################################
+
+## Paramters:
+## x ... univariate time series (vector)
+## gm ... list as produced by function 'arGM' which includes components
+## 'ar' containing the AR(p) coefficient estimates, 'sinnov' containing
+## innovation scale estiamtes from AR(p) fits of orders 1 through p;
+## 'Cx' containing an estimate of the p by p autocovariance matrix,
+## and 'mu', the estimated mean of 'x'.
+## s0 ... scale of nominal Gaussian component of the additive noise
+## psi ... influence function to be used (default: "Hampel",
+## only Hampel's psi function 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
+## the weights; if "deriv", use psi'(t)
+## lagsmo ... logical, if TRUE (default) lag p-1 smoothing is performed;
+## if FALSE filtering from the top of ^X_t is performed
+
+## Variable definitions:
+
+ N <- length(x)
+ phi <- gm$ar
+ p <- length(phi)
+ si <- gm$sinnov[p]
+ Cx <- gm$Cx
+ Phi <- cbind(rbind(phi[-p], diag(rep(1, (p-1)))), c(phi[p], rep(0, (p-1))))
+ Q <- matrix(0, p, p)
+## Q <- diag(rep(0, p))
+ Q[1, 1] <- si^2
+
+ m0 <- rep(0, p)
+ H <- matrix(c(1, rep(0, (p-1))), 1, p)
+ V <- matrix(s0^2)
+ psi <- .psi(psi)
+
+ ## Centering:
+ x <- x - gm$mu
+ 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)]
+ st <- as.numeric(unlist(ACMres$rob1L))
+
+
+ if (!lagsmo) {
+ x.ck <- X.ck[1, ]
+ x <- X[1, ]
+ } else {
+ x.ck <- c(X.ck[p, p:N], X.ck[(p-1):1, N])
+ x <- c(X[p, p:N], X[(p-1):1, N])
+ }
+
+## ARmodel <- .ARmodel(x, p)
+## y <- ARmodel$y
+## Z <- ARmodel$Z
+## r <- resid(lm.fit(Z, y))
+
+ return(list(filt.ck=x.ck +gm$mu, filt=x + gm$mu, st=st)) #,
+## r=c(rep(NA, p), r)))
+
+}
Added: branches/robKalman-Coventry/R/ACMfilter.R
===================================================================
--- branches/robKalman-Coventry/R/ACMfilter.R (rev 0)
+++ branches/robKalman-Coventry/R/ACMfilter.R 2011-08-16 00:05:27 UTC (rev 42)
@@ -0,0 +1,91 @@
+.getcorrCovACM <- function (S1, K, Z, W=diag(nrow(Z)))
+{
+###########################################
+##
+## R-function: .corrCov - computes filtering error covarince matrix
+## (internal function)
+## author: Bernhard Spangl
+## version: 1.0 (2006-05-22)
+##
+###########################################
+
+## Paramters:
+## S1 ... prediction error covariance matrix
+## K ... Kalman gain
+## W ... weight matrix
+## Z ... observation matrix
+
+ S1 - K %*% W %*% Z %*% S1
+}
+
+##steps for classical Kalman filter (cK)
+.ACMinitstep <- function(a, S, ...)
+ {dots <- list(...)
+ if(hasArg("s0"))
+ s0<-dots$"s0"
+ else
+ s0<-NULL
+ list( x0 = a, S0 = S, s0 = s0)}
+
+.ACMpredstep <- function (x0, S0, F, Q, i, rob0, s0, ...) ### S=P F= Phi
+{
+###########################################
+##
+## R-function: .ACMpredstep - prediction step (internal function)
+## author: Bernhard Spangl
+## version: 1.0 (2006-05-22)
+##
+###########################################
+
+## Paramters:
+## x0 ... state vector (filter estimate)
+## F=Phi ... design matrix of state equation
+## S0 ... filtering error covariance matrix
+## 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=1))
+}
+
+.ACMcorrstep <- function (y, x1, S1, Z, V, i, rob1, dum=NULL, psi, apsi, bpsi, cpsi, flag, ...)
+{
+###########################################
+##
+## R-function: .ACMcorrstep - correction step (internal function)
+## author: Bernhard Spangl
+## version: 1.0 (2006-05-22)
+##
+###########################################
+
+## Paramters:
+## y ... univariate time series
+## x1 ... state vector (one-step-ahead predictor)
+## rob1 ... general robust parameter --- here st ... time-dependent scale parameter
+## S1 ... prediction error covariance matrix
+## Z ... observation matrix
+## dum ... dummy variable for compatibility with ... argument of calling function
+## V ... covariance matrix of observation noise
+## psi ... influence function to be used
+## 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)
+ st <- rob1
+
+ K <- .getKG(S1, Z, V)
+
+ rst <- (y - x1[1])/st
+
+ ps <- psi(rst, apsi, bpsi, cpsi)[1,1]
+ dx <- K * st * ps
+ x0 <- x1 + dx
+
+ ind <- (abs(rst-ps)>10^-8)
+
+ 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=Delta, Ind=ind, rob0=rob1, DeltaY = rst))
+}
Added: branches/robKalman-Coventry/R/Psi.R
===================================================================
--- branches/robKalman-Coventry/R/Psi.R (rev 0)
+++ branches/robKalman-Coventry/R/Psi.R 2011-08-16 00:05:27 UTC (rev 42)
@@ -0,0 +1,169 @@
+##############################################################
+# psi functions
+##############################################################
+
+
+psiLS <- function(t, rho=FALSE)
+{
+##################################
+##
+## R-function: psiLS - 'psi-function' for LS estimation
+## author: Bernhard Spangl
+## version: 1.0 (2006-05-26)
+##
+##################################
+
+## t ... input vector
+## rho ... logical, wether the psi- (default) or rho-function is used
+
+ if (!rho) {
+ t
+ } else {
+ t^2/2
+ }
+}
+
+psiHuber <- function (t, k=1.345, rho=FALSE)
+{
+##################################
+##
+## R-function: psiHuber - Huber's psi-function
+## author: Bernhard Spangl
+## version: 1.0 (2006-05-26)
+##
+##################################
+
+## t ... input vector
+## k ... tuning constant (default: k=1.345)
+## rho ... logical, wether the psi- (default) or rho-function is used
+
+ if (!rho) {
+ pmin(k, pmax(-k, t))
+ } else {
+ r <- abs(t)
+ i <- r < k
+ r[i] <- r[i]^2/2
+ r[!i] <- k * (r[!i] - k/2)
+ r
+ }
+}
+
+psiTukey <- function (t, c=4.685, rho=FALSE)
+{
+##################################
+##
+## R-function: psiTukey - Tukey's psi-function
+## author: Bernhard Spangl
+## version: 1.0 (2006-05-26)
+##
+##################################
+
+## t ... input vector
+## c ... tuning constant (default: c=4.685)
+## rho ... logical, wether the psi- (default) or rho-function is used
+
+ if (!rho) {
+ r <- abs(t)
+ i <- r < c
+ t[i] <- t[i]*(1-(t[i]/c)^2)^2
+ t[!i] <- 0
+ t
+ } else {
+ pmin((c^2/6)*(1-(1-(t/c)^2)^3), c^2/6)
+ }
+}
+
+
+psiHampel <- function (t, a=2, b=4, c=8, flag="psi")
+{
+###########################################
+##
+## R-function: psiHampel - Hampel's psi-function
+## author: Bernhard Spangl
+## version: 1.0 (2006-05-21)
+##
+###########################################
+
+## Paramters:
+## t ... input vector
+## a, b, c ... tuning constants
+## flag ... character, either "psi", "weights" or "deriv" to use
+## psi-function (default), weight function psi(t)/t or
+## its derivative respectively
+
+ at <- abs(t)
+ if (flag == "psi") {
+ dummy <- pmin(at, a, a/(c-b)*(c-at))
+ sign(t)*pmax(dummy, 0)
+ } else {
+ T <- pmin(at, c)
+ if (flag == "weights") {
+ ifelse(T <= a, 1, ifelse(T <= b, a/T, a*(c - T)/(c - b)/T))
+ } else if (flag == "deriv") {
+ ifelse(at <= c, ifelse(T <= a, 1, ifelse(T <= b, 0, -a/(c - b))), 0)
+ } else {
+ warning("error in function 'psiHampel': wrong flag type \n")
+ }
+ }
+}
+
+
+.psi <- function (type)
+{
+##################################
+##
+## R-function: .psi - switches to appropriate psi-function
+## (internal function)
+## author: Bernhard Spangl
+## version: 1.0 (2006-05-21)
+##
+##################################
+
+## type ... which psi-function, "Huber", "Tukey", , "Hampel", "Ident"
+
+ switch(type,
+ Ident = get("psiLS", mode="function"),
+ Huber = get("psiHuber", mode="function"),
+ Tukey = get("psiTukey", mode="function"),
+ Hampel = get("psiHampel", mode="function"))
+}
+
+mvpsiHampel <- function (x, a=2, b=4, c=8)
+{
+###########################################
+##
+## R-function: mvpsiHampel - multivariate analogue of
+## Hampel's psi-function
+## author: Bernhard Spangl
+## version: 0.1 (2008-02-23)
+##
+###########################################
+
+## Paramters:
+## x ... vector
+## a, b, c ... tuning constants
+ 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
+ return(x*dummy)
+}
+
+jacobian.Hampel <- function (x, ...)
+{
+###########################################
+##
+## R-function: jacobian.Hampel - Jacobian matrix of multivariate
+## analogue of Hampel's psi-function
+## using R-function 'jacobian' of package 'numDeriv'
+## author: Bernhard Spangl, based on work of Paul Gilbert
+## version: 0.1 (2008-02-23)
+##
+###########################################
+
+## Paramters:
+## x ... vector
+
+ jacobian(mvpsiHampel, x, ...)
+}
+
Added: branches/robKalman-Coventry/R/Util.R
===================================================================
--- branches/robKalman-Coventry/R/Util.R (rev 0)
+++ branches/robKalman-Coventry/R/Util.R 2011-08-16 00:05:27 UTC (rev 42)
@@ -0,0 +1,61 @@
+
+### the Euclidean norm
+
+EuclideanNorm <- function(x) {sqrt(colSums(x^2))}
+
+
+### huberizing a vector to length b
+
+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}
+ {SO0 <- S + 1
+ S0 <- S
+ i <- 0
+# print(S0)
+ while( (sum( (SO0 - S0)^2 ) > tol^2) && (i < itmax) )
+ {i <- i + 1
+ S1 <- .getpredCov(S0, F, Q)
+ SO0 <- S0
+ K <- .getKG(S1, Z, .getDelta(S1, Z, V))
+ S0 <- .getcorrCov(S1, K, Z)
+# print(list(S0,i))
+ }
+ S1
+ }
+
+rootMatrix <- function (X)
+{
+###########################################
+##
+## R-function: rootMatrix - computes the unique square root 'A'
+## of matrix 'X', i.e., A%*%A = X
+## former R-function 'root.matrix' of package 'strucchange'
+## author: Bernhard Spangl, based on work of Achim Zeileis
+## version: 0.2 (2008-02-24)
+##
+###########################################
+
+## Paramters:
+## X ... symmetric and positive semidefinite matrix
+
+ if ((ncol(X) == 1) && (nrow(X) == 1))
+ return(list(X.det=X,
+ X.sqrt=matrix(sqrt(X)), X.sqrt.inv=matrix(1/sqrt(X))))
+ else {
+ X.eigen <- eigen(X, symmetric = TRUE)
+ if (any(X.eigen$values < 0))
+ stop("matrix is not positive semidefinite")
+ sqomega <- sqrt(diag(X.eigen$values))
+ sqomega.inv <- diag(1/sqrt(X.eigen$values))
+ V <- X.eigen$vectors
+ X.sqrt <- V %*% sqomega %*% t(V)
+ X.sqrt.inv <- V %*% sqomega.inv %*% t(V)
+ return(list(X.det=prod(X.eigen$values),
+ X.sqrt=X.sqrt, X.sqrt.inv=X.sqrt.inv))
+ }
+}
+
Added: branches/robKalman-Coventry/R/calibrateRLS.R
===================================================================
--- branches/robKalman-Coventry/R/calibrateRLS.R (rev 0)
+++ branches/robKalman-Coventry/R/calibrateRLS.R 2011-08-16 00:05:27 UTC (rev 42)
@@ -0,0 +1,117 @@
+rLScalibrateB <- function(Z, S, V, repl = 100000, b = NULL, eff = NULL, r = NULL,
+ rlow = 0, rup = NULL, upto = 20, IO = FALSE, seed)#
+# 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
+# with /repl/ replicates; 100000 as default works well for me,
+# but might need a change for slower machines
+# upto is an upper bound where to search for zeroes
+{
+
+ 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(!missing(seed)) set.seed(seed)
+
+ dx <- t(mvrnorm(repl, numeric(pd), S))
+ eps <- t(mvrnorm(repl, numeric(qd), V))
+ dy <- Z %*% dx + eps
+ K <- .getKG(S, Z, .getDelta(S, Z, V))
+
+ trS <- sum(diag(.getcorrCov(S, K, Z)))
+
+ dx0 <- K %*% dy
+
+ if(IO){
+ dx0 <- dy - Z %*% 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(ginv(Z))
+ trSb <- sum( (t(dx) - dxw)^2 )/repl
+ trS / trSb
+ }
+ r.b <- function(b){
+ ex <- mean( pmax(no/b - 1, 0) )
+ 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(ginv(Z))
+ trSb <- sum( (t(dX) - dxw)^2 )/repl0
+ 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){
+ (1 - r0)/r0 * sum(pmax(no0 / b - 1, 0))/repl0 - 1
+ }
+ r1 <- r
+ eff1 <- NULL
+ todo.r <- FALSE
+ }
+
+ b <- 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)$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(10^-6, upto*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-Coventry/R/classEKF.R
===================================================================
--- branches/robKalman-Coventry/R/classEKF.R (rev 0)
+++ branches/robKalman-Coventry/R/classEKF.R 2011-08-16 00:05:27 UTC (rev 42)
@@ -0,0 +1,141 @@
+#######################################################
+##
+## classical extended Kalman filter routines
+## author: Bernhard Spangl & Peter Ruckdeschel
+## version: 0.3 (changed: 2011-08-16, 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)
+}
+
+.cEKFinitstep <- 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))
+
+}
+
+.cEKFpredstep <- function (x0, S0, stateEq, mu.v,
+ # F, Q,
+ t, i, additinfofromCorr, ...)
+ # 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
+
+ F <- stateEq$F$fun(t=t, i=i, x0=x0, mu.v=mu.v,
+ exF = stateEq$exo$fun(t=t, i=i, x1=x0,
+ control=stateEq$F$control,
+ additinfofrompast = additinfofromCorr,
+ ...),
+ control = stateEq$F$control,
+ additinfofrompast = additinfofromCorr,...)
+ x1 <- F$x1
+ A <- F$A
+ B <- F$B
+ Q <- stateEq$Q$fun(t=t, i=i, x0=x0,
+ exQ = stateEq$exo$fun(t=t, i=i, x1=x0,
+ control=stateEq$Q$control,
+ additinfofrompast = additinfofromCorr,
+ ...),
+ control = stateEq$Q$control,
+ additinfofrompast = additinfofromCorr,...)
+ Q.m <- Q$Q
+
+ return(list(x1=x1, S1=.getpredCov(S0=S0, A=A, B=B, Q=Q.m),
+ additinfofromPred=additinfofrompast, F=F, Q=Q))
+
+}
+
+.cEKFcorrstep <- function (y, x1, S1, obsEq, mu.eps,
+ # F, Q,
+ t, i, additinfofromPred, ...)
+ #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
+
+ Z <- obsEq$Z$fun(t=t, i=i, x1=x1, mu.eps=mu.eps,
+ exZ = obsEq$exo$fun(t=t, i=i, x1=x1,
+ control=obsEq$Z$control,
+ additinfofrompast = additinfofromPred,
+ ...),
+ control = obsEq$Z$control,
+ additinfofrompast = additinfofromPred,...)
+ yhat <- Z$y
+ C <- Z$Z
+ D <- Z$Z.s
+
+ V <- obsEq$V$fun(t=t, i=i, x1=x1,
+ exV = obsEq$exo$fun(t=t, i=i, x1=x1,
+ control=obsEq$V$control,
+ additinfofrompast = additinfofromPred,
+ ...),
+ control = obsEq$V$control,
+ additinfofrompast = additinfofromPred,...)
+ V.m <- V$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,
+ additinfofromCorr=additinfofrompast, Z=Z, Q=Q))
+
+}
+
+
Added: branches/robKalman-Coventry/R/classKalman.R
===================================================================
--- branches/robKalman-Coventry/R/classKalman.R (rev 0)
+++ branches/robKalman-Coventry/R/classKalman.R 2011-08-16 00:05:27 UTC (rev 42)
@@ -0,0 +1,65 @@
+################################################################
+#classical Kalman filter
+################################################################
+
+# P. Ruckdeschel 10.06.06
+#
+
+## 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
+# +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
+# +state equation x_t = F x_{t-1} + v_t, t>=1
+# +observation equation y_t = Z x_t + e_t , t>=1
+##
+
+
+
+### Notation for the Kalman filter
+
+# +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
+# +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 )^-
+# error covariance S_{t|t} = Cov(x_t-x_{t|t}) = S_{t|t-1} - K_t Z S_{t|t-1}
+
+#########################################################################################
+#
+#Kalman filter routines
+#
+#
+.getDelta <- function(S1, Z, V)# calculates the Cov(Delta y_t) for S1 = S_{t|t-1}, Z, V as above
+{ H <- S1 %*% t(Z)
+ ginv( Z%*%H + V ) }
+
+.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, F, Q)# calculates S_{t|t-1} for S0 = S_{t-1|t-1}, F, Q as above
+{ F %*% S0 %*% t(F) + Q }
+
+##steps for classical Kalman filter (cK)
+.cKinitstep <- function(a, S, ...)
+ {list( x0 = a, S0 = S )}
+
+.cKpredstep <- function(x0, S0, F, Q, ...)
+ {list( x1 = F %*% x0, S1 = .getpredCov(S0, F, Q), Ind=1)}
+
+.cKcorrstep <- function(y, x1, S1, Z, V, ...)
+ {Delta <- .getDelta(S1, Z, V)
+ K <- .getKG(S1, Z, Delta)
+ DeltaY <- y - Z %*% x1
+ x0 <- x1 + K %*% DeltaY
+ S0 <- .getcorrCov(S1, K, Z)
+ list(x0 = x0, K = K, S0 = S0, Delta = Delta, Ind=1, DeltaY = DeltaY)}
Added: branches/robKalman-Coventry/R/classUKF.R
===================================================================
--- branches/robKalman-Coventry/R/classUKF.R (rev 0)
+++ branches/robKalman-Coventry/R/classUKF.R 2011-08-16 00:05:27 UTC (rev 42)
@@ -0,0 +1,149 @@
+#######################################################
+##
+## classical unscented Kalman filter routines
+## author: Bernhard Spangl
+## version: 0.1 (changed: 2011-06-17, created: 2011-06-17)
+##
+#######################################################
+
+.SigmaPoints <- function (x, gamma, Sigma)
+{
+## x ... vector
+## gamma ... scaling parameter
+## Sigma ... covariance matrix
+ sqrtSigma <- t(chol(Sigma))
+ SP1 <- x %o% rep(1, ncol(Sigma)) + gamma * sqrtSigma
+ SP2 <- x %o% rep(1, ncol(Sigma)) - gamma * sqrtSigma
+ cbind(x, SP1, SP2)
+}
+
+.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 - muX %o% rep(1, ncol(X))
+ X <- X * rep(1, nrow(X)) %o% w
+ Y <- Y - muY %o% rep(1, ncol(Y))
+ X %*% Y
+}
+
+.Eval <- function (vec, G, t, ex, control, dim)
+{
+ v1 <- vec[1:dim]
+ v2 <- vec[-(1:dim)]
+ G <- G(t, v1, v2, 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
+
+ alpha <- controlInit$alpha
+ beta <- controlInit$beta
+ kappa <- controlInit$kappa
+
+ pd <- length(a)
+
+ lambda <- alpha^2*(pd + kappa) - pd
+ Wm <- Wc <- rep(1/(2*(pd + lambda)), 2*pd)
+ Wm <- c(lambda/(pd + lambda), Wm)
+ Wc <- c(lambda/(pd + lambda) + (1 + alpha^2 + beta), Wc)
+ gamma <- sqrt(pd + 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}, 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
+
+ gamma <- controlPred$gamma
+ Wm <- controlPred$Wm
+ Wc <- controlPred$Wc
+
+ Qreturn <- Q(t=i, x0=x0, exQ=exQ, control=controlQ)
+ Q <- Qreturn$Q
+
+ X0x <- .SigmaPoints(x=x0, gamma=gamma, Sigma=S0)
+ Xv <- .SigmaPoints(x=v, gamma=gamma, Sigma=Q)
+
+ X1x <- apply(rbind(X0x, Xv), 2, .Eval,
+ G=F, t=i, ex=u, control=controlF, dim=length(x0))
+ x1 <- X1x %*% Wm
+ S1 <- .CovEst(X=X1x, muX=x1, Y=X1x, muY=x1, w=Wc)
+
+ controlPred$X1x <- X1x
+
+ return(list(x1=x1, S1=S1, controlPred=conrolPred))
+
+}
+
+.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
+ X1x <- controlCorr$X1x
+
+ Vretrun <- V(t=i, x1=x1, exV=exV, control=controlV)
+ V <- Vreturn$V
+
+ Xe <- .SigmaPoints(x=eps, gamma=gamma, Sigma=V)
+
+ Y <- apply(rbind(X1x, Xe), 2, .Eval,
+ G=Z, t=i, ex=w, control=controlZ, dim=nrow(X1x))
+ yhat <- Y %*% Wm
+ Sy <- .CovEst(X=Y, mnX=yhat, Y=Y, muY=yhat, w=Wc)
+ Sxy <- .CovEst(X=X1x, muX=x1, Y=Y, muY=yhat, w=Wc)
+ K <- Sxy %*% ginv( Sy )
+ DeltaY <- y - yhat
+ x0 <- x1 + K %*% DeltaY
+ S0 <- S1 - K %*% Sy %*% t(K)
+
+ return(list(x0=x0, K=K, S0=S0, Delta=Sy, DeltaY=DeltaY,
+ controlCorr=controlCorr))
+
+}
+
+
Added: branches/robKalman-Coventry/R/generate.R
===================================================================
--- branches/robKalman-Coventry/R/generate.R (rev 0)
+++ branches/robKalman-Coventry/R/generate.R 2011-08-16 00:05:27 UTC (rev 42)
@@ -0,0 +1,399 @@
+#######################################################
+##
+## generating functions for the (extended) Kalman filter
+## author: Bernhard Spangl & Peter Ruckdeschel
+## version: 0.3 (changed: 2011-08-15, created: 2011-06-09)
+##
+#######################################################
+
+#######################################################
+##
+## Function: 'F_t'
+## Arguments: t, x_{t-1}, v_t,
+## u_{t-1} (exogenous), control
+## Value: x_t, A_t, B_t (Jacobian matrices),
+## original arguments (t, x_{t-1}, v_t, u_{t-1}, control),
+## call,
+## diagnostics
+
+createF <- function (F, ...)
+{
+ UseMethod("createF")
+}
+
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/robkalman -r 42
More information about the Robkalman-commits
mailing list