[Robkalman-commits] r92 - pkg/robKalman/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat May 17 20:40:00 CEST 2025
Author: bspangl
Date: 2025-05-17 20:39:59 +0200 (Sat, 17 May 2025)
New Revision: 92
Modified:
pkg/robKalman/R/Psi.R
pkg/robKalman/R/Util.R
pkg/robKalman/R/mACMfilter.R
Log:
aktuellsten Stand nach r-forge transformiert
Modified: pkg/robKalman/R/Psi.R
===================================================================
--- pkg/robKalman/R/Psi.R 2025-05-15 20:35:54 UTC (rev 91)
+++ pkg/robKalman/R/Psi.R 2025-05-17 18:39:59 UTC (rev 92)
@@ -149,21 +149,52 @@
return(x*dummy)
}
-jacobian.Hampel <- function (x, ...)
+
+## 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, ...)
+## }
+
+
+jacobian.Hampel <- function (x, a=2, b=4, c=8)
{
###########################################
##
## 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)
+## author: Bernhard Spangl
+## version: 0.1 (2025-05-17)
##
###########################################
## Paramters:
## x ... vector
-
- jacobian(mvpsiHampel, x, ...)
+## a, b, c ... tuning constants
+ x.norm2 <- drop(x%*%x)
+ x.norm <- sqrt(x.norm2)
+ I <- diag(rep(1, length(x)))
+ x.outer <- outer(x, x)/x.norm2
+ if (x.norm < a) {
+ I
+ } else if (x.norm < b) {
+ (I - x.outer)*a/x.norm
+ } else if (x.norm < c) {
+ ((I - x.outer)*c/x.norm - I)*a/(c-b)
+ } else {
+ diag(rep(0, length(x)))
+ }
}
Modified: pkg/robKalman/R/Util.R
===================================================================
--- pkg/robKalman/R/Util.R 2025-05-15 20:35:54 UTC (rev 91)
+++ pkg/robKalman/R/Util.R 2025-05-17 18:39:59 UTC (rev 92)
@@ -35,7 +35,7 @@
## 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)
+## version: 0.3 (2025-02-06)
##
###########################################
@@ -43,7 +43,8 @@
## X ... symmetric and positive semidefinite matrix
if ((ncol(X) == 1) && (nrow(X) == 1))
- return(list(X.det=X,
+ return(list(X.det=X,
+ X.inv=matrix(1/X),
X.sqrt=matrix(sqrt(X)), X.sqrt.inv=matrix(1/sqrt(X))))
else {
X.eigen <- eigen(X, symmetric = TRUE)
@@ -51,11 +52,46 @@
stop("matrix is not positive semidefinite")
sqomega <- sqrt(diag(X.eigen$values))
sqomega.inv <- diag(1/sqrt(X.eigen$values))
+ sqomega.inv2 <- diag(1/X.eigen$values)
V <- X.eigen$vectors
+ X.inv <- V %*% sqomega.inv2 %*% t(V)
X.sqrt <- V %*% sqomega %*% t(V)
X.sqrt.inv <- V %*% sqomega.inv %*% t(V)
- return(list(X.det=prod(X.eigen$values),
+ return(list(X.det=prod(X.eigen$values),
+ X.inv=X.inv,
X.sqrt=X.sqrt, X.sqrt.inv=X.sqrt.inv))
}
}
+
+invMatrix <- function (X)
+{
+###########################################
+##
+## R-function: invMatrix - computes the inverse of matrix 'X',
+## based on R-function 'root.matrix' of package 'strucchange'
+## author: Bernhard Spangl, based on work of Achim Zeileis
+## version: 0.1 (2025-02-06)
+##
+###########################################
+
+## Paramters:
+## X ... symmetric and positive semidefinite matrix
+
+ if ((ncol(X) == 1) && (nrow(X) == 1))
+ return(list(X.det=X,
+ X.inv=matrix(1/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/X.eigen$values)
+ V <- X.eigen$vectors
+ #X.sqrt <- V %*% sqomega %*% t(V)
+ X.inv <- V %*% sqomega.inv %*% t(V)
+ return(list(X.det=prod(X.eigen$values),
+ X.inv=X.inv))
+ }
+}
+
Modified: pkg/robKalman/R/mACMfilter.R
===================================================================
--- pkg/robKalman/R/mACMfilter.R 2025-05-15 20:35:54 UTC (rev 91)
+++ pkg/robKalman/R/mACMfilter.R 2025-05-17 18:39:59 UTC (rev 92)
@@ -17,6 +17,7 @@
S1 - K %*% W %*% t(K)
}
+
.mACMcorrstep <- function (y, x1, S1, Z, V, rob1=NULL,
psi, apsi, bpsi, cpsi, flag, ...)
{
@@ -65,3 +66,56 @@
return(list(x0 = x0, K = K, S0 = S0, Ind=ind, rob0=st, Delta =NULL, DeltaY = yDelta))
}
+
+.mACM2corrstep <- function (y, x1, S1, Z, V, rob1=NULL,
+ psi, apsi, bpsi, cpsi, flag="weights", ...)
+{
+###########################################
+##
+## R-function: .mACM2corrstep - correction step (internal function)
+## author: Bernhard Spangl
+## version: 0.1 (2025-02-06)
+##
+###########################################
+
+## Paramters:
+## y ... observed vector-valued time series
+## x1 ... state vector (one-step-ahead predictor)
+## S1 ... prediction error covariance matrix (formerly: M)
+## Z ... observation matrix (formerly: H)
+## V ... covariance matrix of observation noise (formerly: R)
+## psi ... influence function to be used
+## a, b, c ... tuning constants for Hampel's psi-function
+## (default: a=b=2.5, c=5.0)
+## flag ... character, either "psi", "weights" or "deriv" to use
+## psi-function (default), weight function psi(t)/t or
+## its derivative respectively
+
+ D <- S1 %*% t(Z)
+ Rt <- Z %*% D + V
+ sqrtR <- rootMatrix(Rt)
+ st <- sqrtR$X.sqrt.inv
+ Rt.inv <- sqrtR$X.inv
+
+ K2 <- D %*% st
+ K <- D %*% Rt.inv
+
+ yDelta <- drop(y - Z %*% x1)
+ yDelta.norm <- as.vector(t(yDelta) %*% Rt.inv %*% yDelta)
+
+ yPsi <- psi(sqrt(yDelta.norm),
+ apsi, bpsi, cpsi, flag)
+ xDelta <- K %*% (as.vector(yPsi)*yDelta)
+ x0 <- x1 + xDelta
+
+ ind <- sqrt(t(yDelta) %*% Rt.inv %*% yDelta) > apsi
+
+ #jacobian.psi <- .weighting(flag)
+ #W <- jacobian.psi(yDelta, a=apsi, b=bpsi, c=cpsi)
+ W <- as.vector(yPsi)*diag(rep(1, length(yDelta)))
+
+ S0 <- .getcorrCovmACM(S1, K2, W = W)
+
+ return(list(x0 = x0, K = K, S0 = S0, Ind=ind, rob0=st, Delta =NULL, DeltaY = yDelta))
+}
+
More information about the Robkalman-commits
mailing list