[Robast-commits] r1283 - pkg/ROptEst/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Feb 6 22:54:09 CET 2024
Author: ruckdeschel
Date: 2024-02-06 22:54:09 +0100 (Tue, 06 Feb 2024)
New Revision: 1283
Modified:
pkg/ROptEst/R/LowerCaseMultivariate.R
pkg/ROptEst/R/getAsRisk.R
pkg/ROptEst/R/getInfRobIC_asBias.R
Log:
[ROptEst] [trunk] revised code for lower case solutions
(a) in univariate case for asymmetric bias -> in particular minmaxBias()
(b) in multivariate case to make use of symmetry
Modified: pkg/ROptEst/R/LowerCaseMultivariate.R
===================================================================
--- pkg/ROptEst/R/LowerCaseMultivariate.R 2024-02-06 20:29:39 UTC (rev 1282)
+++ pkg/ROptEst/R/LowerCaseMultivariate.R 2024-02-06 21:54:09 UTC (rev 1283)
@@ -19,14 +19,13 @@
z.comp <- rep(TRUE, nrow(trafo))
force(normtype)
- lA.comp <- sum(A.comp)
A.symm <- (nrow(trafo)==ncol(trafo)) && isTRUE(all.equal(trafo,t(trafo)))
if(A.symm){
A.comp.s <- t(A.comp)|A.comp
- A.comp <- A.com.s[col(A.com.s)>=row(A.com.s)]
- lA.comp <- sum(A.comp.so)
+ A.comp <- A.comp.s[col(A.comp.s)>=row(A.comp.s)]
}
+ lA.comp <- sum(A.comp)
abs.fct <- function(x, L2, stand, cent, normtype.0){
X <- evalRandVar(L2, as.matrix(x))[,,1] - cent
@@ -41,7 +40,7 @@
k <- ncol(trafo)
A <- matrix(0, ncol = k, nrow = p)
- A[A.comp] <- param[1:1A.comp]
+ A[A.comp] <- param[1:lA.comp]
if(A.symm) A[col(A)>row(A)] <- t(A)[col(A)>row(A)]
A.max <- max(abs(A.comp))
A <- A/A.max
Modified: pkg/ROptEst/R/getAsRisk.R
===================================================================
--- pkg/ROptEst/R/getAsRisk.R 2024-02-06 20:29:39 UTC (rev 1282)
+++ pkg/ROptEst/R/getAsRisk.R 2024-02-06 21:54:09 UTC (rev 1283)
@@ -81,8 +81,12 @@
comp <- .getComp(L2deriv, DistrSymm, L2derivSymm, L2derivDistrSymm)
z.comp <- comp$"z.comp"
A.comp <- comp$"A.comp"
- DA.comp <- abs(trafo) %*% A.comp != 0
+ DA.comp <- matrix(TRUE, nrow=nrow(trafo), ncol=ncol(trafo))
+ DA.symm <- (nrow(trafo)==ncol(trafo)) && isTRUE(all.equal(trafo,t(trafo)))
+ if(DA.symm) DA.comp <- A.comp
+
+
eerg <- .LowerCaseMultivariate(L2deriv = L2deriv, neighbor = neighbor,
biastype = biastype, normtype = normtype, Distr = Distr, Finfo = Finfo,
trafo = trafo, z.start = z.start, A.start = A.start, z.comp = z.comp,
@@ -92,6 +96,7 @@
return(list(asBias = bias, normtype = eerg$normtype))
})
+
setMethod("getAsRisk", signature(risk = "asBias",
L2deriv = "RealRandVariable",
neighbor = "TotalVarNeighborhood",
@@ -294,17 +299,15 @@
return(list(asBias = 0, warn = gettext("not attained by IC")))
sign <- sign(biastype)
- w0 <- options("warn")
- on.exit(options(w0))
- options(warn = -1)
- l <- length(support(L2deriv))
- if (sign>0)
- {z0 <- support(L2deriv)[1]; deltahat <- support(L2deriv)[2]-z0}
- else
- {z0 <- support(L2deriv)[l]; deltahat <- z0-support(L2deriv)[l-1]}
-
- bias <- abs(as.vector(trafo))/abs(z0)
+ if(missing(trafo)) trafo <- 1
+ if(length(trafo)>1) stop("Matrix/vector-valued 'trafo' is not (yet) supported.")
+ trafo <- c( trafo )
+ sign <- sign*sign(trafo)
+
+ z0 <- if (sign>0) min(support(L2deriv)) else max(support(L2deriv))
+
+ bias <- abs(trafo)/abs(z0)
return(list(asBias = bias))
})
Modified: pkg/ROptEst/R/getInfRobIC_asBias.R
===================================================================
--- pkg/ROptEst/R/getInfRobIC_asBias.R 2024-02-06 20:29:39 UTC (rev 1282)
+++ pkg/ROptEst/R/getInfRobIC_asBias.R 2024-02-06 21:54:09 UTC (rev 1283)
@@ -202,8 +202,16 @@
if(missing(verbose)|| is.null(verbose))
verbose <- getRobAStBaseOption("all.verbose")
- DA.comp <- abs(trafo) %*% A.comp != 0
- eerg <- .LowerCaseMultivariate(L2deriv = L2deriv, neighbor = neighbor,
+
+ A.test <- trafo * A.comp
+ A.symm <- (nrow(A.test)==ncol(A.test)) && isTRUE(all.equal(A.test,t(A.test)))
+
+ DA.comp <- (abs(A.test) >= 1e-8 * max(abs(A.test)))
+
+ if(A.symm) DA.comp[col(DA.comp)>row(DA.comp)] <- FALSE
+
+
+ eerg <- .LowerCaseMultivariate(L2deriv = L2deriv, neighbor = neighbor,
biastype = biastype, normtype = normtype, Distr = Distr,
Finfo = Finfo, trafo, z.start, A.start = A.start, z.comp = z.comp,
A.comp = DA.comp, maxiter = maxiter, tol = tol, verbose = verbose, ...)
@@ -216,8 +224,9 @@
p <- nrow(trafo)
k <- ncol(trafo)
A <- matrix(0, ncol=k, nrow=p)
- A[DA.comp] <- matrix(param[1:lA.comp], ncol=k, nrow=p)
- A.max <- max(abs(A))
+ A[DA.comp] <- param[1:lA.comp]
+ if(A.symm) A[col(A)>row(A)] <- t(A)[col(A)>row(A)]
+ A.max <- max(abs(A))
A <- A/A.max
z <- numeric(k)
z[z.comp] <- param[(lA.comp+1):length(param)]
@@ -394,28 +403,24 @@
{ if(is.finite(lowerCaseRadius(L2deriv, neighbor, risk = asMSE(), biastype)))
{
sign <- sign(biastype)
- w0 <- options("warn")
- on.exit(options(w0))
- options(warn = -1)
- l <- length(support(L2deriv))
- if (sign>0)
- {z0 <- support(L2deriv)[1]
- deltahat <- support(L2deriv)[2]-z0
- }else{
- z0 <- support(L2deriv)[l]
- deltahat <- z0-support(L2deriv)[l-1]
- }
+ supp.s <- sort(support(L2deriv))
+ if(sign < 0) supp.s <- rev(supp.s)
+ l <- length(supp.s)
+
+ z0 <- supp.s[1]
+ z1 <- supp.s[2]
+
p0 <- d(L2deriv)(z0)
- v1 <- (1-p0)/p0/z0
- v2 <- -1/z0
- c0 <- deltahat*p0/2
- A0 <- abs(1/z0/c0)
- zc <- z0+sign(biastype)*deltahat*(1-p0)/2
+ v0 <- (1-p0)/p0/z0
+ v1 <- -1/z0
+
+ zc <- p0*z0+(1-p0)*z1
+ A0 <- (v1-v0)/(z1-z0)
a0 <- A0*zc
b0 <- abs(1/z0)
d0 <- 0
- asCov <- v1^2*(1-p0)+v2^2*p0
+ asCov <- v0^2*p0+v1^2*(1-p0)
Risk0 <- list(asBias = list(value = b0, biastype = biastype,
normtype = NormType(),
neighbortype = class(neighbor)),
@@ -423,7 +428,7 @@
A0 <- matrix(A0,1,1)
w <- new("HampelWeight")
- cent(w) <- z0
+ cent(w) <- zc
stand(w) <- A0
clip(w) <- b0
weight(w) <- minbiasweight(w, neighbor = neighbor,
More information about the Robast-commits
mailing list