[Robkalman-commits] r77 - in branches: . robKalman_2015 robKalman_2015/R robKalman_2015/chm robKalman_2015/demo robKalman_2015/inst robKalman_2015/man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Mar 6 11:47:19 CET 2017
Author: ruckdeschel
Date: 2017-03-06 11:47:18 +0100 (Mon, 06 Mar 2017)
New Revision: 77
Added:
branches/robKalman_2015/
branches/robKalman_2015/DESCRIPTION
branches/robKalman_2015/NAMESPACE
branches/robKalman_2015/R/
branches/robKalman_2015/R/0AllClass.R
branches/robKalman_2015/R/ACMfilt.R
branches/robKalman_2015/R/ACMfilter.R
branches/robKalman_2015/R/AllGenerics.R
branches/robKalman_2015/R/AllInitialize.R
branches/robKalman_2015/R/AllPlot.R
branches/robKalman_2015/R/AllShow.R
branches/robKalman_2015/R/Psi.R
branches/robKalman_2015/R/Util.R
branches/robKalman_2015/R/calibrateRLS.R
branches/robKalman_2015/R/classEKF.R
branches/robKalman_2015/R/classKalman.R
branches/robKalman_2015/R/classUKF.R
branches/robKalman_2015/R/generate.R
branches/robKalman_2015/R/mACMfilter.R
branches/robKalman_2015/R/mACMinternal.R
branches/robKalman_2015/R/rLSfilter.R
branches/robKalman_2015/R/recEFilter.R
branches/robKalman_2015/R/recFilterInternal.R
branches/robKalman_2015/R/recFilterWrapper.R
branches/robKalman_2015/R/simulateSScont.R
branches/robKalman_2015/chm/
branches/robKalman_2015/chm/00Index.html
branches/robKalman_2015/chm/0robKalman-package.html
branches/robKalman_2015/chm/ACMfilt.html
branches/robKalman_2015/chm/Logo.jpg
branches/robKalman_2015/chm/Rchm.css
branches/robKalman_2015/chm/arGM.html
branches/robKalman_2015/chm/calibrateRLS.html
branches/robKalman_2015/chm/internalACM.html
branches/robKalman_2015/chm/internalKalman.html
branches/robKalman_2015/chm/internalarGM.html
branches/robKalman_2015/chm/internalpsi.html
branches/robKalman_2015/chm/internalrLS.html
branches/robKalman_2015/chm/recFilter.html
branches/robKalman_2015/chm/robKalman.chm
branches/robKalman_2015/chm/robKalman.hhp
branches/robKalman_2015/chm/robKalman.toc
branches/robKalman_2015/chm/simulateSScont.html
branches/robKalman_2015/chm/util.html
branches/robKalman_2015/demo/
branches/robKalman_2015/demo/00Index
branches/robKalman_2015/demo/ACMdemo.R
branches/robKalman_2015/demo/rLSdemo.R
branches/robKalman_2015/inst/
branches/robKalman_2015/inst/NEWS
branches/robKalman_2015/man/
branches/robKalman_2015/man/0robKalman-package.Rd
branches/robKalman_2015/man/ACMfilt.Rd
branches/robKalman_2015/man/arGM.Rd
branches/robKalman_2015/man/calibrateRLS.Rd
branches/robKalman_2015/man/internalACM.Rd
branches/robKalman_2015/man/internalKalman.Rd
branches/robKalman_2015/man/internalarGM.Rd
branches/robKalman_2015/man/internalpsi.Rd
branches/robKalman_2015/man/internalrLS.Rd
branches/robKalman_2015/man/recFilter.Rd
branches/robKalman_2015/man/simulateSScont.Rd
branches/robKalman_2015/man/util.Rd
Log:
committed some changes done offline till Nov. 2015
Added: branches/robKalman_2015/DESCRIPTION
===================================================================
--- branches/robKalman_2015/DESCRIPTION (rev 0)
+++ branches/robKalman_2015/DESCRIPTION 2017-03-06 10:47:18 UTC (rev 77)
@@ -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 uni-oldenburg.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_2015/NAMESPACE
===================================================================
--- branches/robKalman_2015/NAMESPACE (rev 0)
+++ branches/robKalman_2015/NAMESPACE 2017-03-06 10:47:18 UTC (rev 77)
@@ -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_2015/R/0AllClass.R
===================================================================
--- branches/robKalman_2015/R/0AllClass.R (rev 0)
+++ branches/robKalman_2015/R/0AllClass.R 2017-03-06 10:47:18 UTC (rev 77)
@@ -0,0 +1,14 @@
+.onLoad <- function(lib, pkg){
+ require("methods", character = TRUE, quietly = TRUE)
+
+}
+
+
+.onAttach <- function(library, pkg)
+{
+buildStartupMessage(pkg="robKalman", library=library, packageHelp=TRUE #,
+ #MANUAL=""
+ )
+ invisible()
+}
+
Added: branches/robKalman_2015/R/ACMfilt.R
===================================================================
--- branches/robKalman_2015/R/ACMfilt.R (rev 0)
+++ branches/robKalman_2015/R/ACMfilt.R 2017-03-06 10:47:18 UTC (rev 77)
@@ -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_2015/R/ACMfilter.R
===================================================================
--- branches/robKalman_2015/R/ACMfilter.R (rev 0)
+++ branches/robKalman_2015/R/ACMfilter.R 2017-03-06 10:47:18 UTC (rev 77)
@@ -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_2015/R/AllGenerics.R
===================================================================
--- branches/robKalman_2015/R/AllGenerics.R (rev 0)
+++ branches/robKalman_2015/R/AllGenerics.R 2017-03-06 10:47:18 UTC (rev 77)
@@ -0,0 +1,123 @@
+if(!isGeneric("solve")){
+ setGeneric("solve", function(a,b,...) standardGeneric("solve"))
+}
+############################################################################
+# Access methods
+############################################################################
+
+if(!isGeneric("name"))
+ setGeneric("name", function(object) standardGeneric("name"))
+
+if(!isGeneric("getp"))
+ setGeneric("getp", function(object) standardGeneric("getp"))
+if(!isGeneric("getq"))
+ setGeneric("getq", function(object) standardGeneric("getq"))
+
+if(!isGeneric("getF"))
+ setGeneric("getF", function(object) standardGeneric("getF"))
+if(!isGeneric("getZ"))
+ setGeneric("getZ", function(object) standardGeneric("getZ"))
+if(!isGeneric("getQ"))
+ setGeneric("getQ", function(object) standardGeneric("getQ"))
+if(!isGeneric("getV"))
+ setGeneric("getV", function(object) standardGeneric("getV"))
+if(!isGeneric("geta"))
+ setGeneric("geta", function(object) standardGeneric("geta"))
+if(!isGeneric("getS"))
+ setGeneric("getS", function(object) standardGeneric("getS"))
+if(!isGeneric("time"))
+ setGeneric("time", function(x,...) standardGeneric("time"))
+
+if(!isGeneric("SSM"))
+ setGeneric("SSM", function(object, ...) standardGeneric("SSM"))
+if(!isGeneric("Y"))
+ setGeneric("Y", function(object, ...) standardGeneric("Y"))
+if(!isGeneric("X.filtered"))
+ setGeneric("X.filtered", function(object, ...) standardGeneric("X.filtered"))
+if(!isGeneric("X.predicted"))
+ setGeneric("X.predicted", function(object, ...) standardGeneric("X.predicted"))
+if(!isGeneric("Cov.filtered"))
+ setGeneric("Cov.filtered", function(object, ...) standardGeneric("Cov.filtered"))
+if(!isGeneric("Cov.predicted"))
+ setGeneric("Cov.predicted", function(object, ...) standardGeneric("Cov.predicted"))
+if(!isGeneric("Kalman.Gain"))
+ setGeneric("Kalman.Gain", function(object, ...) standardGeneric("Kalman.Gain"))
+if(!isGeneric("X.rob.filtered"))
+ setGeneric("X.rob.filtered", function(object, ...) standardGeneric("X.rob.filtered"))
+if(!isGeneric("X.rob.predicted"))
+ setGeneric("X.rob.predicted", function(object, ...) standardGeneric("X.rob.predicted"))
+if(!isGeneric("Cov.rob.filtered"))
+ setGeneric("Cov.rob.filtered", function(object, ...) standardGeneric("Cov.rob.filtered"))
+if(!isGeneric("Cov.rob.predicted"))
+ setGeneric("Cov.rob.predicted", function(object, ...) standardGeneric("Cov.rob.predicted"))
+if(!isGeneric("Kalman.rob.Gain"))
+ setGeneric("Kalman.rob.Gain", function(object, ...) standardGeneric("Kalman.rob.Gain"))
+if(!isGeneric("rob.correction.ctrl"))
+ setGeneric("rob.correction.ctrl", function(object, ...) standardGeneric("rob.correction.ctrl"))
+if(!isGeneric("rob.prediction.ctrl"))
+ setGeneric("rob.prediction.ctrl", function(object, ...) standardGeneric("rob.prediction.ctrl"))
+if(!isGeneric("IndIO"))
+ setGeneric("IndIO", function(object, ...) standardGeneric("IndIO"))
+if(!isGeneric("IndAO"))
+ setGeneric("IndAO", function(object, ...) standardGeneric("IndAO"))
+if(!isGeneric("nsim"))
+ setGeneric("nsim", function(object, ...) standardGeneric("nsim"))
+if(!isGeneric("RNGstate"))
+ setGeneric("RNGstate", function(object, ...) standardGeneric("RNGstate"))
+if(!isGeneric("Cov.rob.filtered.sim"))
+ setGeneric("Cov.rob.filtered.sim", function(object, ...) standardGeneric("Cov.rob.filtered.sim"))
+if(!isGeneric("Cov.rob.predicted.sim"))
+ setGeneric("Cov.rob.predicted.sim", function(object, ...) standardGeneric("Cov.rob.predicted.sim"))
+if(!isGeneric("init"))
+ setGeneric("init", function(object, ...) standardGeneric("init"))
+if(!isGeneric("predict"))
+ setGeneric("predict", function(object, ...) standardGeneric("predict"))
+if(!isGeneric("correct"))
+ setGeneric("correct", function(object, ...) standardGeneric("correct"))
+if(!isGeneric("init.rob"))
+ setGeneric("init.rob", function(object, ...) standardGeneric("init.rob"))
+if(!isGeneric("name.rob"))
+ setGeneric("name.rob", function(object, ...) standardGeneric("name.rob"))
+if(!isGeneric("predict.rob"))
+ setGeneric("predict.rob", function(object, ...) standardGeneric("predict.rob"))
+if(!isGeneric("correct.rob"))
+ setGeneric("correct.rob", function(object, ...) standardGeneric("correct.rob"))
+if(!isGeneric("controls"))
+ setGeneric("controls", function(object, ...) standardGeneric("controls"))
+
+############################################################################
+# Replacement methods
+############################################################################
+
+if(!isGeneric("name<-"))
+ setGeneric("name<-",
+ function(object, value) standardGeneric("name<-"))
+
+if(!isGeneric("setp<-"))
+ setGeneric("setp<-", function(object, value) standardGeneric("setp<-"))
+if(!isGeneric("setq<-"))
+ setGeneric("setq<-", function(object, value) standardGeneric("setq<-"))
+
+if(!isGeneric("setF<-"))
+ setGeneric("setF<-", function(object, value) standardGeneric("setF<-"))
+if(!isGeneric("setZ<-"))
+ setGeneric("setZ<-", function(object, value) standardGeneric("setZ<-"))
+if(!isGeneric("setQ<-"))
+ setGeneric("setQ<-", function(object, value) standardGeneric("setQ<-"))
+if(!isGeneric("setV<-"))
+ setGeneric("setV<-", function(object, value) standardGeneric("setV<-"))
+if(!isGeneric("seta<-"))
+ setGeneric("seta<-", function(object, value) standardGeneric("seta<-"))
+if(!isGeneric("setS<-"))
+ setGeneric("setS<-", function(object, value) standardGeneric("setS<-"))
+if(!isGeneric("time<-"))
+ setGeneric("time<-", function(x, value) standardGeneric("time<-"))
+
+if(!isGeneric(".make.project"))
+setGeneric(".make.project",function(object, ...) standardGeneric(".make.project"))
+
+if(!isGeneric("kalman"))
+setGeneric("kalman",function(smooth, ...) standardGeneric("kalman"))
+
+if(!isGeneric("kalmanRob"))
+setGeneric("kalmanRob",function(method, smooth, ...) standardGeneric("kalmanRob"))
Added: branches/robKalman_2015/R/AllInitialize.R
===================================================================
--- branches/robKalman_2015/R/AllInitialize.R (rev 0)
+++ branches/robKalman_2015/R/AllInitialize.R 2017-03-06 10:47:18 UTC (rev 77)
@@ -0,0 +1,3 @@
+#### as to whether to use Generating functions or to use initialize methods:
+#### http://tolstoy.newcastle.edu.au/R/e2/devel/07/01/1976.html
+
Added: branches/robKalman_2015/R/AllPlot.R
===================================================================
--- branches/robKalman_2015/R/AllPlot.R (rev 0)
+++ branches/robKalman_2015/R/AllPlot.R 2017-03-06 10:47:18 UTC (rev 77)
@@ -0,0 +1,7 @@
+#setMethod("plot", "ParamFamily",
+# function(x,y=NULL,...){
+# e1 <- x at distribution
+# if(!is(e1, "UnivariateDistribution")) stop("not yet implemented")
+#
+# plot(e1)
+# })
Added: branches/robKalman_2015/R/AllShow.R
===================================================================
--- branches/robKalman_2015/R/AllShow.R (rev 0)
+++ branches/robKalman_2015/R/AllShow.R 2017-03-06 10:47:18 UTC (rev 77)
@@ -0,0 +1,161 @@
+#setMethod("show", "ParamFamParameter",
+# function(object){
+# cat(paste("An object of class", dQuote(class(object)), "\n"))
+# cat("name:\t", object at name, "\n")
+# cat("main:\t", object at main, "\n")
+# if(!is.null(object at nuisance))
+# cat("nuisance:\t", object at nuisance, "\n")
+# if(!identical(all.equal(object at trafo, diag(length(object)),
+# tolerance = .Machine$double.eps^0.5), TRUE)){
+# cat("trafo:\n")
+# print(object at trafo)
+# }
+# })
+#setMethod("show", "Symmetry",
+# function(object){
+# cat("type of symmetry:\t", object at type, "\n")
+# if(!is.null(object at SymmCenter))
+# cat("center of symmetry:\n")
+# print(object at SymmCenter)
+# })
+#setMethod("show", "ParamFamily",
+# function(object){
+# cat(paste("An object of class", dQuote(class(object)), "\n"))
+# cat("### name:\t", object at name, "\n")
+# cat("\n### distribution:\t")
+# print(object at distribution)
+# cat("\n### param:\t")
+# show(object at param)
+# if(length(object at props) != 0){
+# cat("\n### props:\n")
+# show(object at props)
+# }
+# })
+#setMethod("show", "Neighborhood",
+# function(object){
+# cat(paste("An object of class", dQuote(class(object)), "\n"))
+# cat("type:\t", object at type, "\n")
+# cat("radius:\t", object at radius, "\n")
+# })
+#setMethod("show", "FixRobModel",
+# function(object){
+# cat(paste("An object of class", dQuote(class(object)), "\n"))
+# cat("###### center:\t")
+# show(object at center)
+# cat("\n###### neighborhood:\t")
+# show(object at neighbor)
+# })
+#setMethod("show", "InfRobModel",
+# function(object){
+# cat(paste("An object of class", dQuote(class(object)), "\n"))
+# cat("###### center:\t")
+# show(object at center)
+# cat("\n###### neighborhood:\t")
+# show(object at neighbor)
+# })
+#setMethod("show", "RiskType",
+# function(object){
+# cat(paste("An object of class", dQuote(class(object)), "\n"))
+# cat("risk type:\t", object at type, "\n")
+# })
+#setMethod("show", "asUnOvShoot",
+# function(object){
+# cat(paste("An object of class", dQuote(class(object)), "\n"))
+# cat("risk type:\t", object at type, "\n")
+# cat("width:\t", object at width, "\n")
+# })
+#setMethod("show", "asHampel",
+# function(object){
+# cat(paste("An object of class", dQuote(class(object)), "\n"))
+# cat("risk type:\t", object at type, "\n")
+# cat("bound:\t", object at bound, "\n")
+# })
+#setMethod("show", "fiUnOvShoot",
+# function(object){
+# cat(paste("An object of class", dQuote(class(object)), "\n"))
+# cat("risk type:\t", object at type, "\n")
+# cat("width:\t", object at width, "\n")
+# })
+#setMethod("show", "fiHampel",
+# function(object){
+# cat(paste("An object of class", dQuote(class(object)), "\n"))
+# cat("risk type:\t", object at type, "\n")
+# cat("bound:\t", object at bound, "\n")
+# })
+#setMethod("show", "InfluenceCurve",
+# function(object){
+# cat(paste("An object of class", dQuote(class(object)), "\n"))
+# cat("### name:\t", object at name, "\n")
+# cat("\n### 'Curve':\t")
+# show(object at Curve)
+## cat("\n### Risks:\n")
+## print(object at Risks)
+# cat("\n### Infos:\n")
+# print(object at Infos)
+# })
+#setMethod("show", "IC",
+# function(object){
+# cat(paste("An object of class", dQuote(class(object)), "\n"))
+# cat("### name:\t", object at name, "\n")
+# L2Fam <- eval(object at CallL2Fam)
+# cat("### L2-differentiable parametric family:\t", L2Fam at name, "\n")
+# cat("\n### 'Curve':\t")
+# show(object at Curve)
+## cat("\n### Risks:\n")
+## print(object at Risks)
+# cat("\n### Infos:\n")
+# print(object at Infos)
+# })
+#setMethod("show", "ContIC",
+# function(object){
+# cat(paste("An object of class", dQuote(class(object)), "\n"))
+# cat("### name:\t", object at name, "\n")
+# L2Fam <- eval(object at CallL2Fam)
+# cat("\n### L2-differentiable parametric family:\t", L2Fam at name, "\n")
+# cat("### param:\t")
+# show(L2Fam at param)
+# cat("\n### neighborhood radius:\t", object at neighborRadius, "\n")
+## cat("\n### 'Curve':\t")
+## show(object at Curve)
+# cat("\n### clip:\t")
+# show(object at clip)
+# cat("### cent:\t")
+# show(object at cent)
+# cat("### stand:\n")
+# show(object at stand)
+# if(!is.null(object at lowerCase)){
+# cat("### lowerCase:\t")
+# show(object at lowerCase)
+# }
+## cat("\n### Risks:\n")
+## show(object at Risks)
+# cat("\n### Infos:\n")
+# show(object at Infos)
+# })
+#setMethod("show", "TotalVarIC",
+# function(object){
+# cat(paste("An object of class", dQuote(class(object)), "\n"))
+# cat("### name:\t", object at name, "\n")
+# L2Fam <- eval(object at CallL2Fam)
+# cat("\n### L2-differentiable parametric family:\t", L2Fam at name, "\n")
+# cat("### param:\t")
+# show(L2Fam at param)
+# cat("\n### neighborhood radius:\t", object at neighborRadius, "\n")
+## cat("\n### 'Curve':\t")
+## show(object at Curve)
+# cat("\n### clipLo:\t")
+# show(object at clipLo)
+# cat("### clipUp:\t")
+# show(object at clipUp)
+# cat("### stand:\n")
+# show(object at stand)
+## cat("\n### Risks:\n")
+## show(object at Risks)
+# cat("\n### Infos:\n")
+# show(object at Infos)
+# })
+#
+#
+#
+#
+#
Added: branches/robKalman_2015/R/Psi.R
===================================================================
--- branches/robKalman_2015/R/Psi.R (rev 0)
+++ branches/robKalman_2015/R/Psi.R 2017-03-06 10:47:18 UTC (rev 77)
@@ -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_2015/R/Util.R
===================================================================
--- branches/robKalman_2015/R/Util.R (rev 0)
+++ branches/robKalman_2015/R/Util.R 2017-03-06 10:47:18 UTC (rev 77)
@@ -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_2015/R/calibrateRLS.R
===================================================================
--- branches/robKalman_2015/R/calibrateRLS.R (rev 0)
+++ branches/robKalman_2015/R/calibrateRLS.R 2017-03-06 10:47:18 UTC (rev 77)
@@ -0,0 +1,126 @@
+rLScalibrateB <- function(Z, S, V, repl = 100000, b = NULL, eff = NULL, r = NULL,
+ rlow = 0, rup = NULL, upto = 20, IO = FALSE, seed = NULL)#
+# 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(!is.null(seed)) set.seed(seed)
+
+ dx <- t(mvrnorm(repl, numeric(pd), S))
+ eps <- t(mvrnorm(repl, numeric(qd), V))
+ dy <- Z %*% dx + eps
+ gD <- .getDelta(S, Z, V)
+ K <- .getKG(S, Z, gD)
+
+ trS <- sum(diag(.getcorrCov(S, K, Z)))
+
+ dx0 <- K %*% dy
+
+ 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 - dxq)^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(Zs)
+ trSb <- sum( (t(dX) - dxw - dxq)^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_2015/R/classEKF.R
===================================================================
--- branches/robKalman_2015/R/classEKF.R (rev 0)
+++ branches/robKalman_2015/R/classEKF.R 2017-03-06 10:47:18 UTC (rev 77)
@@ -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)
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/robkalman -r 77
More information about the Robkalman-commits
mailing list