[Robkalman-commits] r84 - pkg/robKalman/R

[noreply at r-forge.r-project.org]

Author: ruckdeschel
Date: 2025-05-15 16:35:58 +0200 (Thu, 15 May 2025)
New Revision: 84

Modified:
   pkg/robKalman/R/calibrateRLS.R
Log:
calibrateRLS neue Fassung

Modified: pkg/robKalman/R/calibrateRLS.R
===================================================================
--- pkg/robKalman/R/calibrateRLS.R	2025-05-15 14:34:21 UTC (rev 83)
+++ pkg/robKalman/R/calibrateRLS.R	2025-05-15 14:35:58 UTC (rev 84)
@@ -1,5 +1,5 @@
 rLScalibrateB <- function(Z, S, V, repl = 100000, b = NULL, eff = NULL, r = NULL,
-                          rlow = 0, rup = NULL, upto = 20, IO = FALSE, seed)#
+                          rlow = 0, rup = NULL, upto = 20, IO = FALSE, seed = NULL)#
 # calibrates clipping height b to given Z, V, and S_{t|t-1} 
 # --- 
 #      either to given efficiency in the ideal model
@@ -16,12 +16,13 @@
  qd <- ifelse(length(Z)==1, 1, (dim(Z))[1])
  pd <- ifelse(length(S)==1, 1, (dim(Z))[2])
  
- if(!missing(seed)) set.seed(seed)
- 
+ if(!is.null(seed)) set.seed(seed)
+
  dx <- t(mvrnorm(repl, numeric(pd), S))
  eps <- t(mvrnorm(repl, numeric(qd), V))
  dy <- Z %*% dx + eps
- K  <- .getKG(S, Z, .getDelta(S, Z, V))
+ gD <- .getDelta(S, Z, V)
+ K  <- .getKG(S, Z, gD)
 
  trS <- sum(diag(.getcorrCov(S, K, Z)))
 
@@ -29,14 +30,22 @@
 
  if(IO){
     dx0 <- dy - Z %*% dx0
- }
+    Zs0 <- S %*% t(Z)
+    Zs1 <- Z%*%Zs0
+    Zsm <- ginv(Zs1)
+    Zs <- Zs0 %*% Zsm
+    piq <- diag(qd)-Zs1%*%Zsm
+    Kq <- Zs0 %*% piq %*% gD
+    dxq <- t(Kq %*% dy)
+ }else dxq <- 0 * t(dx0)
+ 
  no  <- sqrt(colSums(dx0^2))
 
  eff.b <- function(b){
           w <- ifelse(no < b, 1, b / no)
           dxw  <- as.vector(w) * t(dx0)
-          if(IO) dxw <-  (t(dy) - dxw) %*% t(ginv(Z))
-          trSb <- sum( (t(dx) - dxw)^2 )/repl
+          if(IO) dxw <-  (t(dy) - dxw) %*% t(Zs)
+          trSb <- sum( (t(dx) - dxw - dxq)^2 )/repl
           trS / trSb
  }
  r.b <- function(b){
@@ -58,8 +67,8 @@
                          eff0 = eff, trS0 = trS, repl0 = repl, dY = dy){
               w <- ifelse(no0 < b, 1, b/no0)
               dxw <- as.vector(w) * t(dX0)
-              if(IO) dxw <-  (t(dY) - dxw) %*% t(ginv(Z))
-              trSb <- sum( (t(dX) - dxw)^2 )/repl0
+              if(IO) dxw <-  (t(dY) - dxw) %*% t(Zs)
+              trSb <- sum( (t(dX) - dxw - dxq)^2 )/repl0
               trS0 / trSb - eff0
           }
           r1 <- NULL