[Robast-commits] r337 - in branches/robast-0.7/pkg: ROptEst ROptEst/R ROptEst/chm ROptEst/man RobAStBase/R RobAStBase/chm RobAStBase/man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Aug 12 03:04:02 CEST 2009
Author: ruckdeschel
Date: 2009-08-12 03:04:02 +0200 (Wed, 12 Aug 2009)
New Revision: 337
Modified:
branches/robast-0.7/pkg/ROptEst/NAMESPACE
branches/robast-0.7/pkg/ROptEst/R/getFixRobIC_fiUnOvShoot.R
branches/robast-0.7/pkg/ROptEst/R/getIneffDiff.R
branches/robast-0.7/pkg/ROptEst/R/getInfCent.R
branches/robast-0.7/pkg/ROptEst/R/getInfClip.R
branches/robast-0.7/pkg/ROptEst/R/getInfGamma.R
branches/robast-0.7/pkg/ROptEst/R/getInfRobIC_asGRisk.R
branches/robast-0.7/pkg/ROptEst/R/getInfRobIC_asHampel.R
branches/robast-0.7/pkg/ROptEst/R/getInfRobIC_asUnOvShoot.R
branches/robast-0.7/pkg/ROptEst/R/getInfStand.R
branches/robast-0.7/pkg/ROptEst/R/getInfV.R
branches/robast-0.7/pkg/ROptEst/R/optIC.R
branches/robast-0.7/pkg/ROptEst/R/radiusMinimaxIC.R
branches/robast-0.7/pkg/ROptEst/R/roptest.R
branches/robast-0.7/pkg/ROptEst/chm/00Index.html
branches/robast-0.7/pkg/ROptEst/chm/ROptEst.chm
branches/robast-0.7/pkg/ROptEst/chm/ROptEst.hhp
branches/robast-0.7/pkg/ROptEst/chm/ROptEst.toc
branches/robast-0.7/pkg/ROptEst/chm/getFixRobIC.html
branches/robast-0.7/pkg/ROptEst/chm/getIneffDiff.html
branches/robast-0.7/pkg/ROptEst/chm/getInfCent.html
branches/robast-0.7/pkg/ROptEst/chm/getInfClip.html
branches/robast-0.7/pkg/ROptEst/chm/getInfGamma.html
branches/robast-0.7/pkg/ROptEst/chm/getInfRobIC.html
branches/robast-0.7/pkg/ROptEst/chm/getInfStand.html
branches/robast-0.7/pkg/ROptEst/chm/getInfV.html
branches/robast-0.7/pkg/ROptEst/chm/optIC.html
branches/robast-0.7/pkg/ROptEst/chm/radiusMinimaxIC.html
branches/robast-0.7/pkg/ROptEst/chm/roptest.html
branches/robast-0.7/pkg/ROptEst/man/getFixRobIC.Rd
branches/robast-0.7/pkg/ROptEst/man/getIneffDiff.Rd
branches/robast-0.7/pkg/ROptEst/man/getInfCent.Rd
branches/robast-0.7/pkg/ROptEst/man/getInfClip.Rd
branches/robast-0.7/pkg/ROptEst/man/getInfGamma.Rd
branches/robast-0.7/pkg/ROptEst/man/getInfRobIC.Rd
branches/robast-0.7/pkg/ROptEst/man/getInfStand.Rd
branches/robast-0.7/pkg/ROptEst/man/getInfV.Rd
branches/robast-0.7/pkg/ROptEst/man/optIC.Rd
branches/robast-0.7/pkg/ROptEst/man/radiusMinimaxIC.Rd
branches/robast-0.7/pkg/RobAStBase/R/TotalVarIC.R
branches/robast-0.7/pkg/RobAStBase/R/bALEstimate.R
branches/robast-0.7/pkg/RobAStBase/chm/00Index.html
branches/robast-0.7/pkg/RobAStBase/chm/RobAStBase.chm
branches/robast-0.7/pkg/RobAStBase/chm/RobAStBase.toc
branches/robast-0.7/pkg/RobAStBase/chm/TotalVarIC-class.html
branches/robast-0.7/pkg/RobAStBase/man/TotalVarIC-class.Rd
Log:
--------------------------------------------------------------------------
multivariate solution Hampel/G-Risk --- new 10-08-09
--------------------------------------------------------------------------
getInfRob_asGRisk.R now have both an iterative and an optimization
getInfRob_asHampel.R algorithm to obtain Lagrange Multipliers;
--- i.e. depending on argument OptOrIter (default: iteration as before...)
they call getLagrangeMultByOptim resp. getLagrangeMultByIter
(unfortunately: optimization is slow; only L-BFGS-B works satisfactory)
+several steps outsourced in internal, non-exported functions
.getLowUpB() .getLowerSol() .getUpperSol() .checkUpLow()
+remains to be done: self standardization does not yet work...
+total variation nbd: case 1=p<k works
but still some checking needed.
+to be checked: symm for total variation
+getInfCent uses getLow/getUp now
Modified: branches/robast-0.7/pkg/ROptEst/NAMESPACE
===================================================================
--- branches/robast-0.7/pkg/ROptEst/NAMESPACE 2009-08-11 02:00:54 UTC (rev 336)
+++ branches/robast-0.7/pkg/ROptEst/NAMESPACE 2009-08-12 01:04:02 UTC (rev 337)
@@ -27,4 +27,4 @@
"cniperCont", "cniperPoint", "cniperPointPlot")
exportMethods("updateNorm")
export("getL2normL2deriv")
-export("roptest")
+export("roptest","getLagrangeMultByOptim","getLagrangeMultByIter")
Modified: branches/robast-0.7/pkg/ROptEst/R/getFixRobIC_fiUnOvShoot.R
===================================================================
--- branches/robast-0.7/pkg/ROptEst/R/getFixRobIC_fiUnOvShoot.R 2009-08-11 02:00:54 UTC (rev 336)
+++ branches/robast-0.7/pkg/ROptEst/R/getFixRobIC_fiUnOvShoot.R 2009-08-12 01:04:02 UTC (rev 337)
@@ -4,7 +4,7 @@
setMethod("getFixRobIC", signature(Distr = "Norm",
risk = "fiUnOvShoot",
neighbor = "UncondNeighborhood"),
- function(Distr, risk, neighbor, sampleSize, upper, maxiter, tol, warn,
+ function(Distr, risk, neighbor, sampleSize, upper, lower, maxiter, tol, warn,
Algo, cont){
radius <- neighbor at radius
if(identical(all.equal(radius, 0), TRUE)){
Modified: branches/robast-0.7/pkg/ROptEst/R/getIneffDiff.R
===================================================================
--- branches/robast-0.7/pkg/ROptEst/R/getIneffDiff.R 2009-08-11 02:00:54 UTC (rev 336)
+++ branches/robast-0.7/pkg/ROptEst/R/getIneffDiff.R 2009-08-12 01:04:02 UTC (rev 337)
@@ -6,14 +6,15 @@
neighbor = "UncondNeighborhood",
risk = "asMSE"),
function(radius, L2Fam, neighbor, risk, loRad, upRad, loRisk, upRisk,
- z.start = NULL, A.start = NULL, upper.b = NULL, MaxIter, eps, warn,
- loNorm = NULL, upNorm = NULL, verbose = FALSE){
+ z.start = NULL, A.start = NULL, upper.b = NULL, lower.b = NULL,
+ MaxIter, eps, warn,
+ loNorm = NULL, upNorm = NULL, verbose = FALSE, ...){
L2derivDim <- numberOfMaps(L2Fam at L2deriv)
if(L2derivDim == 1){
neighbor at radius <- radius
res <- getInfRobIC(L2deriv = L2Fam at L2derivDistr[[1]], neighbor = neighbor,
risk = risk, symm = L2Fam at L2derivDistrSymm[[1]],
- Finfo = L2Fam at FisherInfo, upper = upper.b,
+ Finfo = L2Fam at FisherInfo, upper = upper.b, lower = lower.b,
trafo = trafo(L2Fam at param), maxiter = MaxIter, tol = eps,
warn = warn, verbose = verbose)
trafo <- as.vector(trafo(L2Fam at param))
@@ -52,10 +53,12 @@
Distr = L2Fam at distribution, DistrSymm = L2Fam at distrSymm,
L2derivSymm = L2derivSymm, L2derivDistrSymm = L2derivDistrSymm,
Finfo = L2Fam at FisherInfo, trafo = trafo, z.start = z.start,
- A.start = A.start, upper = upper.b, maxiter = MaxIter,
- tol = eps, warn = warn, verbose = verbose)
+ A.start = A.start, upper = upper.b, lower = lower.b,
+ maxiter = MaxIter,
+ tol = eps, warn = warn, verbose = verbose, ...)
normtype(risk) <- res$normtype
- std <- if(is(normtype(risk),"QFNorm")) QuadForm(normtype(risk)) else diag(p)
+ std <- if(is(normtype(risk),"QFNorm"))
+ QuadForm(normtype(risk)) else diag(p)
biasLo <- biasUp <- res$b
Modified: branches/robast-0.7/pkg/ROptEst/R/getInfCent.R
===================================================================
--- branches/robast-0.7/pkg/ROptEst/R/getInfCent.R 2009-08-11 02:00:54 UTC (rev 336)
+++ branches/robast-0.7/pkg/ROptEst/R/getInfCent.R 2009-08-12 01:04:02 UTC (rev 337)
@@ -11,8 +11,8 @@
z.fct <- function(z, c0, D1){
return(c0 + (z-c0)*p(D1)(z-c0) - (z+c0)*p(D1)(z+c0) + m1df(D1, z+c0) - m1df(D1, z-c0))
}
- lower <- q(L2deriv)(getdistrOption("TruncQuantile"))
- upper <- q(L2deriv)(1-getdistrOption("TruncQuantile"))
+ lower <- getLow(L2deriv)
+ upper <- getUp(L2deriv)
return(uniroot(z.fct, lower = lower, upper = upper, tol = tol.z,
c0=clip, D1=L2deriv)$root)
@@ -26,19 +26,46 @@
D1 <- sign(as.vector(trafo))*L2deriv
g.fct <- function(g, c0, D1){
- return(g*p(D1)(g) + (g+c0)*(1-p(D1)(g+c0)) - m1df(D1, g) + m1df(D1, g+c0))
+ return(g*p(D1)(g) + (g+c0)*(p(D1)(g+c0, lower.tail = FALSE)) - m1df(D1, g) + m1df(D1, g+c0))
}
- lower <- q(L2deriv)(getdistrOption("TruncQuantile"))
- upper <- q(L2deriv)(1-getdistrOption("TruncQuantile"))
+ lower <- getLow(L2deriv)
+ upper <- getUp(L2deriv)
return(uniroot(g.fct, lower = lower, upper = upper, tol = tol.z,
c0 = clip, D1 = D1)$root)
})
+
setMethod("getInfCent", signature(L2deriv = "RealRandVariable",
+ neighbor = "TotalVarNeighborhood",
+ biastype = "BiasType"),
+ function(L2deriv, neighbor, biastype, Distr, z.comp, w,
+ tol.z = .Machine$double.eps^.5){
+ stand <- stand(w)
+ clip <- clip(w)
+ b <- clip[2]-clip[1]
+ ### if(symm) return(b/2)
+
+ g.fct <- function(g, c0, D1){
+ fct <- function(x){
+ Lx <- evalRandVar(L2deriv, as.matrix(x)) [,,1]
+ D1 <- as.numeric(stand%*%Lx)
+ pmin(pmax(g,D1),g+c0)
+ }
+ return(E(object = Distr, fun = fct, useApply = FALSE))
+ }
+ lower <- -b
+ upper <- 0
+
+ return(uniroot(g.fct, lower = lower, upper = upper, tol = tol.z,
+ c0 = b, D1 = D1)$root)
+ })
+
+setMethod("getInfCent", signature(L2deriv = "RealRandVariable",
neighbor = "ContNeighborhood",
biastype = "BiasType"),
- function(L2deriv, neighbor, biastype, Distr, z.comp, w){
+ function(L2deriv, neighbor, biastype, Distr, z.comp, w,
+ tol.z = .Machine$double.eps^.5){
integrand1 <- function(x){
weight(w)(evalRandVar(L2deriv, as.matrix(x)) [,,1])
}
@@ -71,14 +98,14 @@
z.fct <- function(z, c0, D1){
return(c0 - (z+c0)*p(D1)(z+c0) + m1df(D1, z+c0))
}
- lower <- q(L2deriv)(getdistrOption("TruncQuantile"))
+ lower <- getLow(L2deriv)
upper <- 0
}else{
z.fct <- function(z, c0, D1){
return(- z + (z-c0)*p(D1)(z-c0) - m1df(D1, z-c0))
}
lower <- 0
- upper <- q(L2deriv)(1-getdistrOption("TruncQuantile"))
+ upper <- getUp(L2deriv)
}
return(uniroot(z.fct, lower = lower, upper = upper, tol = tol.z,
c0=clip, D1=L2deriv)$root)
@@ -96,8 +123,8 @@
(z+c0/nu2)*p(D1)(z+c0/nu2) + m1df(D1, z+c0/nu2) -
m1df(D1, z-c0/nu1))
}
- lower <- q(L2deriv)(getdistrOption("TruncQuantile"))
- upper <- q(L2deriv)(1-getdistrOption("TruncQuantile"))
+ lower <- getLow(L2deriv)
+ upper <- getUp(L2deriv)
return(uniroot(z.fct, lower = lower, upper = upper, tol = tol.z,
c0=clip, D1=L2deriv)$root)
Modified: branches/robast-0.7/pkg/ROptEst/R/getInfClip.R
===================================================================
--- branches/robast-0.7/pkg/ROptEst/R/getInfClip.R 2009-08-11 02:00:54 UTC (rev 336)
+++ branches/robast-0.7/pkg/ROptEst/R/getInfClip.R 2009-08-12 01:04:02 UTC (rev 337)
@@ -30,7 +30,7 @@
setMethod("getInfClip", signature(clip = "numeric",
L2deriv = "EuclRandVariable",
risk = "asMSE",
- neighbor = "ContNeighborhood"),
+ neighbor = "UncondNeighborhood"),
function(clip, L2deriv, risk, neighbor, biastype,
Distr, stand, cent, trafo){
return(neighbor at radius^2*clip +
Modified: branches/robast-0.7/pkg/ROptEst/R/getInfGamma.R
===================================================================
--- branches/robast-0.7/pkg/ROptEst/R/getInfGamma.R 2009-08-11 02:00:54 UTC (rev 336)
+++ branches/robast-0.7/pkg/ROptEst/R/getInfGamma.R 2009-08-12 01:04:02 UTC (rev 337)
@@ -28,19 +28,36 @@
neighbor = "ContNeighborhood",
biastype = "BiasType"),
function(L2deriv, risk, neighbor, biastype, Distr,
- stand, cent, clip){
- integrandG <- function(x, L2, stand, cent, clip){
+ stand, cent, clip, power = 1L){
+ integrandG <- function(x, L2, stand, cent, clip){
X <- evalRandVar(L2, as.matrix(x))[,,1] - cent
Y <- stand %*% X
res <- norm(risk)(Y) - clip
- return((res > 0)*res)
+ return((res > 0)*res^power)
}
return(-E(object = Distr, fun = integrandG, L2 = L2deriv,
stand = stand, cent = cent, clip = clip, useApply = FALSE))
})
+setMethod("getInfGamma", signature(L2deriv = "RealRandVariable",
+ risk = "asMSE",
+ neighbor = "TotalVarNeighborhood",
+ biastype = "BiasType"),
+ function(L2deriv, risk, neighbor, biastype, Distr,
+ stand, cent, clip, power = 1L){
+ integrandG <- function(x, L2, stand, cent, clip){
+ X <- evalRandVar(L2, as.matrix(x))[,,1] - cent
+ Y <- stand %*% X
+ res <- Y - clip
+
+ return((res > 0)*res^power)
+ }
+
+ return(-E(object = Distr, fun = integrandG, L2 = L2deriv,
+ stand = stand, cent = cent, clip = clip, useApply = FALSE))
+ })
###############################################################################
## gamma in case of asymptotic under-/overshoot risk
###############################################################################
Modified: branches/robast-0.7/pkg/ROptEst/R/getInfRobIC_asGRisk.R
===================================================================
--- branches/robast-0.7/pkg/ROptEst/R/getInfRobIC_asGRisk.R 2009-08-11 02:00:54 UTC (rev 336)
+++ branches/robast-0.7/pkg/ROptEst/R/getInfRobIC_asGRisk.R 2009-08-12 01:04:02 UTC (rev 337)
@@ -4,7 +4,8 @@
setMethod("getInfRobIC", signature(L2deriv = "UnivariateDistribution",
risk = "asGRisk",
neighbor = "UncondNeighborhood"),
- function(L2deriv, risk, neighbor, symm, Finfo, trafo, upper = NULL, maxiter, tol,
+ function(L2deriv, risk, neighbor, symm, Finfo, trafo, upper = NULL,
+ lower = NULL, maxiter, tol,
warn, noLow = FALSE, verbose = FALSE){
biastype <- biastype(risk)
radius <- neighbor at radius
@@ -46,12 +47,16 @@
## new
L1n <- getL1normL2deriv(L2deriv = L2deriv, cent = z)
lower0 <- L1n/(1 + radius^2)
- if(!is(neighbor,"ContNeighborhood")) lower0 <- lower0/2
+ if(is(neighbor,"TotalVarNeighborhood")) {
+ lower0 <- (L1n-z)/(1 + radius^2)/2}
upper0 <- max(L1n/radius,
sqrt( as.numeric( Finfo + z^2 )/(( 1 + radius^2)^2 - 1) ))
- if (!is.null(upper)|(iter == 1))
- {lower <- .Machine$double.eps^0.6
- }else{ lower <- lower0; upper <- upper0}
+ if (is.null(lower)|(iter == 1))
+ lower <- .Machine$double.eps^0.6
+ else {if(iter>1) lower <- max(lower0,lower)}
+ if (is.null(upper)|(iter == 1))
+ upper <- 5* max(abs(trafo))*max(Finfo)
+ else {if(iter>1) upper <- min(upper,upper0)}
##
c0 <- try(uniroot(getInfClip,
## new
@@ -123,9 +128,9 @@
if(is(neighbor,"ContNeighborhood")){
w <- new("HampelWeight")
clip(w) <- b
- cent(w) <- z
+ cent(w) <- as.numeric(z)
stand(w) <- A
- }else{
+ }else if (is(neighbor,"TotalVarNeighborhood")){
w <- new("BdStWeight")
clip(w) <- c(0,b)+a
stand(w) <- A
@@ -140,200 +145,621 @@
-################################################################################
+###################################################################################
+# multivariate solution G-Risk --- new 10-08-09
+###################################################################################
+setMethod("getInfRobIC", signature(L2deriv = "RealRandVariable",
+ risk = "asGRisk",
+ neighbor = "UncondNeighborhood"),
+ function(L2deriv, risk, neighbor, Distr, DistrSymm, L2derivSymm,
+ L2derivDistrSymm, Finfo, trafo, onesetLM = FALSE,
+ z.start, A.start, upper = NULL, lower = NULL,
+ OptOrIter = "iterate",
+ maxiter, tol, warn, verbose = FALSE, ...){
-setMethod("getInfRobIC", signature(L2deriv = "RealRandVariable",
- risk = "asGRisk",
- neighbor = "UncondNeighborhood"),
- function(L2deriv, risk, neighbor, Distr, DistrSymm, L2derivSymm,
- L2derivDistrSymm, Finfo, trafo, onesetLM = FALSE,
- z.start, A.start, upper = NULL, maxiter, tol, warn, verbose = FALSE){
+ mc <- match.call()
+
+ ## some abbreviations / checks
+ radius <- neighbor at radius
biastype <- biastype(risk)
normtype <- normtype(risk)
+
p <- nrow(trafo)
- if(! is(neighbor,"ContNeighborhood") && p>1)
+ k <- ncol(trafo)
+
+ if( is(neighbor,"TotalVarNeighborhood") && p>1)
stop("Not yet implemented")
+ ## non-standard norms
FI <- solve(trafo%*%solve(Finfo)%*%t(trafo))
- if(is(normtype,"InfoNorm") || is(normtype,"SelfNorm") )
- {QuadForm(normtype) <- PosSemDefSymmMatrix(FI);
- normtype(risk) <- normtype}
- QF <- if(is(normtype,"QFNorm")) QuadForm(normtype) else diag(nrow(trafo))
+ if(is(normtype,"InfoNorm") || is(normtype,"SelfNorm") ){
+ QuadForm(normtype) <- PosSemDefSymmMatrix(FI)
+ normtype(risk) <- normtype
+ }
+ std <- if(is(normtype,"QFNorm"))
+ QuadForm(normtype) else diag(p)
- if(is.null(z.start)) z.start <- numeric(ncol(trafo))
+ ## starting values
+ if(is.null(z.start)) z.start <- numeric(k)
if(is.null(A.start)) A.start <- trafo %*% solve(Finfo)
- radius <- neighbor at radius
- if(identical(all.equal(radius, 0), TRUE)){
- if(warn) cat("'radius == 0' => (classical) optimal IC\n",
- "in sense of Cramer-Rao bound is returned\n")
- res <- getInfRobIC(L2deriv = L2deriv, risk = asCov(), neighbor = neighbor,
- Distr = Distr, Finfo = Finfo, trafo = trafo,
- QuadForm = QF, verbose = verbose)
- res <- c(res, list(biastype = biastype, normtype = normtype))
- if(!is(risk, "asMSE")){
- Risk <- getAsRisk(risk = risk, L2deriv = L2deriv, neighbor = neighbor,
- biastype = biastype, cent = res$a,
- stand = res$A, trafo = trafo)
- res$risk <- c(Risk, res$risk)
- }
- trAsCov <- sum(diag(QF%*%res$risk$asCov));
- res$risk$trAsCov <- list(value = trAsCov, normtype = normtype)
- res$risk$asBias <- list(value = b, biastype = biastype,
- normtype = normtype,
- neighbortype = class(neighbor))
- res$risk$asMSE <- list(value = trAsCov + radius^2*b^2,
- r = radius,
- at = neighbor)
- return(res)
- }
+ ## sort out upper solution if radius = 0
+ if(identical(all.equal(radius, 0), TRUE))
+ return(.getUpperSol(L2deriv = L2deriv, b = b, radius = radius,
+ risk = risk, neighbor = neighbor,
+ biastype = biastype, normtype = normtype,
+ Distr = Distr, Finfo = Finfo, trafo = trafo,
+ QF = std, verbose = verbose, warn = warn))
+
+ ## determine which entries must be computed
comp <- .getComp(L2deriv, DistrSymm, L2derivSymm, L2derivDistrSymm)
-
z.comp <- comp$"z.comp"
A.comp <- comp$"A.comp"
+ ## selection of the algorithm
+ pM <- pmatch(tolower(OptOrIter),c("optimize","iterate", "doubleiterate"))
+ OptOrIter <- pM
+ if (is.na(pM)) OptOrIter <- 1
+
+ ## initialize
if(is(neighbor,"ContNeighborhood")){
w <- new("HampelWeight")
- }else{
+ }else{ if (is(neighbor,"TotalVarNeighborhood"))
w <- new("BdStWeight")
}
z <- z.start
A <- A.start
b <- 0
+ a <- 0
iter <- 0
prec <- 1
- repeat{
- iter <- iter + 1
- z.old <- z
- b.old <- b
- A.old <- A
- ##
- if(is(neighbor,"ContNeighborhood"))
- cent(w) <- z
- stand(w) <- A
+ iter.In <- 0
- ## new
- L1n <- getL1normL2deriv(L2deriv = L2deriv, cent = z, stand = A,
- Distr = Distr, normtype = normtype)
- lower0 <- L1n/(1+radius^2)
+ if(is(neighbor,"ContNeighborhood")){
+ w <- new("HampelWeight")
+ }else if(is(neighbor,"TotalVarNeighborhood")){
+ w <- new("BdStWeight")
+ }
- QF <- if(is(normtype,"QFNorm")) QuadForm(normtype) else diag(nrow(A))
- upper0 <- max(L1n/radius,
- sqrt( (sum( diag(QF%*%A%*%Finfo%*%t(A))) + t(A%*%z)%*%QF%*%(A%*%z)) /
- ((1 + radius^2)^2-1)))
+ ## determining A,a,b with either optimization of iteration:
+ if(OptOrIter == 1){
+ if(is.null(lower)){
+ lowBerg <- .getLowerSol(L2deriv = L2deriv, risk = risk,
+ neighbor = neighbor, Distr = Distr,
+ DistrSymm = DistrSymm,
+ L2derivSymm = L2derivSymm,
+ L2derivDistrSymm = L2derivDistrSymm,
+ z.start = z.start, A.start = A.start,
+ trafo = trafo, maxiter = maxiter, tol = tol,
+ warn = FALSE, Finfo = Finfo, verbose = verbose)
+ lower <- lowBerg$b}
+ #if(is.null(upper))
+ upper <- 5*max(solve(Finfo))
- if (!is.null(upper)|(iter == 1))
- {lower <- .Machine$double.eps^0.6;
- if(is.null(upper)) upper <- 10*upper0
- }else{ lower <- lower0; upper <- upper0}
+ iter.In <- 0
+ prec.In <- 0
+ OptIterCall <- numeric(1)
+ Cov <- 0
+ Risk <- 1e10
+ normtype.old <- normtype
+ z.opt <- z
+ A.opt <- A
+ a.opt <- as.numeric(A.opt %*% z.opt)
+ w.opt <- w
+ std.opt <- std
+ normtype.opt <- normtype
- if(!is(neighbor,"ContNeighborhood")) lower0 <- lower0/2
- ##
- b <- try(uniroot(getInfClip,
- ## new
- lower = lower, upper = upper,
- ##
- tol = tol, L2deriv = L2deriv, risk = risk,
- biastype = biastype, Distr = Distr, neighbor = neighbor,
- stand = A, cent = z, trafo = trafo)$root, silent = TRUE)
- if(!is.numeric(b)){
- if(warn) cat("Could not determine optimal clipping bound!\n",
+ asGRiskb <- function(b0){
+ iter <<- iter + 1
+ erg <- getLagrangeMultByOptim(b = b0, L2deriv = L2deriv, risk = risk,
+ FI = Finfo, trafo = trafo, neighbor = neighbor,
+ biastype = biastype, normtype = normtype.opt, Distr = Distr,
+ z.start = z.opt, A.start = A.opt, w.start = w.opt, std = std.opt,
+ z.comp = z.comp, A.comp = A.comp,
+ maxiter = round(maxiter/50*iter^5), tol = tol^(iter^5/40),
+ onesetLM = onesetLM, verbose = verbose, ...)
+
+ w0 <- erg$w
+ A0 <- erg$A
+ a0 <- erg$a
+ z0 <- erg$z
+ std0 <- erg$std
+ biastype0 <- erg$biastype
+ normtype.old0 <- normtype
+ normtype0 <- erg$normtype
+ risk0 <- erg$risk
+ iter.In <<- iter.In + erg$iter
+# cat("Iteration ",iter ,"--innen:", erg$iter,"\n")
+# print(round(c(maxiter = maxiter/30*iter, tol = tol^(iter/15)),4))
+ Cov0 <- getInfV(L2deriv = L2deriv, neighbor = neighbor,
+ biastype = biastype, Distr = Distr,
+ V.comp = A.comp, cent = a0,
+ stand = A0, w = w0)
+
+ if(!is(risk, "asMSE")){
+ Risk0 <- getAsRisk(risk = risk0, L2deriv = L2deriv, neighbor = neighbor,
+ biastype = biastype0, clip = b0, cent = a0, stand = A0,
+ trafo = trafo)
+ }else{
+ Risk0 <- sum(diag(std0%*%Cov0)) + radius^2 * b0^2
+ }
+ Risk <<- Risk0
+# print("A")
+# print(list(Risk=Risk0,A=A0,a=a0,z=z0,b=b0))
+# print("...")
+ w.opt <<- w <<- w0
+ A.opt <<- A <<- A0
+ a.opt <<- a <<- a0
+ z.opt <<- z <<- z0
+ b <<- b0
+ std.opt <<- std <<- std0
+ biastype <<- biastype0
+ normtype.old <<- normtype.old0
+ normtype.opt <<- normtype <<- normtype0
+ risk <<- risk0
+ prec.In <<- erg$prec
+ OptIterCall <<- erg$call
+ Cov <<- Cov0
+
+ return(Risk0-sum(diag(std0%*%A0%*%t(trafo))))
+ }
+ tol0 <- tol^.5
+ f.l <- asGRiskb(lower)
+ f.u <- asGRiskb(upper)
+
+ du <- uniroot(asGRiskb, interval = c(lower,upper), tol = tol0,
+ f.lower=f.l, f.upper=f.u)
+# du <- optimize(asGRiskb, interval = c(lower,upper), tol = tol0,
+# A.o = A.opt, z.o = z.opt, a.o = a.opt, w.o = w.opt,
+# std.o = std.opt, normtype.o = normtype.opt)
+ }else{
+ repeat{
+ iter <- iter + 1
+ z.old <- z
+ b.old <- b
+ A.old <- A
+ ##
+ ## get interval for b-search:
+ LUB <- .getLowUpB(L2deriv = L2deriv, Finfo = Finfo, Distr = Distr,
+ normtype = normtype, z = z, A = A, radius = radius,
+ iter = iter)
+ lower <- LUB$lower
+ upper <- if(is.null(upper)) LUB$upper else min(upper,LUB$upper)
+
+# print(c(lower,upper))
+ ## solve for b
+ b <- try(uniroot(getInfClip,
+ lower = lower, upper = upper,
+ tol = tol, L2deriv = L2deriv, risk = risk,
+ biastype = biastype, Distr = Distr, neighbor = neighbor,
+ stand = A, cent = z, trafo = trafo)$root, silent = TRUE)
+
+ ## if no solution return lower case:
+ if(!is.numeric(b))
+ return(.getLowerSol(L2deriv = L2deriv, risk = risk,
+ neighbor = neighbor, Distr = Distr,
+ DistrSymm = DistrSymm,
+ L2derivSymm = L2derivSymm,
+ L2derivDistrSymm = L2derivDistrSymm,
+ z.start = z.start, A.start = A.start,
+ trafo = trafo, maxiter = maxiter, tol = tol,
+ warn = warn, Finfo = Finfo, verbose = verbose)
+ )
+
+ maxit2 <- if(OptOrIter==2) 0 else maxiter
+ erg <- getLagrangeMultByIter(b = b, L2deriv = L2deriv, risk = risk,
+ trafo = trafo, neighbor = neighbor, biastype = biastype,
+ normtype = normtype, Distr = Distr,
+ z.start = z, A.start = A, w.start = w,
+ std = std, z.comp = z.comp,
+ A.comp = A.comp, maxiter = maxit2, tol = tol,
+ onesetLM = onesetLM, verbose = verbose, warnit = (OptOrIter!=2))
+
+ ## read out solution
+ w <- erg$w
+ A <- erg$A
+ a <- erg$a
+ z <- erg$z
+ biastype <- erg$biastype
+ normtype.old <- normtype
+ normtype <- erg$normtype
+ risk <- erg$risk
+ iter.In <- iter.In + erg$iter
+ prec.In <- erg$prec
+ OptIterCall <- erg$call
+ std <- erg$std
+# print(list(z=z,A=A,b=b))
+
+ if (onesetLM&&maxiter>1){
+ if(is(neighbor,"ContNeighborhood"))
+ cent(w) <- as.numeric(z)
+ if(is(neighbor,"TotalVarNeighborhood"))
+ clip(w) <- c(0,b)+a
+ stand(w) <- A
+ weight(w) <- getweight(w, neighbor = neighbor, biastype = biastype,
+ normW = normtype)
+ }
+ else normtype <- normtype.old
+
+ ## check precision and number of iterations in outer b-loop
+ prec.old <- prec
+ prec <- max(abs(b-b.old), max(abs(A-A.old)), max(abs(z-z.old)))
+ if(verbose)
+ cat("current precision in IC algo:\t", prec, "\n")
+ if(prec < tol) break
+ if(abs(prec.old - prec) < 1e-10){
+ cat("algorithm did not converge!\n", "achieved precision:\t", prec, "\n")
+ break
+ }
+ if(iter > maxiter){
+ cat("maximum iterations reached!\n", "achieved precision:\t", prec, "\n")
+ break
+ }
+ }
+ ### issue some diagnostics if wanted
+ if(verbose){
+ cat("Iterations needed: outer (b-loop):",
+ iter," inner (A,a-loop):", iter.In,"\n")
+ cat("Precision achieved: all in all (b+A,a-loop):",
+ prec," inner (A,a-loop):", prec.In,"\n")
+ }
+ ### determine Covariance of pIC
+ Cov <- getInfV(L2deriv = L2deriv, neighbor = neighbor,
+ biastype = biastype, Distr = Distr,
+ V.comp = A.comp, cent = a,
+ stand = A, w = w)
+ if(!is(risk, "asMSE")){
+ Risk <- getAsRisk(risk = risk, L2deriv = L2deriv, neighbor = neighbor,
+ biastype = biastype, clip = b, cent = a, stand = A,
+ trafo = trafo)
+ }else{
+ Risk <- NULL
+ }
+ }
+
+
+ ### add some further informations for the pIC-slots info and risk
+ info <- paste("optimally robust IC for", sQuote(class(risk)[1]))
+
+ trAsCov <- sum(diag(std%*%Cov))
+ Risk <- c(Risk, list(asCov = Cov,
+ asBias = list(value = b, biastype = biastype,
+ normtype = normtype,
+ neighbortype = class(neighbor)),
+ trAsCov = list(value = trAsCov,
+ normtype = normtype),
+ asMSE = list(value = trAsCov + radius^2*b^2,
+ r = radius,
+ at = neighbor)))
+
+ return(list(A = A, a = a, b = b, d = NULL, risk = Risk, info = info, w = w,
+ biastype = biastype, normtype = normtype,
+ call = mc, iter = iter, prec = prec, OIcall = OptIterCall,
+ iter.In = iter.In, prec.In = prec.In))
+ })
+
+
+### helper function to return the upper case solution if r=0
+.getUpperSol <- function(L2deriv, b, radius, risk, neighbor, biastype,
+ normtype, Distr, Finfo, trafo,
+ QuadForm, verbose, warn){
+
+ if(warn) cat("'radius == 0' => (classical) optimal IC\n",
+ "in sense of Cramer-Rao bound is returned\n")
+ res <- getInfRobIC(L2deriv = L2deriv, risk = asCov(), neighbor = neighbor,
+ Distr = Distr, Finfo = Finfo, trafo = trafo,
+ QuadForm = QuadForm, verbose = verbose)
+ res <- c(res, list(biastype = biastype, normtype = normtype))
+ if(!is(risk, "asMSE")){
+ Risk <- getAsRisk(risk = risk, L2deriv = L2deriv, neighbor = neighbor,
+ biastype = biastype, cent = res$a,
+ stand = res$A, trafo = trafo)
+ res$risk <- c(Risk, res$risk)
+ }
+ trAsCov <- sum(diag(QuadForm%*%res$risk$asCov));
+ res$risk$trAsCov <- list(value = trAsCov, normtype = normtype)
+ res$risk$asBias <- list(value = b, biastype = biastype,
+ normtype = normtype,
+ neighbortype = class(neighbor))
+ res$risk$asMSE <- list(value = trAsCov + radius^2*b^2,
+ r = radius,
+ at = neighbor)
+ return(res)
+}
+
+### helper function to return the lower case solution if b-search was not successful
+.getLowerSol <- function(L2deriv, risk, neighbor, Distr, DistrSymm,
+ L2derivSymm, L2derivDistrSymm,
+ z.start, A.start, trafo,
+ maxiter, tol, warn, Finfo, verbose){
+ if(warn) cat("Could not determine optimal clipping bound!\n",
"'radius >= maximum radius' for the given risk?\n",
"=> the minimum asymptotic bias (lower case) solution is returned\n",
"If 'no' => Try again with modified starting values ",
"'z.start' and 'A.start'\n")
- res <- getInfRobIC(L2deriv = L2deriv,
+ res <- getInfRobIC(L2deriv = L2deriv,
risk = asBias(biastype = biastype(risk),
- normtype = normtype(risk)),
+ normtype = normtype(risk)),
neighbor = neighbor, Distr = Distr, DistrSymm = DistrSymm,
L2derivSymm = L2derivSymm, L2derivDistrSymm = L2derivDistrSymm,
- z.start = z.start, A.start = A.start, trafo = trafo,
+ z.start = z.start, A.start = A.start, trafo = trafo,
maxiter = maxiter, tol = tol, warn = warn, Finfo = Finfo,
verbose = verbose)
normtype(risk) <- res$normtype
if(!is(risk, "asMSE")){
- Risk <- getAsRisk(risk = risk, L2deriv = L2deriv, neighbor = neighbor,
- biastype = biastype, clip = NULL,
+ Risk <- getAsRisk(risk = risk, L2deriv = L2deriv, neighbor = neighbor,
+ biastype = biastype, clip = NULL,
cent = res$a, stand = res$A, trafo = trafo)
res$risk <- c(Risk, res$risk)
}
return(res)
- }
- if(is(neighbor,"ContNeighborhood")){
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/robast -r 337
More information about the Robast-commits
mailing list