[Robkalman-commits] r80 - in branches/robkalman_pr/pkg/robKalman: R R/R-files from trunc R/vonBernhard demo
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu May 15 16:08:21 CEST 2025
Author: ruckdeschel
Date: 2025-05-15 16:08:20 +0200 (Thu, 15 May 2025)
New Revision: 80
Added:
branches/robkalman_pr/pkg/robKalman/R/R-files from trunc/0AllClass.R
branches/robkalman_pr/pkg/robKalman/R/R-files from trunc/AllGenerics.R
branches/robkalman_pr/pkg/robKalman/R/R-files from trunc/AllInitialize.R
branches/robkalman_pr/pkg/robKalman/R/R-files from trunc/AllPlot.R
branches/robkalman_pr/pkg/robKalman/R/R-files from trunc/AllShow.R
branches/robkalman_pr/pkg/robKalman/R/vonBernhard/
branches/robkalman_pr/pkg/robKalman/R/vonBernhard/classEKF-1.R
branches/robkalman_pr/pkg/robKalman/R/vonBernhard/classEKF-2.R
branches/robkalman_pr/pkg/robKalman/R/vonBernhard/classEKF.R
branches/robkalman_pr/pkg/robKalman/R/vonBernhard/classUKF-1.R
branches/robkalman_pr/pkg/robKalman/R/vonBernhard/classUKF.R
branches/robkalman_pr/pkg/robKalman/R/vonBernhard/generate-1.R
branches/robkalman_pr/pkg/robKalman/R/vonBernhard/generate-2.R
branches/robkalman_pr/pkg/robKalman/R/vonBernhard/generate-3.R
branches/robkalman_pr/pkg/robKalman/R/vonBernhard/generate.R
branches/robkalman_pr/pkg/robKalman/R/vonBernhard/recEFilter.R
Removed:
branches/robkalman_pr/pkg/robKalman/R/AllGeneric.R
Modified:
branches/robkalman_pr/pkg/robKalman/R/recFilter.R
branches/robkalman_pr/pkg/robKalman/demo/rLSdemo.R
Log:
parts Bernhard rLSdemo
Deleted: branches/robkalman_pr/pkg/robKalman/R/AllGeneric.R
===================================================================
--- branches/robkalman_pr/pkg/robKalman/R/AllGeneric.R 2025-05-15 14:06:59 UTC (rev 79)
+++ branches/robkalman_pr/pkg/robKalman/R/AllGeneric.R 2025-05-15 14:08:20 UTC (rev 80)
@@ -1,10 +0,0 @@
-#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: branches/robkalman_pr/pkg/robKalman/R/R-files from trunc/0AllClass.R
===================================================================
--- branches/robkalman_pr/pkg/robKalman/R/R-files from trunc/0AllClass.R (rev 0)
+++ branches/robkalman_pr/pkg/robKalman/R/R-files from trunc/0AllClass.R 2025-05-15 14:08:20 UTC (rev 80)
@@ -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_pr/pkg/robKalman/R/R-files from trunc/AllGenerics.R
===================================================================
--- branches/robkalman_pr/pkg/robKalman/R/R-files from trunc/AllGenerics.R (rev 0)
+++ branches/robkalman_pr/pkg/robKalman/R/R-files from trunc/AllGenerics.R 2025-05-15 14:08:20 UTC (rev 80)
@@ -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_pr/pkg/robKalman/R/R-files from trunc/AllInitialize.R
===================================================================
--- branches/robkalman_pr/pkg/robKalman/R/R-files from trunc/AllInitialize.R (rev 0)
+++ branches/robkalman_pr/pkg/robKalman/R/R-files from trunc/AllInitialize.R 2025-05-15 14:08:20 UTC (rev 80)
@@ -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_pr/pkg/robKalman/R/R-files from trunc/AllPlot.R
===================================================================
--- branches/robkalman_pr/pkg/robKalman/R/R-files from trunc/AllPlot.R (rev 0)
+++ branches/robkalman_pr/pkg/robKalman/R/R-files from trunc/AllPlot.R 2025-05-15 14:08:20 UTC (rev 80)
@@ -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_pr/pkg/robKalman/R/R-files from trunc/AllShow.R
===================================================================
--- branches/robkalman_pr/pkg/robKalman/R/R-files from trunc/AllShow.R (rev 0)
+++ branches/robkalman_pr/pkg/robKalman/R/R-files from trunc/AllShow.R 2025-05-15 14:08:20 UTC (rev 80)
@@ -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)
+# })
+#
+#
+#
+#
+#
Modified: branches/robkalman_pr/pkg/robKalman/R/recFilter.R
===================================================================
--- branches/robkalman_pr/pkg/robKalman/R/recFilter.R 2025-05-15 14:06:59 UTC (rev 79)
+++ branches/robkalman_pr/pkg/robKalman/R/recFilter.R 2025-05-15 14:08:20 UTC (rev 80)
@@ -1,4 +1,4 @@
-recursiveFilter <- function(Y, a, S, F, Q, Z, V,
+recursiveFilter.old <- function(Y, a, S, F, Q, Z, V,
initSc = .cKinitstep, predSc = .cKpredstep,
corrSc = .cKcorrstep,
initSr = NULL, predSr = NULL, corrSr = NULL,
Added: branches/robkalman_pr/pkg/robKalman/R/vonBernhard/classEKF-1.R
===================================================================
--- branches/robkalman_pr/pkg/robKalman/R/vonBernhard/classEKF-1.R (rev 0)
+++ branches/robkalman_pr/pkg/robKalman/R/vonBernhard/classEKF-1.R 2025-05-15 14:08:20 UTC (rev 80)
@@ -0,0 +1,105 @@
+#######################################################
+##
+## classical extended Kalman filter routines
+## author: Bernhard Spangl
+## version: 0.1 (changed: 2011-06-10, 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, ...)
+{
+## 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
+
+ return(list(x0=a, S0=S))
+
+}
+
+.cEKFpredstep <- function (x0, S0, F, Q, i,
+ v, u, controlF, # arguments of F
+ exQ, controlQ, ...) # arguments of Q
+{
+## 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
+
+ F <- F(t=i, x0=x0, v=v, u=u, control=controlF)
+ x1 <- F$x1
+ A <- F$A
+ B <- F$B
+
+ Qreturn <- Q(t=i, x0=x0, exQ=exQ, control=controlQ)
+ Q <- Qreturn$Q
+
+ return(list(x1=x1, S1=.getpredCov(S0=S0, A=A, B=B, Q=Q)))
+
+}
+
+.cEKFcorrstep <- function (y, x1, S1, Z, V, i,
+ eps, w, controlZ, # arguments of Z
+ exV, controlV, ...) # arguments of V
+{
+## 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
+
+ Z <- Z(t=i, x1=x1, eps=eps, w=w, control=controlZ)
+ yhat <- Z$y
+ C <- Z$C
+ D <- Z$D
+
+ Vretrun <- V(t=i, x1=x1, exV=exV, control=controlV)
+ V <- Vreturn$V
+
+ Delta <- .getDelta(S1=S1, C=C, D=D, V=V)
+ K <- .getKG(S1=S1, Z=C, Delta=Delta)
+ DeltaY <- y - yhat
+ x0 <- x1 + K %*% DeltaY
+ S0 <- .getcorrCov(S1=S1, K=K, Z=C)
+
+ return(list(x0=x0, K=K, S0=S0, Delta=Delta, DeltaY=DeltaY))
+
+}
+
Added: branches/robkalman_pr/pkg/robKalman/R/vonBernhard/classEKF-2.R
===================================================================
--- branches/robkalman_pr/pkg/robKalman/R/vonBernhard/classEKF-2.R (rev 0)
+++ branches/robkalman_pr/pkg/robKalman/R/vonBernhard/classEKF-2.R 2025-05-15 14:08:20 UTC (rev 80)
@@ -0,0 +1,105 @@
+#######################################################
+##
+## classical extended Kalman filter routines
+## author: Bernhard Spangl
+## version: 0.1 (changed: 2011-06-10, 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, ...)
+{
+## 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
+
+ return(list(x0=a, S0=S))
+
+}
+
+.cEKFpredstep <- function (x0, S0, F, Q, i,
+ v, u, controlF, # arguments of F
+ exQ, controlQ, ...) # arguments of Q
+{
+## 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
+
+ F <- F(t=i, x0=x0, v=v, u=u, control=controlF)
+ x1 <- F$x1
+ A <- F$A
+ B <- F$B
+
+ Qreturn <- Q(t=i, x0=x0, exQ=exQ, control=controlQ)
+ Q <- Qreturn$Q
+
+ return(list(x1=x1, S1=.getpredCov(S0=S0, A=A, B=B, Q=Q)))
+
+}
+
+.cEKFcorrstep <- function (y, x1, S1, Z, V, i,
+ eps, w, controlZ, # arguments of Z
+ exV, controlV, ...) # arguments of V
+{
+## 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
+
+ Z <- Z(t=i, x1=x1, eps=eps, w=w, control=controlZ)
+ yhat <- Z$y
+ C <- Z$C
+ D <- Z$D
+
+ Vretrun <- V(t=i, x1=x1, exV=exV, control=controlV)
+ V <- Vreturn$V
+
+ Delta <- .getDelta(S1=S1, C=C, D=D, V=V)
+ K <- .getKG(S1=S1, Z=C, Delta=Delta)
+ DeltaY <- y - yhat
+ x0 <- x1 + K %*% DeltaY
+ S0 <- .getcorrCov(S1=S1, K=K, Z=C)
+
+ return(list(x0=x0, K=K, S0=S0, Delta=Delta, DeltaY=DeltaY))
+
+}
+
Added: branches/robkalman_pr/pkg/robKalman/R/vonBernhard/classEKF.R
===================================================================
--- branches/robkalman_pr/pkg/robKalman/R/vonBernhard/classEKF.R (rev 0)
+++ branches/robkalman_pr/pkg/robKalman/R/vonBernhard/classEKF.R 2025-05-15 14:08:20 UTC (rev 80)
@@ -0,0 +1,113 @@
+#######################################################
+##
+## classical extended Kalman filter routines
+## author: Bernhard Spangl
+## version: 0.2 (changed: 2011-06-21, 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, 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
+
+ F <- F(t=i, x0=x0, v=v, u=u, control=controlF)
+ x1 <- F$x1
+ A <- F$A
+ B <- F$B
+
+ Qreturn <- Q(t=i, x0=x0, exQ=exQ, control=controlQ)
+ Q <- Qreturn$Q
+
+ return(list(x1=x1, S1=.getpredCov(S0=S0, A=A, B=B, Q=Q),
+ controlPred=controlPred))
+
+}
+
+.cEKFcorrstep <- 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
+
+ Z <- Z(t=i, x1=x1, eps=eps, w=w, control=controlZ)
+ yhat <- Z$y
+ C <- Z$C
+ D <- Z$D
+
+ Vretrun <- V(t=i, x1=x1, exV=exV, control=controlV)
+ V <- Vreturn$V
+
+ Delta <- .getDelta(S1=S1, C=C, D=D, V=V)
+ K <- .getKG(S1=S1, Z=C, Delta=Delta)
+ DeltaY <- y - yhat
+ x0 <- x1 + K %*% DeltaY
+ S0 <- .getcorrCov(S1=S1, K=K, Z=C)
+
+ return(list(x0=x0, K=K, S0=S0, Delta=Delta, DeltaY=DeltaY,
+ controlCorr=controlCorr))
+
+}
+
+
Added: branches/robkalman_pr/pkg/robKalman/R/vonBernhard/classUKF-1.R
===================================================================
--- branches/robkalman_pr/pkg/robKalman/R/vonBernhard/classUKF-1.R (rev 0)
+++ branches/robkalman_pr/pkg/robKalman/R/vonBernhard/classUKF-1.R 2025-05-15 14:08:20 UTC (rev 80)
@@ -0,0 +1,148 @@
+#######################################################
+##
+## 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, ...)
+{
+## 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
+## control ... list containing alpha, beta and kappa
+
+ alpha <- control$alpha
+ beta <- control$beta
+ kappa <- control$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)
+
+ control$gamma <-gamma
+ control$Wm <- Wm
+ control$Wc <- Wc
+
+ return(list(x0=a, S0=S, control=control))
+
+}
+
+.cUKFpredstep <- function (x0, S0, F, Q, i,
+ v, u, controlF, # arguments of F
+ exQ, controlQ, # arguments of Q
+ control, ...) # 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
+## control ... control parameters of used filter
+
+ gamma <- control$gamma
+ Wm <- control$Wm
+ Wc <- control$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)
+
+ control$X1x <- X1x
+
+ return(list(x1=x1, S1=S1, control=conrol))
+
+}
+
+.cUKFcorrstep <- function (y, x1, S1, Z, V, i,
+ eps, w, controlZ, # arguments of Z
+ exV, controlV, # arguments of V
+ control, ...) # 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
+## control ... control parameters of used filter
+
+ gamma <- control$gamma
+ Wm <- control$Wm
+ Wc <- control$Wc
+ X1x <- control$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))
+
+}
+
+
Added: branches/robkalman_pr/pkg/robKalman/R/vonBernhard/classUKF.R
===================================================================
--- branches/robkalman_pr/pkg/robKalman/R/vonBernhard/classUKF.R (rev 0)
+++ branches/robkalman_pr/pkg/robKalman/R/vonBernhard/classUKF.R 2025-05-15 14:08:20 UTC (rev 80)
@@ -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)
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/robkalman -r 80
More information about the Robkalman-commits
mailing list