[Robast-commits] r409 - in branches/robast-0.8/pkg/ROptEst: . R inst/scripts man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Aug 25 05:02:37 CEST 2010
Author: ruckdeschel
Date: 2010-08-25 05:02:36 +0200 (Wed, 25 Aug 2010)
New Revision: 409
Added:
branches/robast-0.8/pkg/ROptEst/R/AllShow.R
branches/robast-0.8/pkg/ROptEst/R/asAnscombe.R
branches/robast-0.8/pkg/ROptEst/R/getInfRobIC_asAnscombe.R
Modified:
branches/robast-0.8/pkg/ROptEst/NAMESPACE
branches/robast-0.8/pkg/ROptEst/R/AllClass.R
branches/robast-0.8/pkg/ROptEst/R/AllGeneric.R
branches/robast-0.8/pkg/ROptEst/R/getAsRisk.R
branches/robast-0.8/pkg/ROptEst/R/getInfRobIC_asGRisk.R
branches/robast-0.8/pkg/ROptEst/R/leastFavorableRadius.R
branches/robast-0.8/pkg/ROptEst/R/radiusMinimaxIC.R
branches/robast-0.8/pkg/ROptEst/inst/scripts/BinomialModel.R
branches/robast-0.8/pkg/ROptEst/inst/scripts/GammaModel.R
branches/robast-0.8/pkg/ROptEst/inst/scripts/NbinomModel.R
branches/robast-0.8/pkg/ROptEst/inst/scripts/NormalLocationScaleModel.R
branches/robast-0.8/pkg/ROptEst/inst/scripts/NormalScaleModel.R
branches/robast-0.8/pkg/ROptEst/inst/scripts/PoissonModel.R
branches/robast-0.8/pkg/ROptEst/man/getAsRisk.Rd
branches/robast-0.8/pkg/ROptEst/man/getInfRobIC.Rd
branches/robast-0.8/pkg/ROptEst/man/internals.Rd
Log:
+ added new Risk class asAnscombe for determining OBRE to given ARE in the ideal model (see examples in inst folder (Binomial, Gamma, Nbinom, NormalLocationScale, NormalScale, Poisson);
+ getAsRisk now has (almost) the same signature for different methods (so that in can easily be called when the actual risk is not clear)
Modified: branches/robast-0.8/pkg/ROptEst/NAMESPACE
===================================================================
--- branches/robast-0.8/pkg/ROptEst/NAMESPACE 2010-06-10 22:04:52 UTC (rev 408)
+++ branches/robast-0.8/pkg/ROptEst/NAMESPACE 2010-08-25 03:02:36 UTC (rev 409)
@@ -4,6 +4,7 @@
import("distrMod")
import("RobAStBase")
+exportClasses("asAnscombe")
exportMethods("optIC",
"getInfRobIC",
"getFixRobIC",
@@ -25,6 +26,6 @@
"getL1normL2deriv",
"getModifyIC",
"cniperCont", "cniperPoint", "cniperPointPlot")
-exportMethods("updateNorm", "scaleUpdateIC")
-export("getL2normL2deriv")
+exportMethods("updateNorm", "scaleUpdateIC", "eff")
+export("getL2normL2deriv","asAnscombe")
export("roptest","getLagrangeMultByOptim","getLagrangeMultByIter")
Modified: branches/robast-0.8/pkg/ROptEst/R/AllClass.R
===================================================================
--- branches/robast-0.8/pkg/ROptEst/R/AllClass.R 2010-06-10 22:04:52 UTC (rev 408)
+++ branches/robast-0.8/pkg/ROptEst/R/AllClass.R 2010-08-25 03:02:36 UTC (rev 409)
@@ -8,3 +8,14 @@
}
+## asymptotic Anscombe risk
+setClass("asAnscombe", representation(eff = "numeric"),
+ prototype = prototype(eff = .95,
+ type = "optimal bias robust IC for given ARE in the ideal model"),
+ contains = "asRiskwithBias",
+ validity = function(object){
+ if(any(object at eff <= 0|object at eff > 1))
+ stop("'eff' has to be in (0,1]")
+ else TRUE
+ })
+
Modified: branches/robast-0.8/pkg/ROptEst/R/AllGeneric.R
===================================================================
--- branches/robast-0.8/pkg/ROptEst/R/AllGeneric.R 2010-06-10 22:04:52 UTC (rev 408)
+++ branches/robast-0.8/pkg/ROptEst/R/AllGeneric.R 2010-08-25 03:02:36 UTC (rev 409)
@@ -86,3 +86,6 @@
if(!isGeneric("cniperPointPlot")){
setGeneric("cniperPointPlot", function(L2Fam, neighbor, risk, ...) standardGeneric("cniperPointPlot"))
}
+if(!isGeneric("eff")){
+ setGeneric("eff", function(object) standardGeneric("eff"))
+}
Added: branches/robast-0.8/pkg/ROptEst/R/AllShow.R
===================================================================
--- branches/robast-0.8/pkg/ROptEst/R/AllShow.R (rev 0)
+++ branches/robast-0.8/pkg/ROptEst/R/AllShow.R 2010-08-25 03:02:36 UTC (rev 409)
@@ -0,0 +1,6 @@
+setMethod("show", "asAnscombe",
+ function(object){
+ cat(paste("An object of class", dQuote(class(object)), "\n"))
+ cat("risk type:\t", object at type, "\n")
+ cat("ARE in the ideal model:\t", object at eff, "\n")
+ })
Added: branches/robast-0.8/pkg/ROptEst/R/asAnscombe.R
===================================================================
--- branches/robast-0.8/pkg/ROptEst/R/asAnscombe.R (rev 0)
+++ branches/robast-0.8/pkg/ROptEst/R/asAnscombe.R 2010-08-25 03:02:36 UTC (rev 409)
@@ -0,0 +1,8 @@
+## generating function
+asAnscombe <- function(eff = .95, biastype = symmetricBias(),
+ normtype = NormType()){
+ new("asAnscombe", eff = eff, biastype = biastype,
+ normtype = normtype) }
+
+## access method
+setMethod("eff", "asAnscombe", function(object) object at eff)
Modified: branches/robast-0.8/pkg/ROptEst/R/getAsRisk.R
===================================================================
--- branches/robast-0.8/pkg/ROptEst/R/getAsRisk.R 2010-06-10 22:04:52 UTC (rev 408)
+++ branches/robast-0.8/pkg/ROptEst/R/getAsRisk.R 2010-08-25 03:02:36 UTC (rev 409)
@@ -5,7 +5,7 @@
L2deriv = "UnivariateDistribution",
neighbor = "Neighborhood",
biastype = "ANY"),
- function(risk, L2deriv, neighbor, biastype, clip = NULL, cent = NULL, stand, trafo){
+ function(risk, L2deriv, neighbor, biastype, normtype = NULL, clip = NULL, cent = NULL, stand, trafo, ...){
if(!is.finite(neighbor at radius))
mse <- Inf
else
@@ -17,7 +17,7 @@
L2deriv = "EuclRandVariable",
neighbor = "Neighborhood",
biastype = "ANY"),
- function(risk, L2deriv, neighbor, biastype, clip = NULL, cent = NULL, stand, trafo){
+ function(risk, L2deriv, neighbor, biastype, normtype = NULL, clip = NULL, cent = NULL, stand, trafo, ...){
if(!is.finite(neighbor at radius))
mse <- Inf
else{
@@ -37,7 +37,7 @@
L2deriv = "UnivariateDistribution",
neighbor = "ContNeighborhood",
biastype = "ANY"),
- function(risk, L2deriv, neighbor, biastype, trafo){
+ function(risk, L2deriv, neighbor, biastype, normtype = NULL, clip = NULL, cent = NULL, stand = NULL, trafo, ...){
z <- q(L2deriv)(0.5)
bias <- abs(as.vector(trafo))/E(L2deriv, function(x, z){abs(x - z)},
useApply = FALSE, z = z)
@@ -48,7 +48,7 @@
L2deriv = "UnivariateDistribution",
neighbor = "TotalVarNeighborhood",
biastype = "ANY"),
- function(risk, L2deriv, neighbor, biastype, trafo){
+ function(risk, L2deriv, neighbor, biastype, normtype = NULL, clip = NULL, cent = NULL, stand = NULL, trafo, ...){
bias <- abs(as.vector(trafo))/(-m1df(L2deriv, 0))
return(list(asBias = bias))
@@ -57,9 +57,10 @@
L2deriv = "RealRandVariable",
neighbor = "ContNeighborhood",
biastype = "ANY"),
- function(risk, L2deriv, neighbor, biastype, Distr, DistrSymm, L2derivSymm,
+ function(risk, L2deriv, neighbor, biastype, normtype = NULL,
+ clip = NULL, cent = NULL, stand = NULL, Distr, DistrSymm, L2derivSymm,
L2derivDistrSymm, Finfo, trafo, z.start, A.start, maxiter, tol, warn,
- verbose = NULL){
+ verbose = NULL, ...){
if(missing(verbose)|| is.null(verbose))
verbose <- getRobAStBaseOption("all.verbose")
@@ -95,9 +96,10 @@
L2deriv = "RealRandVariable",
neighbor = "TotalVarNeighborhood",
biastype = "ANY"),
- function(risk, L2deriv, neighbor, biastype, Distr, DistrSymm, L2derivSymm,
+ function(risk, L2deriv, neighbor, biastype, normtype = NULL,
+ clip = NULL, cent = NULL, stand = NULL, Distr, DistrSymm, L2derivSymm,
L2derivDistrSymm, Finfo, trafo, z.start, A.start, maxiter, tol, warn,
- verbose = NULL){
+ verbose = NULL, ...){
if(missing(verbose)|| is.null(verbose))
verbose <- getRobAStBaseOption("all.verbose")
@@ -124,44 +126,49 @@
L2deriv = "UnivariateDistribution",
neighbor = "ContNeighborhood",
biastype = "ANY"),
- function(risk, L2deriv, neighbor, biastype, clip, cent, stand){
+ function(risk, L2deriv, neighbor, biastype, normtype = NULL,clip, cent, stand, trafo = NULL, ...){
# c0 <- clip/abs(as.vector(stand))
# D1 <- L2deriv - cent/as.vector(stand)
# Cov <- (clip^2*(p(D1)(-c0) + 1 - p(D1)(c0))
# + as.vector(stand)^2*(m2df(D1, c0) - m2df(D1, -c0)))
- return(list(asCov =
- getInfV(L2deriv, neighbor, biastype(risk), clip/abs(as.vector(stand)),
+ Cov <- getInfV(L2deriv, neighbor, biastype(risk), clip/abs(as.vector(stand)),
cent/abs(as.vector(stand)), abs(as.vector(stand)))
- ))
+
+ if(!is.null(trafo)) Cov <- trafo%*%Cov%*%t(trafo)
+ return(list(asCov = Cov ))
})
setMethod("getAsRisk", signature(risk = "asCov",
L2deriv = "UnivariateDistribution",
neighbor = "TotalVarNeighborhood",
biastype = "ANY"),
- function(risk, L2deriv, neighbor, biastype, clip, cent, stand){
+ function(risk, L2deriv, neighbor, biastype, normtype = NULL, clip, cent, stand, trafo = NULL, ...){
+ if(missing(biastype)||is.null(biastype)) biastype <- biastype(risk)
# g0 <- cent/abs(as.vector(stand))
# c0 <- clip/abs(as.vector(stand))
# Cov <- (abs(as.vector(stand))^2*(g0^2*p(L2deriv)(g0)
# + (g0+c0)^2*(1 - p(L2deriv)(g0+c0))
# + m2df(L2deriv, g0+c0) - m2df(L2deriv, g0)))
# return(list(asCov = Cov))
- return(list(asCov =
- getInfV(L2deriv, neighbor, biastype(risk), clip/abs(as.vector(stand)),
+ Cov <- getInfV(L2deriv, neighbor, biastype, clip/abs(as.vector(stand)),
cent/abs(as.vector(stand)), abs(as.vector(stand)))
- ))
+ if(!is.null(trafo)) Cov <- trafo%*%Cov%*%t(trafo)
+ return(list(asCov = Cov))
})
setMethod("getAsRisk", signature(risk = "asCov",
L2deriv = "RealRandVariable",
neighbor = "ContNeighborhood",
biastype = "ANY"),
- function(risk, L2deriv, neighbor, biastype, Distr, cent,
- stand, V.comp = matrix(TRUE, ncol = nrow(stand), nrow = nrow(stand)),
- w){
- return(getInfV(L2deriv = L2deriv, neighbor = neighbor,
- biastype = biastype(risk), Distr = Distr,
+ function(risk, L2deriv, neighbor, biastype, normtype = NULL, clip=NULL, cent,
+ stand, Distr, trafo = NULL, V.comp = matrix(TRUE, ncol = nrow(stand), nrow = nrow(stand)),
+ w, ...){
+ if(missing(biastype)||is.null(biastype)) biastype <- biastype(risk)
+ Cov <- getInfV(L2deriv = L2deriv, neighbor = neighbor,
+ biastype = biastype, Distr = Distr,
V.comp = V.comp, cent = cent,
- stand = stand, w = w))
+ stand = stand, w = w)
+ if(!is.null(trafo)) Cov <- trafo%*%Cov%*%t(trafo)
+ return(list(asCov = Cov))
})
# Y <- as(stand %*% L2deriv - cent, "EuclRandVariable")
# absY <- sqrt(Y %*% Y)
@@ -181,6 +188,7 @@
+
###############################################################################
## trace of asymptotic covariance
###############################################################################
@@ -188,35 +196,76 @@
L2deriv = "UnivariateDistribution",
neighbor = "UncondNeighborhood",
biastype = "ANY"),
- function(risk, L2deriv, neighbor, biastype, clip, cent, stand){
+ function(risk, L2deriv, neighbor, biastype, normtype = NULL, clip, cent, stand, trafo=NULL, ...){
+ if(missing(biastype)||is.null(biastype)) biastype <- biastype(risk)
+ if(missing(normtype)||is.null(normtype)) normtype <- normtype(risk)
Cov <- getAsRisk(risk = asCov(), L2deriv = L2deriv, neighbor = neighbor,
- biastype = biastype(risk), clip = clip, cent = cent, stand = stand)$asCov
+ biastype = biastype, clip = clip, cent = cent, stand = stand)$asCov
+ std <- if(is(normtype,"QFNorm"))
+ QuadForm(normtype) else diag(nrow(as.matrix(Cov)))
- return(list(trAsCov = as.vector(Cov)))
+ return(list(trAsCov = sum(diag(std%*%as.matrix(Cov)))))
})
setMethod("getAsRisk", signature(risk = "trAsCov",
L2deriv = "RealRandVariable",
neighbor = "ContNeighborhood",
biastype = "ANY"),
- function(risk, L2deriv, neighbor, biastype, Distr, clip, cent, stand, normtype){
+ function(risk, L2deriv, neighbor, biastype, normtype, clip, cent, stand, Distr, trafo = NULL,
+ V.comp = matrix(TRUE, ncol = nrow(stand), nrow = nrow(stand)), w,...){
+ if(missing(biastype)||is.null(biastype)) biastype <- biastype(risk)
+ if(missing(normtype)||is.null(normtype)) normtype <- normtype(risk)
Cov <- getAsRisk(risk = asCov(), L2deriv = L2deriv, neighbor = neighbor,
- biastype = biastype(risk), Distr = Distr, clip = clip,
- cent = cent, stand = stand)$asCov
+ biastype = biastype, Distr = Distr, clip = clip,
+ cent = cent, stand = stand, trafo = trafo,
+ V.comp = V.comp, w = w)$asCov
p <- nrow(stand)
std <- if(is(normtype,"QFNorm")) QuadForm(normtype) else diag(p)
+ print(std)
return(list(trAsCov = sum(diag(std%*%Cov))))
})
###############################################################################
+## Anscombe criterion
+###############################################################################
+setMethod("getAsRisk", signature(risk = "asAnscombe",
+ L2deriv = "UnivariateDistribution",
+ neighbor = "UncondNeighborhood",
+ biastype = "ANY"),
+ function(risk, L2deriv, neighbor, biastype, normtype = NULL, clip, cent, stand, trafo = NULL, FI, ... ){
+ if(missing(biastype)||is.null(biastype)) biastype <- biastype(risk)
+ if(missing(normtype)||is.null(normtype)) normtype <- normtype(risk)
+ trAsCov.0 <- getAsRisk(risk = trAsCov(), L2deriv = L2deriv, neighbor = neighbor,
+ biastype = biastype, normtype = normtype,
+ clip = clip, cent = cent, stand = stand)$trAsCov
+ return(list(asAnscombe = FI/trAsCov.0))
+ })
+setMethod("getAsRisk", signature(risk = "asAnscombe",
+ L2deriv = "RealRandVariable",
+ neighbor = "ContNeighborhood",
+ biastype = "ANY"),
+ function(risk, L2deriv, neighbor, biastype, normtype, clip, cent, stand, Distr, trafo = NULL,
+ V.comp = matrix(TRUE, ncol = nrow(stand), nrow = nrow(stand)), FI,
+ w, ...){
+ if(missing(biastype)||is.null(biastype)) biastype <- biastype(risk)
+ if(missing(normtype)||is.null(normtype)) normtype <- normtype(risk)
+ trAsCov.0 <- getAsRisk(risk = trAsCov(), L2deriv = L2deriv, neighbor = neighbor,
+ biastype = biastype, normtype = normtype,
+ Distr = Distr, clip = clip,
+ cent = cent, stand = stand, V.comp = V.comp,
+ w = w)$trAsCov
+ return(list(asAnscombe = FI/trAsCov.0))
+ })
+
+###############################################################################
## asymptotic under-/overshoot risk
###############################################################################
setMethod("getAsRisk", signature(risk = "asUnOvShoot",
L2deriv = "UnivariateDistribution",
neighbor = "UncondNeighborhood",
biastype = "ANY"),
- function(risk, L2deriv, neighbor, biastype, clip, cent, stand, trafo){
+ function(risk, L2deriv, neighbor, biastype, normtype = NULL,clip, cent, stand, trafo, ...){
if(identical(all.equal(neighbor at radius, 0), TRUE))
return(list(asUnOvShoot = pnorm(-risk at width/sqrt(as.vector(stand)))))
@@ -236,7 +285,8 @@
L2deriv = "UnivariateDistribution",
neighbor = "ContNeighborhood",
biastype = "onesidedBias"),
- function(risk, L2deriv, neighbor, biastype, trafo){
+ function(risk, L2deriv, neighbor, biastype, normtype = NULL,
+ clip = NULL, cent = NULL, stand = NULL, trafo, ...){
D1 <- L2deriv
if(!is(D1, "DiscreteDistribution"))
@@ -265,7 +315,8 @@
L2deriv = "UnivariateDistribution",
neighbor = "ContNeighborhood",
biastype = "asymmetricBias"),
- function(risk, L2deriv, neighbor, biastype, trafo){
+ function(risk, L2deriv, neighbor, biastype, normtype = NULL,
+ clip = NULL, cent = NULL, stand = NULL, trafo, ...){
nu1 <- nu(biastype)[1]
nu2 <- nu(biastype)[2]
num <- nu2/(nu1+nu2)
@@ -284,8 +335,8 @@
L2deriv = "UnivariateDistribution",
neighbor = "Neighborhood",
biastype = "onesidedBias"),
- function(risk, L2deriv, neighbor, biastype,
- clip, cent, stand, trafo){
+ function(risk, L2deriv, neighbor, biastype, normtype = NULL,
+ clip, cent, stand, trafo, ...){
A <- as.vector(stand)*as.vector(trafo)
r <- neighbor at radius
b <- clip*A
Added: branches/robast-0.8/pkg/ROptEst/R/getInfRobIC_asAnscombe.R
===================================================================
--- branches/robast-0.8/pkg/ROptEst/R/getInfRobIC_asAnscombe.R (rev 0)
+++ branches/robast-0.8/pkg/ROptEst/R/getInfRobIC_asAnscombe.R 2010-08-25 03:02:36 UTC (rev 409)
@@ -0,0 +1,231 @@
+###############################################################################
+## IC algorithm for asymptotic Anscombe risk
+###############################################################################
+setMethod("getInfRobIC", signature(L2deriv = "UnivariateDistribution",
+ risk = "asAnscombe",
+ neighbor = "UncondNeighborhood"),
+ function(L2deriv, risk, neighbor, symm, Finfo, trafo,
+ upper = NULL, lower = NULL, maxiter, tol, warn, noLow = FALSE,
+ verbose = NULL, checkBounds = TRUE){
+
+ if(missing(verbose)|| is.null(verbose))
+ verbose <- getRobAStBaseOption("all.verbose")
+
+ eff <- eff(risk)
+
+ i <- 0
+
+ maxi <- min(5,maxiter%/%4)
+ toli <- min(tol*100,1e-3)
+
+ FI0 <- trafo%*%matrix(1/Finfo)%*%t(trafo)
+ normtype <- normtype(risk)
+ std <- if(is(normtype(risk),"QFNorm"))
+ QuadForm(normtype(risk)) else diag(nrow(trafo))
+ FI <- sum(diag(std%*%FI0))
+
+
+ if(is.null(upper))
+ upper <- q(L2deriv)(eff^.5)*3
+ e.up <- 0
+ while(e.up < eff){
+ risk.b <- asHampel(bound = upper, biastype = biastype(risk),
+ normtype = normtype(risk))
+ upBerg <- getInfRobIC(L2deriv, risk.b, neighbor, symm, Finfo, trafo,
+ upper = upper, lower = lower, maxiter = maxi,
+ tol = toli, warn, noLow = noLow,
+ verbose = FALSE, checkBounds = FALSE)
+ trV <- upBerg$risk$trAsCov$value
+ e.up <- FI/trV
+ upper <- upper * 3
+ }
+ upper <- upper / 3
+
+ lowBerg <- minmaxBias(L2deriv = L2deriv, neighbor = neighbor,
+ biastype = biastype(risk), symm = symm,
+ trafo = trafo, maxiter = maxi,
+ tol = toli, warn = warn, Finfo = Finfo)
+ V <- lowBerg$risk$asCov
+ trV <- sum(diag(std%*%V))
+
+
+ if(1/(Vbmin*Finfo)>eff){
+ lowBerg$eff <- FI/trV
+ return(lowBerg)
+ }
+
+ it.erg <- 0
+ erg <- 0
+ if(is.null(lower) || lower < lowBerg$risk$asBias$value)
+ { lower <- lowBerg$risk$asBias$value
+ f.low <- 1/(Vbmin*Finfo)-eff
+ } else f.low <- NULL
+
+ funb <- function(b0){
+ risk.b <- asHampel(bound = b0, biastype = biastype(risk),
+ normtype = normtype(risk))
+ it.erg <<- it.erg + 1
+ maxi <- min(5,maxiter%/%4^(1/it.erg))
+ toli <- min(tol*100^(1/it.erg),1e-3)
+ erg <<- getInfRobIC(L2deriv, risk.b, neighbor, symm, Finfo, trafo,
+ upper = upper, lower = lower, maxiter = maxi, tol = toli, warn, noLow = noLow,
+ verbose = verbose, checkBounds = checkBounds)
+ trV <- erg$risk$trAsCov$value
+ if(verbose) cat("Outer iteration:", it.erg," b_0=", round(b0,3),
+ " eff=", round(1/(v0*Finfo),3), "\n")
+ return(1/(v0*Finfo)-eff)
+ }
+
+ if(is.null(f.low)) f.low <- fun(lower)
+
+ print(c(lower,upper, f.lower=f.low, f.upper=e.up-eff))
+
+ b <- uniroot(funb, interval=c(lower,upper), f.lower=f.low,
+ f.upper=e.up-eff,tol=tol,maxiter=maxiter)
+ erg$info <- c(erg$info,
+ paste("optimally bias-robust IC for ARE", eff, " in the ideal model"))
+
+ erg$risk$eff <- b$f.root+eff
+ return(erg)
+ })
+
+###################################################################################
+# multivariate solution Anscombe --- new 24-08-10
+###################################################################################
+
+setMethod("getInfRobIC", signature(L2deriv = "RealRandVariable",
+ risk = "asAnscombe",
+ 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 = NULL, checkBounds = TRUE, ...){
+
+ if(missing(verbose)|| is.null(verbose))
+ verbose <- getRobAStBaseOption("all.verbose")
+
+ mc <- match.call()
+
+ eff <- eff(risk)
+
+ ## some abbreviations / checks
+ biastype <- biastype(risk)
+ normtype <- normtype(risk)
+
+ p <- nrow(trafo)
+ k <- ncol(trafo)
+
+ maxi <- min(5,maxiter%/%4)
+ toli <- min(tol*100,1e-3)
+
+ std <- if(is(normtype,"QFNorm")) QuadForm(normtype) else diag(p)
+
+ if(! is(neighbor,"ContNeighborhood") && p>1)
+ stop("Not yet implemented")
+
+ ## non-standard norms
+ FI1 <- trafo%*%solve(Finfo)
+ FI0 <- FI1%*%t(trafo)
+ FI <- solve(FI0)
+ if(is(normtype,"InfoNorm") || is(normtype,"SelfNorm") ){
+ QuadForm(normtype) <- PosSemDefSymmMatrix(FI)
+ normtype(risk) <- normtype
+ }
+ std <- if(is(normtype,"QFNorm"))
+ QuadForm(normtype) else diag(p)
+ print(std)
+ trV.ML <- sum(diag(std%*%FI0))
+
+ if(is.null(upper))
+ upper <- sqrt(eff*max(diag(std%*%FI0)))*3
+
+
+ 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 = maxi,
+ tol = toli,
+ warn = FALSE, Finfo = Finfo,
+ QuadForm = std, verbose = verbose)
+
+ if(is.null(lower)||(lower< lowBerg$b))
+ {lower <- lowBerg$b
+ print(lowBerg$risk$asAnscombe)
+ f.low <- lowBerg$risk$asAnscombe - eff
+ } else {
+ risk.b <- asHampel(bound = lower, biastype = biastype,
+ normtype = normtype)
+ lowBerg <- getInfRobIC(L2deriv, risk.b, neighbor,
+ Distr, DistrSymm, L2derivSymm,
+ L2derivDistrSymm, Finfo, trafo, onesetLM = onesetLM,
+ z.start, A.start, upper = upper, lower = lower,
+ OptOrIter = OptOrIter, maxiter=maxi,
+ tol=toli, warn,
+ verbose = FALSE, checkBounds = FALSE, ...)
+ trV <- lowBerg$risk$trAsCov$value
+ f.low <- trV.ML/trV -eff
+ }
+
+ if(f.low > 0){
+ lowBerg$call <- mc
+ lowBerg$eff <- f.low + eff
+ return(lowBerg)
+ }
+
+ e.up <- 0
+ if(lower>=upper) upper <- lower*3
+ while(e.up < eff){
+ risk.b <- asHampel(bound = upper, biastype = biastype,
+ normtype = normtype)
+ upBerg <- getInfRobIC(L2deriv, risk.b, neighbor,
+ Distr, DistrSymm, L2derivSymm,
+ L2derivDistrSymm, Finfo, trafo, onesetLM = onesetLM,
+ z.start, A.start, upper = upper, lower = lower,
+ OptOrIter = OptOrIter, maxiter=maxi,
+ tol=toli, warn,
+ verbose = FALSE, checkBounds = FALSE, ...)
+ trV <- upBerg$risk$trAsCov$value
+ e.up <- trV.ML/trV
+ upper <- upper * 3
+ }
+ upper <- upper / 3
+
+
+ erg <- 0
+ it.erg <- 0
+ funb <- function(b0){
+ risk.b <- asHampel(bound = b0, biastype = biastype(risk),
+ normtype = normtype(risk))
+ it.erg <<- it.erg + 1
+ maxi <- min(5,maxiter%/%4^(1/it.erg))
+ toli <- min(tol*100^(1/it.erg),1e-3)
+ chkbd <- if(it.erg<25) FALSE else checkBounds
+ verbL <- if(it.erg<25) FALSE else verbose
+
+ erg <<- getInfRobIC(L2deriv, risk.b, neighbor,
+ Distr, DistrSymm, L2derivSymm,
+ L2derivDistrSymm, Finfo, trafo, onesetLM = onesetLM,
+ z.start, A.start, upper = upper, lower = lower,
+ OptOrIter = OptOrIter, maxiter = maxi, tol = toli , warn,
+ verbose = verbL, checkBounds = chkbd, ...)
+ trV <- erg$risk$trAsCov$value
+ if(verbose) cat("Outer iteration:", it.erg," b_0=", round(b0,3),
+ " eff=", round(trV.ML/trV,3), "\n")
+ return(trV.ML/trV-eff)
+ }
+ print(c(lower,upper, f.lower=f.low, f.upper=e.up-eff))
+ b <- uniroot(funb, interval=c(lower,upper), f.lower=f.low,
+ f.upper=e.up-eff,tol=tol,maxiter=maxiter)
+ erg$info <- c(erg$info,
+ paste("optimally bias-robust IC for ARE", eff, " in the ideal model"))
+
+ erg$risk$eff <- b$f.root+eff
+ erg$call <- mc
+ return(erg)
+ }
+ )
+
Modified: branches/robast-0.8/pkg/ROptEst/R/getInfRobIC_asGRisk.R
===================================================================
--- branches/robast-0.8/pkg/ROptEst/R/getInfRobIC_asGRisk.R 2010-06-10 22:04:52 UTC (rev 408)
+++ branches/robast-0.8/pkg/ROptEst/R/getInfRobIC_asGRisk.R 2010-08-25 03:02:36 UTC (rev 409)
@@ -14,6 +14,19 @@
biastype <- biastype(risk)
radius <- neighbor at radius
+ p <- nrow(trafo)
+ k <- ncol(trafo)
+
+ ## non-standard norms
+ FI <- solve(trafo%*%matrix(1/Finfo,1,1)%*%t(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)
+ FI0 <- sum(diag(std%*%as.matrix(FI)))
+
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")
@@ -23,8 +36,9 @@
res <- c(res, list(biastype = biastype, normtype = NormType()))
if(!is(risk, "asMSE")){
Risk <- getAsRisk(risk = risk, L2deriv = L2deriv, neighbor = neighbor,
- biastype = biastype, clip = res$b, cent = res$a,
- stand = res$A, trafo = trafo)
+ biastype = biastype, normtype = normtype,
+ clip = res$b, cent = res$a,
+ stand = res$A, trafo = trafo, FI = FI0)
res$risk <- c(Risk, res$risk)
}
Cov <- res$risk$asCov
@@ -87,8 +101,9 @@
maxiter = maxiter, tol = tol, warn = warn,
verbose = verbose)
Risk <- getAsRisk(risk = risk, L2deriv = L2deriv, neighbor = neighbor,
- biastype = biastype, clip = res$b, cent = res$a,
- stand = res$A, trafo = trafo)
+ biastype = biastype, normtype = normtype,
+ clip = res$b, cent = res$a,
+ stand = res$A, trafo = trafo, FI = FI0)
res$risk <- c(Risk, res$risk)
return(res)
}
@@ -126,8 +141,9 @@
b <- abs(as.vector(A))*c0
if(!is(risk, "asMSE")){
Risk <- getAsRisk(risk = risk, L2deriv = L2deriv, neighbor = neighbor,
- biastype = biastype, clip = b, cent = a, stand = A,
- trafo = trafo)
+ biastype = biastype, normtype = normtype,
+ clip = b, cent = a, stand = A,
+ trafo = trafo, FI = FI0)
}else{
Risk <- NULL
}
@@ -201,6 +217,7 @@
}
std <- if(is(normtype,"QFNorm"))
QuadForm(normtype) else diag(p)
+ FI0 <- sum(diag(std%*%as.matrix(FI)))
## starting values
if(is.null(z.start)) z.start <- numeric(k)
@@ -256,7 +273,8 @@
L2derivDistrSymm = L2derivDistrSymm,
z.start = z.start, A.start = A.start,
trafo = trafo, maxiter = maxiter, tol = tol,
- warn = FALSE, Finfo = Finfo, verbose = verbose)
+ warn = FALSE, Finfo = Finfo, QuadForm = std,
+ verbose = verbose)
lower <- lowBerg$b}
#if(is.null(upper))
upper <- 5*max(solve(Finfo))
@@ -288,6 +306,11 @@
biastype <- erg$biastype
normtype.old <- erg$normtype.old
normtype <- erg$normtype
+
+ std.n <- if(is(normtype,"QFNorm"))
+ QuadForm(normtype) else diag(p)
+ FI0 <- sum(diag(std.n%*%as.matrix(FI)))
+
risk <- erg$risk
iter <- erg$iter
prec <- prec.In <- iter.In <- NULL
@@ -298,8 +321,9 @@
if(!is(risk, "asMSE")){
Risk <- getAsRisk(risk = risk, L2deriv = L2deriv, neighbor = neighbor,
- biastype = biastype, clip = b, cent = a, stand = A,
- trafo = trafo)
+ biastype = biastype, normtype = normtype,
+ clip = b, cent = a, stand = A,
+ trafo = trafo, FI = FI0)
}else{
Risk <- sum(diag(std%*%Cov)) + radius^2 * b^2
}
@@ -343,7 +367,8 @@
L2derivDistrSymm = L2derivDistrSymm,
z.start = z.start, A.start = A.start,
trafo = trafo, maxiter = maxiter, tol = tol,
- warn = warn, Finfo = Finfo, verbose = verbose)
+ warn = warn, Finfo = Finfo, QuadForm = std,
+ verbose = verbose)
)
maxit2 <- if(OptOrIter==2) 0 else maxiter
@@ -369,6 +394,7 @@
OptIterCall <- erg$call
std <- if(is.null(erg$std)) std else erg$std
# print(list(z=z,A=A,b=b))
+ FI0 <- sum(diag(std%*%as.matrix(FI)))
## check precision and number of iterations in outer b-loop
@@ -418,8 +444,9 @@
if(verbose && withPICcheck) print(list(Cov=Cov,A=A,a=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)
+ biastype = biastype, normtype = normtype,
+ clip = b, cent = a, stand = A,
+ trafo = trafo, FI = FI0)
}else{
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/robast -r 409
More information about the Robast-commits
mailing list