[Robkalman-commits] r41 - in branches: robKalman_itwm/pkg/robKalman/R robKalman_itwm/pkg/robKalman/man robkalman_pr/pkg/robKalman/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Sep 29 16:07:42 CEST 2009
Author: ruckdeschel
Date: 2009-09-29 16:07:42 +0200 (Tue, 29 Sep 2009)
New Revision: 41
Added:
branches/robKalman_itwm/pkg/robKalman/R/EM_Algo.R
branches/robKalman_itwm/pkg/robKalman/R/Estep.R
branches/robKalman_itwm/pkg/robKalman/R/KalmanMethods.R
branches/robKalman_itwm/pkg/robKalman/R/MakePositive.R
branches/robKalman_itwm/pkg/robKalman/R/Mstep.R
branches/robKalman_itwm/pkg/robKalman/R/SSM.R
branches/robKalman_itwm/pkg/robKalman/R/getSmooth.R
branches/robKalman_itwm/pkg/robKalman/R/getSmoothFixInterval.R
branches/robKalman_itwm/pkg/robKalman/R/getSmoothFixPoint.R
branches/robKalman_itwm/pkg/robKalman/R/getSmootherR.txt
branches/robKalman_itwm/pkg/robKalman/R/initEMstepclass.R
branches/robKalman_itwm/pkg/robKalman/R/recSmoother.R
branches/robKalman_itwm/pkg/robKalman/man/EM_Algo.Rd
branches/robKalman_itwm/pkg/robKalman/man/MStepShST.Rd
Modified:
branches/robKalman_itwm/pkg/robKalman/R/0AllClass.R
branches/robKalman_itwm/pkg/robKalman/R/ACMfilt.R
branches/robKalman_itwm/pkg/robKalman/R/ACMfilter.R
branches/robKalman_itwm/pkg/robKalman/R/Util.R
branches/robKalman_itwm/pkg/robKalman/R/calibrateRLS.R
branches/robKalman_itwm/pkg/robKalman/R/classKalman.R
branches/robKalman_itwm/pkg/robKalman/R/mACMfilter.R
branches/robKalman_itwm/pkg/robKalman/R/rLSfilter.R
branches/robKalman_itwm/pkg/robKalman/R/recFilter.R
branches/robKalman_itwm/pkg/robKalman/R/simulateSScont.R
branches/robKalman_itwm/pkg/robKalman/man/0robKalman-package.Rd
branches/robKalman_itwm/pkg/robKalman/man/calibrateRLS.Rd
branches/robKalman_itwm/pkg/robKalman/man/internalACM.Rd
branches/robKalman_itwm/pkg/robKalman/man/internalKalman.Rd
branches/robKalman_itwm/pkg/robKalman/man/internalrLS.Rd
branches/robKalman_itwm/pkg/robKalman/man/recFilter.Rd
branches/robKalman_itwm/pkg/robKalman/man/simulateSScont.Rd
branches/robkalman_pr/pkg/robKalman/R/getSmoothFixPoint.R
Log:
prepared branches/robKalman_itwm for work by Cezar
Modified: branches/robKalman_itwm/pkg/robKalman/R/0AllClass.R
===================================================================
--- branches/robKalman_itwm/pkg/robKalman/R/0AllClass.R 2009-09-03 09:35:13 UTC (rev 40)
+++ branches/robKalman_itwm/pkg/robKalman/R/0AllClass.R 2009-09-29 14:07:42 UTC (rev 41)
@@ -1,3 +1,5 @@
+################ general package preparation code, nothing to do with classes:
+
.onLoad <- function(lib, pkg){
require("methods", character = TRUE, quietly = TRUE)
@@ -12,3 +14,381 @@
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)
+ })
+
+
+
+#
+
+### 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 = "timeSeries"), ## 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 = timeSeries(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",
+ Delta = "array",
+ DeltaY = "array",
+ time = "timeSeries"),
+ 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)),
+ Delta = "array",
+ DeltaY = "array",
+ time = timeSeries(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",
+ Deltar = "array",
+ DeltaYr = "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")
+
+###############################################################################
+#
+# Smoother classes
+#
+###############################################################################
+
+setClass("recSmoother", representation(
+ X.smoothed = "array",
+ Cov.smoothed = "array",
+ J = "array"),
+ prototype = prototype(name="classical Kalman Smoother",
+ 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 = timeSeries(1)),
+
+ contains="recFilter")
+
+
+setClass("robrecSmoother", representation(
+ X.rob.smoothed = "array",
+ Cov.rob.smoothed = "array",
+ J.rob = "array"),
+ prototype = prototype(
+ name="rLS Smoother",
+ X.rob.smoothed = array(1,dim=c(1,1,1,1)),
+ Cov.rob.smoothed = array(1,dim = c(1,1,1,1)),
+ J.rob = array(1,dim = c(1,1,1,1))),
+ contains = "recSmoother")
+
+###############################################################################
+#
+# Control classes
+#
+###############################################################################
+#setClass("ACMcontrol",representation())
+#setClass("rLScontrol",representation())
+
+
+setClass("RecFiltControl", representation(
+ name = "character",
+ init = "function",
+ predict = "function",
+ correct = "function"),
+ contains = "VIRTUAL")
+
+setClass("RecSmoothControl", representation(
+ smooth = "function"),
+ contains = c("VIRTUAL","RecFiltControl"))
+
+setClass("KalmanFiltControl",
+ prototype = prototype(
+ name = "classical Kalman Filter",
+ init = .cKinitstep,
+ predict = .cKpredstep,
+ correct = .cKcorrstep),
+ contains="RecFiltControl"
+ )
+setClass("KalmanSmoothControl",
+ prototype = prototype(
+ name = "classical Kalman Smoother",
+ smooth = .cKsmoothstep),
+ contains="KalmanFiltControl"
+ )
+
+setClass("RobRecFilterControl", 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"
+ )
+setClass("RobRecSmoothControl", representation(
+ smooth.rob = "function"),
+ prototype = prototype(
+ name.rob = "rLS Filter",
+ init.rob = .cKinitstep,
+ predict.rob = .cKpredstep,
+ correct.rob = .rLScorrstep,
+ smooth.rob = .cKsmoothstep,
+ controls = list(b = 2, norm = EuclideanNorm)
+ ),
+ contains="RecFiltControl"
+ )
+# = paste(gettext(
+# "Control set and init, prediction and"
+# ),gettext(
+# "correction step for the classical Kalman Filter\n"
+# ),gettext(
+# "and the rLS Filter"
+# )),
+
+###############################################################################
+#
+# multivariate Distribution classes
+#
+###############################################################################
+
+setClass("SSMDistribution.f", representation(
+ r.init = "function",
+ r.innov = "function",
+ r.obs = "function"),
+ prototype = prototype(r.init = mvrnorm,
+ r.innov = mvrnorm,
+ r.obs = mvrnorm))
+
+setClass("SSMellDistribution.f", representation(
+ m.init = "sHyperparamtype",
+ S.init = "Hyperparamtype",
+ m.innov = "sHyperparamtype",
+ S.innov = "Hyperparamtype",
+ m.obs = "sHyperparamtype",
+ S.obs = "Hyperparamtype"),
+ prototype = prototype(m.init = 0, S.init = matrix(1),
+ m.innov = 0, S.innov = matrix(1),
+ m.obs = 0, S.obs = matrix(1)),
+ contains = "SSMDistribution.f")
+
+setClass("SSMConvDistribution.f", representation(
+ ideal = "SSMDistribution.f",
+ cont = "SSMellDistribution.f",
+ r.IO = "numeric",
+ r.AO = "numeric"),
+ prototype = prototype(ideal = new("SSMellDistribution.f"),
+ cont = new("SSMellDistribution.f"),
+ r.IO = 0, r.AO = 0),
+ validity = function(object){
+ if (object at r.IO<0||object at r.AO<0||object at r.IO>1||object at r.AO>1)
+ stop("Radii must be between 0 and 1")
+ return(TRUE)})
+
+
+###############################################################################
+#
+# SSM + Distribution classes
+#
+###############################################################################
+setClass("SSMwithDistribution", representation(
+ SSM = "SSM",
+ Distribution = "SSMDistribution.f"),
+ prototype = prototype(SSM = new("SSM"),
+ Distribution = SSMDistribution.f(new("SSM"))))
+setClass("SSMwithConvDistribution", representation(
+ SSM = "SSM",
+ Distribution = "SSMConvDistribution.f"),
+ prototype = prototype(SSM = new("SSM"),
+ Distribution = SSMContDistribution.f(new("SSM"))))
+
+setClassUnion("SSMDistr", c("SSMDistribution.f",
+ "SSMConvDistribution.f"))
+
+###############################################################################
+#
+# SSM - Simulation classes
+#
+###############################################################################
+setClass("SSMsimulation", representation(
+ SSM = "SSM",
+ Distr = "SSMDistr",
+ RNGstate = "numeric",
+ states = "ArrayOrMatrix",
+ obs = "ArrayOrMatrix"),
+ prototype = prototype(SSM = new("SSM"),
+ Distr = new("SSMDistribution.f"),
+ RNGstate = structure(1, kind = as.list(RNGkind())),
+ states = matrix(1), obs = matrix(1)))
+
+setClass("SSMcontSimulation", representation(
+ states.id = "ArrayOrMatrix",
+ obs.id = "ArrayOrMatrix",
+ Ind.IO = "MatrixOrLogical",
+ Ind.AO = "MatrixOrLogical"
+ ),
+ prototype = prototype(
+ Distr = new("SSMConvDistribution.f"),
+ states.id = matrix(1), obs.id = matrix(1),
+ Ind.IO = FALSE, Ind.AO = FALSE
+ ),
+ validity = function(object){
+ fct <- function(m){
+ mt <- paste(deparse(substitute(m)),sep="",collapse="")
+ if(is.matrix(m)){
+ if(!all(is.logical(m)))
+ stop(gettextf("Matrix %s has to have logical entries.", mt))
+ return(TRUE)
+ }
+ else return(TRUE)
+ }
+ return(fct(object at Ind.IO)&&fct(object at Ind.AO))
+ },
+ contains = "SSMsimulation")
+
+
+
\ No newline at end of file
Modified: branches/robKalman_itwm/pkg/robKalman/R/ACMfilt.R
===================================================================
--- branches/robKalman_itwm/pkg/robKalman/R/ACMfilt.R 2009-09-03 09:35:13 UTC (rev 40)
+++ branches/robKalman_itwm/pkg/robKalman/R/ACMfilt.R 2009-09-29 14:07:42 UTC (rev 41)
@@ -52,10 +52,10 @@
## 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)
+ 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)
### from version ... the return value $X[r]f is of
# dimension p x runs x (N+1)
Modified: branches/robKalman_itwm/pkg/robKalman/R/ACMfilter.R
===================================================================
--- branches/robKalman_itwm/pkg/robKalman/R/ACMfilter.R 2009-09-03 09:35:13 UTC (rev 40)
+++ branches/robKalman_itwm/pkg/robKalman/R/ACMfilter.R 2009-09-29 14:07:42 UTC (rev 41)
@@ -19,16 +19,16 @@
}
##steps for classical Kalman filter (cK)
-.ACMinitstep <- function(a, S, ...)
+.ACMinitstep <- function(a, S, i, ...)
{dots <- list(...)
if(hasArg("s0"))
- s0<-dots$"s0"
+ s0 <- dots$"s0"
else
- s0<-NULL
+ s0 <- NULL
list( x0 = a, S0 = S, s0 = s0)}
-.ACMpredstep <- function (x0, S0, F, Q, rob0, s0, ...) ### S=P F= Phi
+.ACMpredstep <- function (x0, S0, F, Q, i, rob0, s0, ...) ### S=P F= Phi
{
###########################################
##
@@ -45,10 +45,11 @@
## 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))
+ return(list(x1 = F %*% x0, S1 = S1, rob1 = sqrt(S1[1, 1] + s0), Ind = FALSE))
}
-.ACMcorrstep <- function (y, x1, S1, Z, V, rob1, dum=NULL, psi, apsi, bpsi, cpsi, flag, ...)
+.ACMcorrstep <- function (y, x1, S1, Z, V, i, rob1, dum = NULL,
+ psi, apsi, bpsi, cpsi, flag, ...)
{
###########################################
##
@@ -67,18 +68,22 @@
## 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)
+## 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)
+ ps <- psi(rst, apsi, bpsi, cpsi)[1,1]
dx <- K * st * ps
x0 <- x1 + dx
@@ -88,5 +93,7 @@
S0 <- .getcorrCovACM(S1, K, Z, W = w*diag(rep(1, nrow(Z))))
- return(list(x0 = x0, K = K, S0 = S0, Ind=ind, rob0=rob1))
+ 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_itwm/pkg/robKalman/R/EM_Algo.R
===================================================================
--- branches/robKalman_itwm/pkg/robKalman/R/EM_Algo.R (rev 0)
+++ branches/robKalman_itwm/pkg/robKalman/R/EM_Algo.R 2009-09-29 14:07:42 UTC (rev 41)
@@ -0,0 +1,99 @@
+EMAlgoSSM <- function(y, a, Z, initStepEM = initStepSSMclass,
+ EStepEM = EStepSSMclass,
+ MStepEM = MStepSSMclass,
+ eabs, erel, maxit,
+ verbose = FALSE, ...)
+{
+
+##main/frame-function for the EM Algorithm(à la Shumway&Stoffer,1981)
+
+ eabs <- rep(eabs, length.out = 3)
+ erel <- rep(erel, length.out = 3)
+ if(is.null(names(eabs))) names(eabs) <- c("F","Q","V")
+ if(is.null(names(erel))) names(erel) <- c("F","Q","V")
+
+ qd <- if(length(dim(Z)) 1 else (dim(Y))[1]
+ pd <- if(length(dim(Z)) 1 else (dim(Z))[2]
+ tt <- (dim(Y))[3]
+ runs <- (dim(Y))[2]
+
+
+ Fa <- array(0,dim=c(pd,pd,runs))
+ Qa <- array(0,dim=c(pd,pd,runs))
+ Va <- array(0,dim=c(qd,qd,runs))
+ ka <- numeric(runs)
+
+ for(i in 1:runs){
+ resinitStepEM <- initStepEM(y,Z)
+
+ F <- resinitStepEM$F
+ Q <- resinitStepEM$Q
+ V <- resinitStepEM$V
+
+
+ log.F <- TRUE
+ log.Q <- TRUE
+ log.V <- TRUE
+ k <- 0
+
+ while((k<maxit) || log.F || log.V || log.Q )
+ {
+
+
+##NOTE: In order to be able to run the EM Algorithm recursively, we need to compute
+##those values first, so we just need a starting value for the differences of the
+##norms for the parameters V,F and Q. Hence we just run the E and M-steps once.
+##For the E-Step, one also needs input parameters a and S, so:
+ ##a<-miu_init;
+ ##S<-sigma_init;
+
+ resEStepEM <- EStepEM(y, a, S=Q, F, Q, Z, V, dropRuns = FALSE, ...)
+
+ F_old <- F
+ Q_old <- Q
+ V_old <- V
+
+## define now the parameters needed for the M-Step:
+ x <- resEStepEM$x
+ P <- resEStepEM$P
+ Pdep <- resEStepEM$Pdep
+
+##Call function MStepEM()...
+
+ resMStepEM <- MStepEM(M = Z, y = y, # matrix(y[,1,],qd,T),
+ x = x, P = P, Pdep = Pdep, ...)
+
+
+ F <- resMStepEM$F
+ Q <- resMStepEM$Q
+ V <- resMStepEM$V
+
+ log.F <- .isDeltaSmall(F_old,F,e.rel["F"],e.abs["F"])
+ log.Q <- .isDeltaSmall(Q_old,Q,e.rel["Q"],e.abs["Q"])
+ log.V <- .isDeltaSmall(V_old,V,e.rel["V"],e.abs["V"])
+
+
+ k <- k+1
+ if(verbose)
+ print(list(k=k,x=x,F=F,Q=Q,V=V));
+ }
+
+ Fa[,,i] <- F
+ Qa[,,i] <- Q
+ Va[,,i] <- V
+ ka[i] <- k
+ }
+#############END of Algorithm###################
+##OUTPUT:
+list(F=Fa, Q=Qa, V=Va, Z=Z, k=ka)
+
+}
+
+##Call this function(see how it works)
+##EmAlgo.intern(Z,y,Q_init=Q,F_init=F,V_init=V,EStepEM=recursiveFilter,MStepEM=MStep,eabs,erel,maxit,noise=NULL)
+
+.isDeltaSmall <- function(x.old, x.new, e.rel, e.abs){
+ no.a <- (sum((x.old-x.new)^2))^.5
+ no.r <- no.a / (sum((x.old)^2))^.5
+ return(no.a < e.abs && no.r < e.rel)
+ }
Added: branches/robKalman_itwm/pkg/robKalman/R/Estep.R
===================================================================
--- branches/robKalman_itwm/pkg/robKalman/R/Estep.R (rev 0)
+++ branches/robKalman_itwm/pkg/robKalman/R/Estep.R 2009-09-29 14:07:42 UTC (rev 41)
@@ -0,0 +1,50 @@
+EStepSSM <- function(y, a, S = Q, F, Z, Q, V, initSc = .cKinitstep, predSc = .cKpredstep,
+ corrSc = .cKcorrstep, smoothSc = .cKsmoothstep,
+ initSr = NULL, predSr = NULL, corrSr = NULL, smoothSr = NULL,
+ calibrateFilter = NULL,
+ nsim = 0, seed = NULL, ..., dropRuns = FALSE)
+{Y <- array(y, dim=c(nrow(y),1,ncol(y)))
+
+ dots <- list(...)
+ if(!is.null(calibrateFilter)){
+ SS <- limitS(S = S, F = F, Q = Q, Z = Z, V = V)
+ calibF <- calibrateFilter(..., S = SS, Z = Z, V = V)
+ }
+ resFilter <- do.call(recursiveFilter, args = c( list( Y=Y, a=a, S=Q, F=F, Q=Q,
+ Z=Z, V=V, initSc = initSc, predSc = predSc,
+ corrSc = corrSc,
+ initSr = initSr, predSr = predSr, corrSr = corrSr,
+ nsim = nsim, seed = seed), calibF, dots,
+ dropRuns = dropRuns))
+ res <- recursiveFPSmoother( Y = Y,
+ F = F, Q = Q, V = V, Z = Z,
+ smoothSc = smoothSc, smoothSr = smoothSr,
+ KG = resFilter$KG, xf.c = resFilter$xf,
+ S0.c = resFilter$S0, S1.c = resFilter$S1,
+ xf.r = resFilter$xrf,
+ S0.r = resFilter$Sr0, S1.r = resFilter$Sr1,
+ KG.r = resFilter$KGr, ..., dropRuns = dropRuns)
+
+ if(is.null(initSr)&&is.null(predSr)&&is.null(corrSr)&&is.null(smoothSr))
+ return(list(M = resFilter$KG, x = res$xS, P = res$SS, Pdep = res$SS1))
+ else
+ return(list(M = resFilter$KGr, x = res$xSr, P = res$SSr,
+ Pdep = res$SS1r))
+}
+
+EStepSSMclass <- function(y, a, S = Q, F, Z, Q, V, ...)
+{EStepSSM(y = y, a = a, S = S, F = F, Z = Z , Q = Q, V = V,
+ initSc = .cKinitstep, predSc = .cKpredstep,
+ corrSc = .cKcorrstep, smoothSc = .cKsmoothstep,
+ initSr = NULL, predSr = NULL, corrSr = NULL, smoothSr = NULL,
+ calibrateFilter = NULL,
+ nsim = 0, seed = NULL, ..., dropRuns = FALSE)}
+
+EStepSSMrLS.AO <- function(y, a, S = Q, F, Z, Q, V, ...)
+{EStepSSM(y = y, a = a, S = S, F = F, Z = Z , Q = Q, V = V,
+ initSc = .cKinitstep, predSc = .cKpredstep,
+ corrSc = .cKcorrstep, smoothSc = .cKsmoothstep,
+ initSr = .cKinitstep, predSr = .cKpredstep,
+ corrSr = .rLScorrstep, smoothSr = .cKsmoothstep,
+ calibrateFilter = rLScalibrateB,
+ nsim = 0, seed = NULL, ..., dropRuns = FALSE)}
Added: branches/robKalman_itwm/pkg/robKalman/R/KalmanMethods.R
===================================================================
--- branches/robKalman_itwm/pkg/robKalman/R/KalmanMethods.R (rev 0)
+++ branches/robKalman_itwm/pkg/robKalman/R/KalmanMethods.R 2009-09-29 14:07:42 UTC (rev 41)
@@ -0,0 +1,41 @@
+setMethod("kalmanRob", signature(method = "robrecControl", smooth = "missing"),
+ function(method = rLSControl(), Y, SSM, nsim = 0, seed = NULL){
+ SSMa <- makeArrayRepresentation(SSM)
+ erg <- do.call(recursiveFilter,
+ args = c(list(Y = Y, a = SSM at a, S = SSMa at S,
+ F = SSMa at F, Q = SSMa at Q, Z = SSMa at Z, V = SSMa at V,
+ nsim = nsim, seed = seed,
+ initSc = init(method),
+ predSc = predict(method),
+ corrSc = correct(method),
+ initSr = init.rob(method),
+ predSr = predict.rob(method),
+ corrSr = correct.rob(method)),
+ controls(method), dropRuns = FALSE)
+ )
+
+ return(generateRobRecFilter(
+ name = name(method), name.rob = name.rob(method),
+ SSM = SSM, Y = Y, time = SSM at time,
+ Xf = erg$Xf, Xp = erg$Xp, S0 = erg$S0, S1 = erg$S1, KG = erg$KG,
+ Xrf = erg$Xrf, Xrp = erg$Xrp, Sr0 = erg$Sr0, Sr1 = erg$Sr1,
+ KGr = erg$KGr, IndIO = erg$IndIO, IndAO = erg$IndAO,
+ rob0L = erg$rob0L, rob1L = erg$rob1L, St0s = erg$St0s,
+ St1s = erg$St1s, nsim = erg$nsim, RNGstate = erg$RNGstate
+ ))
+ })
+
+setMethod("kalman", signature(smooth = "missing"),
+ function(Y, SSM){
+ method <- KalmanControl()
+ SSMa <- makeArrayRepresentation(SSM)
+ erg <- recursiveFilter(Y = Y, a = SSM at a, S = SSMa at S,
+ F = SSMa at F, Q = SSMa at Q, Z = SSMa at Z, V = SSMa at V,
+ initSc = init(method),
+ predSc = predict(method),
+ corrSc = correct(method), dropRuns = FALSE)
+ return(generateRecFilter(
+ name = name(method), SSM = SSM, Y = Y, time = SSM at time,
+ Xf = erg$Xf, Xp = erg$Xp, S0 = erg$S0, S1 = erg$S1, KG = erg$KG
+ ))})
+
Added: branches/robKalman_itwm/pkg/robKalman/R/MakePositive.R
===================================================================
--- branches/robKalman_itwm/pkg/robKalman/R/MakePositive.R (rev 0)
+++ branches/robKalman_itwm/pkg/robKalman/R/MakePositive.R 2009-09-29 14:07:42 UTC (rev 41)
@@ -0,0 +1,10 @@
+MakePositive <- function(matrix){
+### forces negative eigen values
+### of a symmetric matrix to 0
+if(length(matrix)==1) return ((matrix>0)*matrix)
+X.eigen <- eigen(matrix, symmetric=TRUE)
+V <- X.eigen$vectors
+D <- X.eigen$values
+D[X.eigen$values<0] <- 0
+return(V%*%diag(D)%*%t(V))
+}
\ No newline at end of file
Added: branches/robKalman_itwm/pkg/robKalman/R/Mstep.R
===================================================================
--- branches/robKalman_itwm/pkg/robKalman/R/Mstep.R (rev 0)
+++ branches/robKalman_itwm/pkg/robKalman/R/Mstep.R 2009-09-29 14:07:42 UTC (rev 41)
@@ -0,0 +1,93 @@
+MStepSSMclass <- function(Y,M,x,P,Pdep)
+{
+
+######################################################################################
+##FIRST APPROACH:
+##########################################################################
+##This algorithm is based on the paper by Shumway&Stoffer,1981: #
+##"An approach to time series smoothing and forecasting using the #
+##EM algorithm". #
+##If one knows the values for the parameters miu, sigma, fi, Q and R, #
+##one can calculate the Kalman smoothing estimators as described in #
+##the *.Rd-documentation file. This is actually the goal of this routine.#
+##########################################################################
+
+##SOME MORE EXPLANATION REGARDING THE IMPLEMENTATION:
+###################################################################################
+##What we would like to do is to compute equations (12),(13), (14). #
+##For this purpose, we assume(just in this beginner step), that #
+##we can calculate A,B and C as in equations (9), (10), (11) respectively. #
+##Hence, let us assume, that we already performed the estimation step, and that #
+##this step is realized within the function "Estimation", whose output parameters:#
+##x,P, and Pdep are calculated via equations (A3)-(A12) from Appendix A. So: #
+###################################################################################
+##
+##number of observations:
+## to be drawn from input
+##
+runs <- dim(Y)[2]
+n <- dim(x)[2]-1
+p <- dim(P)[1]
+q <- if(is.null(dim(y))) 1 else dim(y)[1]
+
+##Like mentioned above, we use the values resulted by the "E"-step, i.e.
+##Estimation-step:
+
+## comment P.R.: only works for one run at a time --- hence arguments are now:
+
+## M observation matrix (Z in our notation); q x q x n
+## y observations; q x n
+## x smoothed value x_t|n computed with getSmoother(); p x (n+1)
+## P smoother cov's Sigma_{t|n} comp. with getSmoother(); p x p x (n+1)
+## Pdep smoother cov Cov(X_t,X_t-1|y_1:n) comp. with getSmoother(); p x p x n
+
+## output: list of estimated F,Q,V
+
+##for(t in 1:n )
+##{
+##x[,,t]<-Estimation$x[,,t];
+##P[,,t]<-Estimation$P[,,t];
+##Pdep[,,t]<-Estimation$Pdep[,,t];
+##}
+##This part just gives an outline about the way we will procede, when we will also have
+##estimation step(Step 1 in the paper) done. For the calculations presented now,
+##we asssume, we have some input x, P, Pdep.The file MStep.Rd explains the general setting of
+##the algorithm and its goal.
+########################################################################################
+
+##initialize matrix values:
+
+ R0<-array(0,dim=c(q,q,n))
+
+ null <- (1:n)
+ plus <- null +1
+
+ x0 <- x[,null,drop=FALSE]
+ x1 <- x[,plus,drop=FALSE]
+ P0 <- P[,,null,drop=FALSE]
+ P1 <- P[,,plus,drop=FALSE]
+
+##Now: calculate the values of A, B and C as in eq. (9), (10), (11) respectively.
+ A <- apply(P0,c(1,2),sum) + x0 %*% t(x0)
+
+ C <- apply(P1,c(1,2),sum) + x1 %*% t(x1)
+
+ B <- apply(Pdep,c(1,2),sum) + x1 %*% t(x0)
+
+##having computed these values, one can compute now:
+ F <- B %*% ginv(A)
+ Q <- (C-F%*%t(B))/n
+
+ yt <- 0 *y
+ for (i in 1:n){
+ yt[,i] <- M[,,i]%*%x[,i+1]
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/robkalman -r 41
More information about the Robkalman-commits
mailing list