[Robkalman-commits] r25 - in pkg/robKalman: . R chm demo inst man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Mar 18 16:41:02 CET 2009
Author: ruckdeschel
Date: 2009-03-18 16:41:02 +0100 (Wed, 18 Mar 2009)
New Revision: 25
Added:
pkg/robKalman/DESCRIPTION
pkg/robKalman/NAMESPACE
pkg/robKalman/R/
pkg/robKalman/R/0AllClass.R
pkg/robKalman/R/ACMfilt.R
pkg/robKalman/R/ACMfilter.R
pkg/robKalman/R/AllGeneric.R
pkg/robKalman/R/AllGenerics.R
pkg/robKalman/R/AllInitialize.R
pkg/robKalman/R/AllPlot.R
pkg/robKalman/R/AllShow.R
pkg/robKalman/R/Psi.R
pkg/robKalman/R/Util.R
pkg/robKalman/R/arGM.R
pkg/robKalman/R/arGMinternal.R
pkg/robKalman/R/calibrateRLS.R
pkg/robKalman/R/classKalman.R
pkg/robKalman/R/mACMfilter.R
pkg/robKalman/R/mACMinternal.R
pkg/robKalman/R/rLSfilter.R
pkg/robKalman/R/recFilter.R
pkg/robKalman/R/simulateSScont.R
pkg/robKalman/chm/
pkg/robKalman/chm/00Index.html
pkg/robKalman/chm/0robKalman-package.html
pkg/robKalman/chm/ACMfilt.html
pkg/robKalman/chm/Logo.jpg
pkg/robKalman/chm/Rchm.css
pkg/robKalman/chm/arGM.html
pkg/robKalman/chm/calibrateRLS.html
pkg/robKalman/chm/internalACM.html
pkg/robKalman/chm/internalKalman.html
pkg/robKalman/chm/internalarGM.html
pkg/robKalman/chm/internalpsi.html
pkg/robKalman/chm/internalrLS.html
pkg/robKalman/chm/recFilter.html
pkg/robKalman/chm/robKalman.chm
pkg/robKalman/chm/robKalman.hhp
pkg/robKalman/chm/robKalman.toc
pkg/robKalman/chm/simulateSScont.html
pkg/robKalman/chm/util.html
pkg/robKalman/demo/
pkg/robKalman/demo/00Index
pkg/robKalman/demo/ACMdemo.R
pkg/robKalman/demo/rLSdemo.R
pkg/robKalman/inst/
pkg/robKalman/inst/NEWS
pkg/robKalman/man/
pkg/robKalman/man/0robKalman-package.Rd
pkg/robKalman/man/ACMfilt.Rd
pkg/robKalman/man/arGM.Rd
pkg/robKalman/man/calibrateRLS.Rd
pkg/robKalman/man/internalACM.Rd
pkg/robKalman/man/internalKalman.Rd
pkg/robKalman/man/internalarGM.Rd
pkg/robKalman/man/internalpsi.Rd
pkg/robKalman/man/internalrLS.Rd
pkg/robKalman/man/recFilter.Rd
pkg/robKalman/man/simulateSScont.Rd
pkg/robKalman/man/util.Rd
Log:
yet another restructuring:
copied R-package robKalman from trunc, revision 12 into trunc (step 2)
should be an installable version now...
Added: pkg/robKalman/DESCRIPTION
===================================================================
--- pkg/robKalman/DESCRIPTION (rev 0)
+++ pkg/robKalman/DESCRIPTION 2009-03-18 15:41:02 UTC (rev 25)
@@ -0,0 +1,13 @@
+Package: robKalman
+Version: 0.2
+Date: 2008-07-28
+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, dse1, dse2, MASS, limma, robustbase, numDeriv
+Imports: stats, MASS
+SaveImage: no
+LazyLoad: yes
+License: GPL (version 2 or later)
+URL: https://r-forge.r-project.org/projects/robkalman
Added: pkg/robKalman/NAMESPACE
===================================================================
--- pkg/robKalman/NAMESPACE (rev 0)
+++ pkg/robKalman/NAMESPACE 2009-03-18 15:41:02 UTC (rev 25)
@@ -0,0 +1,44 @@
+import("methods")
+import("stats")
+import("startupmsg")
+
+exportClasses("PosSemDefSymmMatrix","PosDefSymmMatrix")
+exportClasses("ArrayOrMatrix",
+ "Hyperparamtype",
+ "sHyperparamtype")
+exportClasses("SSM","TimeInvariantSSM")
+exportClasses("recFilter","robrecFilter")
+exportClasses("RecFiltControl","KalmanControl","robrecControl")
+exportClasses("SSMDistribution.f","SSMellDistribution.f",
+ "SSMConvDistribution.f")
+exportClasses("SSMwithDistribution","SSMwithConvDistribution")
+exportClasses("SSMDistr")
+exportClasses("SSMsimulation","SSMcontSimulation")
+exportMethods("getp","setp<-","getq","setq<-",
+ "getF","setF<-","getZ","setZ<-",
+ "getV","setV<-","getQ","setQ<-",
+ "geta","seta<-","getS","setS<-",
+ "time","time<-","name","name<-")
+exportMethods("SSM", "Y", "X.filtered", "X.predicted",
+ "Cov.filtered", "Cov.predicted","Kalman.Gain",
+ "X.rob.filtered", "X.rob.predicted",
+ "Cov.rob.filtered", "Cov.rob.predicted","Kalman.rob.Gain")
+exportMethods("rob.correction.ctrl", "rob.prediction.ctrl")
+exportMethods("IndIO", "IndAO", "nsim", "RNGstate")
+exportMethods("Cov.rob.filtered.sim", "Cov.rob.predicted.sim")
+exportMethods("init", "predict", "correct")
+exportMethods("init.rob", "predict.rob", "correct.rob", "name.rob")
+exportMethods("controls")
+exportMethods(".make.project", "kalman", "kalmanRob")
+exportMethods("solve", "simulate")
+export("TI.SSM", "makeArrayRepresentation")
+export("SSMellDistribution.f", "SSMContDistribution.f",
+ "SSMConvDistribution.f", "SSMwithDistribution",
+ "SSMwithConvDistribution")
+export("rcvcont")
+export("KalmanControl", "ACMControl", "rLSControl")
+export("generateRecFilter", "generateRobRecFilter")
+export("ACMfilt", "ACMfilter", "arGM", "Euclidnorm",
+ "simulateState", "simulateObs", "rcvmvnorm", "Huberize",
+ "rLScalibrateB", "limitS", "rLSFilter", "KalmanFilter",
+ "recursiveFilter")
Added: pkg/robKalman/R/0AllClass.R
===================================================================
--- pkg/robKalman/R/0AllClass.R (rev 0)
+++ pkg/robKalman/R/0AllClass.R 2009-03-18 15:41:02 UTC (rev 25)
@@ -0,0 +1,335 @@
+################ general package preparation code, nothing to do with classes:
+
+.onLoad <- function(lib, pkg){
+ require("methods", character = TRUE, quietly = TRUE)
+
+}
+
+
+.onAttach <- function(library, pkg)
+{
+buildStartupMessage(pkg="robKalman", library=library, packageHelp=TRUE #,
+ #MANUAL=""
+ )
+ invisible()
+}
+
+#.onUnload <- function(libpath)
+#{
+# library.dynam.unload("distrEx", libpath)
+#}
+#
+
+
+################ Matrix class
+
+#### Code borrowed from package distrMod
+
+if (!(isClass("PosSemDefSymmMatrix")))
+setClass("PosSemDefSymmMatrix", contains = "matrix",
+ prototype = prototype(matrix(1)),
+ validity = function(object){
+ if(nrow(object) != ncol(object))
+ stop("no square matrix")
+ if(any(!is.finite(object)))
+ stop("inifinite or missing values in matrix")
+ if(!isTRUE(all.equal(object, t(object), .Machine$double.eps^0.5)))
+ stop("matrix is not symmetric")
+ if(!all(eigen(object)$values > -100*.Machine$double.eps))
+ stop("matrix is (numerically) not positive semi - definite")
+ return(TRUE)
+ })
+
+## positive definite, symmetric matrices with finite entries
+if (!(isClass("PosDefSymmMatrix")))
+setClass("PosDefSymmMatrix", contains = "PosSemDefSymmMatrix",
+ validity = function(object){
+ if(!all(eigen(object)$values > 100*.Machine$double.eps))
+ stop("matrix is (numerically) not positive definite")
+ valid <- getValidity(getClass("PosSemDefSymmMatrix"))
+ valid(as(object, "PosSemDefSymmMatrix"))
+ return(TRUE)
+ })
+
+
+
+#
+## register zoo as "S4"-class
+setOldClass("zoo")
+#
+
+### infra-structure classes / class unions
+
+setClassUnion("IntegerOrNULL", c("integer", "NULL"))
+setClassUnion("ArrayOrNULL", c("array", "NULL"))
+setClassUnion("ArrayOrMatrix", c("array", "matrix"))
+setClassUnion("Hyperparamtype",
+ c("NULL","ArrayOrMatrix", "OptionalFunction"))
+setClassUnion("sHyperparamtype",
+ c("Hyperparamtype", "numeric"))
+setClassUnion("MatrixOrLogical", c("logical", "matrix"))
+
+
+###############################################################################
+#
+# State Space Model classes (SSMs)
+#
+###############################################################################
+
+
+# class SSM --- State space model
+setClass("SSM",
+ representation = representation(
+ name = "character", ## name of the ssm
+ F = "Hyperparamtype", ## transition matrix/ces or NULL
+ Z = "Hyperparamtype", ## observation matrix/ces or NULL
+ Q = "Hyperparamtype", ## innovation covariance or NULL
+ V = "Hyperparamtype", ## observation error covariance or NULL
+ p = "numeric", ## state dimension
+ q = "numeric", ## observation dimension
+ a = "numeric", ## mean value of starting state
+ S = "Hyperparamtype", ## variance of starting state
+ time = "zoo"), ## time index
+ prototype = prototype(name = gettext("a state space"),
+ F = NULL,
+ Z = NULL,
+ Q = NULL,
+ V = NULL,
+ p = 1,
+ q = 1,
+ a = 0,
+ S = NULL,
+ time = zoo(1)),
+ )
+
+# class TimeInvariantSSM
+setClass("TimeInvariantSSM",
+ prototype = prototype(name = gettext("a time-invariant state space"),
+ F = matrix(1),
+ Z = matrix(1),
+ Q = matrix(1),
+ V = matrix(1),
+ p = 1,
+ q = 1,
+ a = 0,
+ S = matrix(1)),
+ validity = function(object){
+ if(!is.matrix(object at F)|!is.matrix(object at Z)|
+ !is.matrix(object at Q)|!is.matrix(object at V)|
+ !is.matrix(object at S))
+ stop("Hyperparameters have to be matrices")
+ return(TRUE)
+ },
+ contains = "SSM")
+
+###############################################################################
+#
+# Filter classes
+#
+###############################################################################
+
+setClass("recFilter", representation(name = "character",
+ SSM = "SSM",
+ Y = "array",
+ X.filtered = "array",
+ X.predicted = "array",
+ Cov.filtered = "array",
+ Cov.predicted = "array",
+ Kalman.Gain = "array",
+ time = "zoo"),
+ prototype = prototype(name="classical Kalman Filter",
+ SSM = new("TimeInvariantSSM"),
+ Y = array(1,dim=c(1,1,1)),
+ X.filtered = array(1,dim=c(1,1,1)),
+ X.predicted = array(1,dim=c(1,1,1)),
+ Cov.filtered = array(1,dim = c(1,1,1)),
+ Cov.predicted = array(1,dim = c(1,1,1)),
+ Kalman.Gain = array(1,dim = c(1,1,1)),
+ time = zoo(1))
+
+ )
+
+
+setClass("robrecFilter", representation(
+ name.rob = "character",
+ X.rob.filtered = "array",
+ X.rob.predicted = "array",
+ Cov.rob.filtered = "array",
+ Cov.rob.predicted = "array",
+ Kalman.rob.Gain = "array",
+ IndIO = "MatrixOrLogical",
+ IndAO = "MatrixOrLogical",
+ rob.correction.ctrl = "list",
+ rob.prediction.ctrl = "list",
+ nsim = "numeric",
+ RNGstate = "IntegerOrNULL",
+ Cov.rob.filtered.sim = "ArrayOrNULL",
+ Cov.rob.predicted.sim = "ArrayOrNULL"
+ ),
+ prototype = prototype(
+ name="rLS Filter",
+ X.rob.filtered = array(1,dim=c(1,1,1)),
+ X.rob.predicted = array(1,dim=c(1,1,1)),
+ Cov.rob.filtered = array(1,dim = c(1,1,1)),
+ Cov.rob.predicted = array(1,dim = c(1,1,1)),
+ Kalman.rob.Gain = array(1,dim = c(1,1,1)),
+ IndIO = FALSE,
+ IndAO = FALSE,
+ nsim = 0,
+ RNGstate = as.integer(0),
+ rob.correction.ctrl = list(NULL),
+ rob.prediction.ctrl = list(NULL),
+ Cov.rob.filtered.sim = array(1,dim = c(1,1,1)),
+ Cov.rob.predicted.sim = array(1,dim = c(1,1,1))),
+ contains = "recFilter")
+
+###############################################################################
+#
+# Control classes
+#
+###############################################################################
+#setClass("ACMcontrol",representation())
+#setClass("rLScontrol",representation())
+
+
+setClass("RecFiltControl", representation(
+ name = "character",
+ init = "function",
+ predict = "function",
+ correct = "function"),
+ contains = "VIRTUAL")
+
+setClass("KalmanControl",
+ prototype = prototype(
+ name = "classical Kalman Filter",
+ init = .cKinitstep,
+ predict = .cKpredstep,
+ correct = .cKcorrstep),
+ contains="RecFiltControl"
+ )
+
+setClass("robrecControl", representation(
+ controls = "list",
+ name.rob = "character",
+ init.rob = "function",
+ predict.rob = "function",
+ correct.rob = "function"),
+ prototype = prototype(
+ name.rob = "rLS Filter",
+ init.rob = .cKinitstep,
+ predict.rob = .cKpredstep,
+ correct.rob = .rLScorrstep,
+ controls = list(b = 2, norm = EuclideanNorm)
+ ),
+ contains="RecFiltControl"
+ )
+# = paste(gettext(
+# "Control set and init, prediction and"
+# ),gettext(
+# "correction step for the classical Kalman Filter\n"
+# ),gettext(
+# "and the rLS Filter"
+# )),
+
+###############################################################################
+#
+# multivariate Distribution classes
+#
+###############################################################################
+
+setClass("SSMDistribution.f", representation(
+ r.init = "function",
+ r.innov = "function",
+ r.obs = "function"),
+ prototype = prototype(r.init = mvrnorm,
+ r.innov = mvrnorm,
+ r.obs = mvrnorm))
+
+setClass("SSMellDistribution.f", representation(
+ m.init = "sHyperparamtype",
+ S.init = "Hyperparamtype",
+ m.innov = "sHyperparamtype",
+ S.innov = "Hyperparamtype",
+ m.obs = "sHyperparamtype",
+ S.obs = "Hyperparamtype"),
+ prototype = prototype(m.init = 0, S.init = matrix(1),
+ m.innov = 0, S.innov = matrix(1),
+ m.obs = 0, S.obs = matrix(1)),
+ contains = "SSMDistribution.f")
+
+setClass("SSMConvDistribution.f", representation(
+ ideal = "SSMDistribution.f",
+ cont = "SSMellDistribution.f",
+ r.IO = "numeric",
+ r.AO = "numeric"),
+ prototype = prototype(ideal = new("SSMellDistribution.f"),
+ cont = new("SSMellDistribution.f"),
+ r.IO = 0, r.AO = 0),
+ validity = function(object){
+ if (object at r.IO<0||object at r.AO<0||object at r.IO>1||object at r.AO>1)
+ stop("Radii must be between 0 and 1")
+ return(TRUE)})
+
+
+###############################################################################
+#
+# SSM + Distribution classes
+#
+###############################################################################
+setClass("SSMwithDistribution", representation(
+ SSM = "SSM",
+ Distribution = "SSMDistribution.f"),
+ prototype = prototype(SSM = new("SSM"),
+ Distribution = SSMDistribution.f(new("SSM"))))
+setClass("SSMwithConvDistribution", representation(
+ SSM = "SSM",
+ Distribution = "SSMConvDistribution.f"),
+ prototype = prototype(SSM = new("SSM"),
+ Distribution = SSMContDistribution.f(new("SSM"))))
+
+setClassUnion("SSMDistr", c("SSMDistribution.f",
+ "SSMConvDistribution.f"))
+
+###############################################################################
+#
+# SSM - Simulation classes
+#
+###############################################################################
+setClass("SSMsimulation", representation(
+ SSM = "SSM",
+ Distr = "SSMDistr",
+ RNGstate = "numeric",
+ states = "ArrayOrMatrix",
+ obs = "ArrayOrMatrix"),
+ prototype = prototype(SSM = new("SSM"),
+ Distr = new("SSMDistribution.f"),
+ RNGstate = structure(1, kind = as.list(RNGkind())),
+ states = matrix(1), obs = matrix(1)))
+
+setClass("SSMcontSimulation", representation(
+ states.id = "ArrayOrMatrix",
+ obs.id = "ArrayOrMatrix",
+ Ind.IO = "MatrixOrLogical",
+ Ind.AO = "MatrixOrLogical"
+ ),
+ prototype = prototype(
+ Distr = new("SSMConvDistribution.f"),
+ states.id = matrix(1), obs.id = matrix(1),
+ Ind.IO = FALSE, Ind.AO = FALSE
+ ),
+ validity = function(object){
+ fct <- function(m){
+ mt <- paste(deparse(substitute(m)),sep="",collapse="")
+ if(is.matrix(m)){
+ if(!all(is.logical(m)))
+ stop(gettextf("Matrix %s has to have logical entries.", mt))
+ return(TRUE)
+ }
+ else return(TRUE)
+ }
+ return(fct(object at Ind.IO)&&fct(object at Ind.AO))
+ },
+ contains = "SSMsimulation")
+
+
+
\ No newline at end of file
Added: pkg/robKalman/R/ACMfilt.R
===================================================================
--- pkg/robKalman/R/ACMfilt.R (rev 0)
+++ pkg/robKalman/R/ACMfilt.R 2009-03-18 15:41:02 UTC (rev 25)
@@ -0,0 +1,85 @@
+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 = matrix(x,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)
+
+ ### from version ... the return value $X[r]f is of
+ # dimension p x runs x (N+1)
+ # => have to cast it back to dimension p x (N+1)
+
+ X.ck <- matrix(ACMres$Xf[,1,], p, N+1); X.ck <- X.ck[,2:(N+1)]
+ X <- matrix(ACMres$Xrf[,1,], p, N+1); 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: pkg/robKalman/R/ACMfilter.R
===================================================================
--- pkg/robKalman/R/ACMfilter.R (rev 0)
+++ pkg/robKalman/R/ACMfilter.R 2009-03-18 15:41:02 UTC (rev 25)
@@ -0,0 +1,97 @@
+.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, i, ...)
+ {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 = FALSE))
+}
+
+.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
+## apsi, bpsi, cpsi ... tuning constants for Hampel's psi-function
+## (default: apsi=bpsi=2.5, cpsi=5.0)
+## flag ... character, if "weights" (default), use psi(t)/t to calculate
+## the weights; if "deriv", use psi'(t)
+
+ # to be compatible with parallel computing of a bunch of time series
+ y <- y[,1]
+ x1 <- x1[,1]
+
+ st <- rob1
+
+ K <- .getKG(S1, Z, V)
+
+ rst <- (y - x1[1])/st
+
+ ps <- psi(rst, apsi, bpsi, cpsi)
+ 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))))
+
+ return(list(x0 = x0, K = K, S0 = S0, Ind = ind, rob0 = rob1))
+}
Added: pkg/robKalman/R/AllGeneric.R
===================================================================
--- pkg/robKalman/R/AllGeneric.R (rev 0)
+++ pkg/robKalman/R/AllGeneric.R 2009-03-18 15:41:02 UTC (rev 25)
@@ -0,0 +1,10 @@
+#if(!isGeneric("type")){
+# setGeneric("type", function(object) standardGeneric("type"))
+#}
+#if(!isGeneric("center")){
+# setGeneric("center", function(object) standardGeneric("center"))
+#}
+#if(!isGeneric("center<-")){
+# setGeneric("center<-", function(object, value) standardGeneric("center<-"))
+#}
+#
Added: pkg/robKalman/R/AllGenerics.R
===================================================================
--- pkg/robKalman/R/AllGenerics.R (rev 0)
+++ pkg/robKalman/R/AllGenerics.R 2009-03-18 15:41:02 UTC (rev 25)
@@ -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: pkg/robKalman/R/AllInitialize.R
===================================================================
--- pkg/robKalman/R/AllInitialize.R (rev 0)
+++ pkg/robKalman/R/AllInitialize.R 2009-03-18 15:41:02 UTC (rev 25)
@@ -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: pkg/robKalman/R/AllPlot.R
===================================================================
--- pkg/robKalman/R/AllPlot.R (rev 0)
+++ pkg/robKalman/R/AllPlot.R 2009-03-18 15:41:02 UTC (rev 25)
@@ -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: pkg/robKalman/R/AllShow.R
===================================================================
--- pkg/robKalman/R/AllShow.R (rev 0)
+++ pkg/robKalman/R/AllShow.R 2009-03-18 15:41:02 UTC (rev 25)
@@ -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"))
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/robkalman -r 25
More information about the Robkalman-commits
mailing list