[Robast-commits] r548 - in branches/robast-0.9/pkg/ROptEst: R man

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Wed Jan 23 02:46:38 CET 2013


Author: ruckdeschel
Date: 2013-01-23 02:46:37 +0100 (Wed, 23 Jan 2013)
New Revision: 548

Modified:
   branches/robast-0.9/pkg/ROptEst/R/getIneffDiff.R
   branches/robast-0.9/pkg/ROptEst/R/internal.roptest.R
   branches/robast-0.9/pkg/ROptEst/R/leastFavorableRadius.R
   branches/robast-0.9/pkg/ROptEst/R/radiusMinimaxIC.R
   branches/robast-0.9/pkg/ROptEst/R/roptest.new.R
   branches/robast-0.9/pkg/ROptEst/man/getIneffDiff.Rd
   branches/robast-0.9/pkg/ROptEst/man/robest.Rd
   branches/robast-0.9/pkg/ROptEst/man/roptest.Rd
Log:
ROptEst: 
+ cleared bug from ticket #131
-> fixed a minor error in helper function .check.eps
-> some more try-catches for radiusMinimaxIC (to keep this clean rewritten, with call
manipulations) and leastFavorableRadius (only minimal call manipulations -> got longer)
-> to do so, getIneffDiff now returns a vector (if called as from radiusMinimaxIC, leastFavorableRadius)
   with lower and upper inefficiency; difference is computed in local functions now
+ timings in roptest, robest are now optional; at any rate are attached as attribute


Modified: branches/robast-0.9/pkg/ROptEst/R/getIneffDiff.R
===================================================================
--- branches/robast-0.9/pkg/ROptEst/R/getIneffDiff.R	2013-01-21 18:55:57 UTC (rev 547)
+++ branches/robast-0.9/pkg/ROptEst/R/getIneffDiff.R	2013-01-23 01:46:37 UTC (rev 548)
@@ -9,10 +9,10 @@
              z.start = NULL, A.start = NULL, upper.b = NULL, lower.b = NULL,
              OptOrIter = "iterate", MaxIter, eps, warn,
              loNorm = NULL, upNorm = NULL,
-             verbose = NULL, ...){
-
+             verbose = NULL, ..., withRetIneff = FALSE){
         if(missing(verbose)|| is.null(verbose))
            verbose <- getRobAStBaseOption("all.verbose")
+
         L2derivDim <- numberOfMaps(L2Fam at L2deriv)
         if(L2derivDim == 1){
             ##print(radius)
@@ -35,9 +35,10 @@
                 ineffUp <- res$b^2/upRisk
             else
                 ineffUp <- (as.vector(res$A)*trafo - res$b^2*(radius^2-upRad^2))/upRisk
-            assign("ineff", ineffUp, envir = sys.frame(which = -4))
+            ##assign("ineff", ineffUp, envir = sys.frame(which = -5))
             ##print(c(ineffUp,ineffLo,ineffUp - ineffLo))
-            return(ineffUp - ineffLo)
+             if(withRetIneff) return(c(lo= ineffLo, up=ineffUp))
+             else return(ineffUp - ineffLo)
         }else{
             if(is(L2Fam at distribution, "UnivariateDistribution")){
                 if((length(L2Fam at L2deriv) == 1) & is(L2Fam at L2deriv[[1]], "RealRandVariable")){
@@ -94,12 +95,12 @@
                 }else{
                     ineffLo <- (sum(diag(std%*%res$A%*%t(trafo))) - 
                                 biasLo^2*(radius^2-loRad^2))/loRisk
-                if(upRad == Inf)
-                    ineffUp <- biasUp^2/upRisk
-                else
-                    ineffUp <- (sum(diag(std%*%res$A%*%t(trafo))) - 
-                                biasUp^2*(radius^2-upRad^2))/upRisk}
-                assign("ineff", ineffUp, envir = sys.frame(which = -4))
+                   if(upRad == Inf)
+                      ineffUp <- biasUp^2/upRisk
+                   else
+                      ineffUp <- (sum(diag(std%*%res$A%*%t(trafo))) -
+                                  biasUp^2*(radius^2-upRad^2))/upRisk
+                }
                 if(verbose)
                     cat(paste(rep("-",75), sep = "", collapse = ""),"\n",
                         "current radius:   ", round(radius,4),
@@ -109,7 +110,8 @@
                                          collapse = ""),"\n",sep="")
                         )
 
-                return(ineffUp - ineffLo)
+                if(withRetIneff) return(c(lo= ineffLo, up=ineffUp))
+                else return(ineffUp - ineffLo)
             }else{
                 stop("not yet implemented")
             }

Modified: branches/robast-0.9/pkg/ROptEst/R/internal.roptest.R
===================================================================
--- branches/robast-0.9/pkg/ROptEst/R/internal.roptest.R	2013-01-21 18:55:57 UTC (rev 547)
+++ branches/robast-0.9/pkg/ROptEst/R/internal.roptest.R	2013-01-23 01:46:37 UTC (rev 548)
@@ -104,20 +104,20 @@
 }
 
 .check.eps <- function(...){
-   mc <- match.call(expand.dots=TRUE)
+   mc <- as.list(match.call(expand.dots=TRUE)[-1])
 
    eps <- eps.lower <- eps.upper <- NULL
+#   print(names(mc))
 #   lapply(mc, function(xx) print(str(xx)))
 
    ine <- is.null(mc[["eps"]]) || is.symbol(mc[["eps"]])
    inl <- is.null(mc[["eps.lower"]]) || is.symbol(mc[["eps.lower"]])
    inu <- is.null(mc[["eps.upper"]]) || is.symbol(mc[["eps.upper"]])
-
+#   print(c(ine=ine,inl=inl,inu=inu))
    if(ine && inl && inu){
         eps.lower <- 0
         eps.upper <- 0.5
-    }
-    if(ine){
+    }else if(ine){
         if(!inl && !inu)
             eps.lower <- mc[["eps.lower"]]
             eps.upper <- mc[["eps.upper"]]

Modified: branches/robast-0.9/pkg/ROptEst/R/leastFavorableRadius.R
===================================================================
--- branches/robast-0.9/pkg/ROptEst/R/leastFavorableRadius.R	2013-01-21 18:55:57 UTC (rev 547)
+++ branches/robast-0.9/pkg/ROptEst/R/leastFavorableRadius.R	2013-01-23 01:46:37 UTC (rev 548)
@@ -1,8 +1,8 @@
 ###############################################################################
-## radius minimax optimally robust IC 
+## radius minimax optimally robust IC
 ## for L2ParamFamily and asymptotic risks
 ###############################################################################
-setMethod("leastFavorableRadius", signature(L2Fam = "L2ParamFamily", 
+setMethod("leastFavorableRadius", signature(L2Fam = "L2ParamFamily",
                                             neighbor = "UncondNeighborhood",
                                             risk = "asGRisk"),
     function(L2Fam, neighbor, risk, rho, upRad = 1,
@@ -22,13 +22,13 @@
         trafo <- trafo(L2Fam at param)
         FI0 <- trafo%*%solve(L2Fam at FisherInfo)%*%t(trafo)
         FI <- solve(FI0)
-        if(is(normtype,"InfoNorm") || is(normtype,"SelfNorm") ) 
-           {QuadForm(normtype) <- PosSemDefSymmMatrix(FI); 
+        if(is(normtype,"InfoNorm") || is(normtype,"SelfNorm") )
+           {QuadForm(normtype) <- PosSemDefSymmMatrix(FI);
             normtype(risk) <- normtype}
 
         L2derivDim <- numberOfMaps(L2Fam at L2deriv)
         if(L2derivDim == 1){
-            leastFavFct <- function(r, L2Fam, neighbor, risk, rho, 
+            leastFavFct <- function(r, L2Fam, neighbor, risk, rho,
                                     upper.b, MaxIter, eps, warn){
                 loRad <- r*rho
                 upRad <- r/rho
@@ -37,61 +37,89 @@
                 ow <- options("warn")
                 on.exit(options(ow))
                 options(warn = -1)
+                .getRisk <- function(rad, fac = 1){
+                    neighbor at radius <- rad
+                    res <- getInfRobIC(L2deriv = L2Fam at L2derivDistr[[1]],
+                                neighbor = neighbor,
+                                risk = risk, symm = L2Fam at L2derivDistrSymm[[1]],
+                                Finfo = L2Fam at FisherInfo, upper = upper.b,
+                                trafo = trafo, maxiter = MaxIter*fac, tol = eps,
+                                warn = warn, verbose = verbose)
+                   return(getAsRisk(risk = risk, L2deriv = L2Fam at L2derivDistr[[1]],
+                                        neighbor = neighbor, biastype = biastype,
+                                        normtype = normtype,
+                                        clip = res$b, cent = res$a,
+                                        stand = res$A, trafo = trafo)[[1]])
+                 }
+
                 if(identical(all.equal(loRad, 0), TRUE)){
                     loRad <- 0
                     loRisk <- 1/as.vector(L2Fam at FisherInfo)
                 }else{
-                    neighbor at radius <- loRad
-                    resLo <- getInfRobIC(L2deriv = L2Fam at L2derivDistr[[1]], neighbor = neighbor, 
-                                risk = risk, symm = L2Fam at L2derivDistrSymm[[1]],
-                                Finfo = L2Fam at FisherInfo, upper = upper.b,
-                                trafo = trafo, maxiter = MaxIter, tol = eps, 
-                                warn = warn, verbose = verbose)
-                    loRisk <- getAsRisk(risk = risk, L2deriv = L2Fam at L2derivDistr[[1]], 
-                                        neighbor = neighbor, biastype = biastype, 
-                                        normtype = normtype,
-                                        clip = resLo$b, cent = resLo$a, 
-                                        stand = resLo$A, trafo = trafo)[[1]]
+                    loRisk <- .getRisk(loRad,6)
                 }
 
                 if(upRad == Inf){
-                    bmin <- getAsRisk(risk = asBias(biastype = biastype), 
-                                L2deriv = L2Fam at L2derivDistr[[1]], 
-                                neighbor = neighbor, biastype = biastype, 
+                    bmin <- getAsRisk(risk = asBias(biastype = biastype),
+                                L2deriv = L2Fam at L2derivDistr[[1]],
+                                neighbor = neighbor, biastype = biastype,
                                 normtype = normtype,
                                 trafo = trafo, symm = L2Fam at L2derivSymm[[1]])
                     upRisk <- bmin^2
                 }else{
-                    neighbor at radius <- upRad
-                    resUp <- getInfRobIC(L2deriv = L2Fam at L2derivDistr[[1]], neighbor = neighbor, 
-                                risk = risk, symm = L2Fam at L2derivDistrSymm[[1]],
-                                Finfo = L2Fam at FisherInfo, upper = upper.b,
-                                trafo = trafo, maxiter = MaxIter, tol = eps, 
-                                warn = warn, verbose = verbose)
-                    upRisk <- getAsRisk(risk = risk, L2deriv = L2Fam at L2derivDistr[[1]], 
-                                        neighbor = neighbor, biastype = biastype, 
-                                        normtype = normtype,
-                                        clip = resUp$b, cent = resUp$a, 
-                                        stand = resUp$A, trafo = trafo)[[1]]
+                    upRisk <- .getRisk(upRad)
                 }
                 loNorm<- upNorm <- NormType()
-                leastFavR <- uniroot(getIneffDiff, lower = lower, upper = upper, 
-                                tol = .Machine$double.eps^0.25, L2Fam = L2Fam, neighbor = neighbor, 
-                                risk = risk, loRad = loRad, upRad = upRad, loRisk = loRisk, 
-                                upRisk = upRisk, upper.b = upper.b, eps = eps, MaxIter = MaxIter, 
-                                warn = warn, 
-                                loNorm = loNorm, upNorm = upNorm)$root
-            options(ow)
+                args.Ie <- list(L2Fam = L2Fam, neighbor = neighbor, risk = risk,
+                                loRad = loRad, upRad = upRad, loRisk = loRisk,
+                                upRisk = upRisk, upper.b = upper.b, eps = eps,
+                                MaxIter = MaxIter, warn = warn,
+                                loNorm = loNorm, upNorm = upNorm,
+                                withRetIneff = TRUE)
+                ineff <- NULL
+                fct.Ie <- function(x){
+                  args.Ie$radius  <- x
+                  res <- do.call(getIneffDiff,args.Ie)
+                  ineff <<- res["up"]
+                  res["up"] - res["lo"]
+                }
+                leastFavR <- try(uniroot(fct.Ie, lower = lower, upper = upper,
+                                 tol = .Machine$double.eps^0.25)$root, silent =TRUE)
+                if(is(leastFavR, "try-error")){
+                   warnRund <- 1; isE <- TRUE
+                   while(warnRund < 7 && isE ){
+                     warnRund <- warnRund + 1
+                     lower <- lower * 2;  upper <- upper / 2
+                     if( warnRund == 4 ) min(upper, 1.5)
+                     if(is.finite(upRad)){
+                        args.Ie$upRad <- upper; args.Ie$upRisk <- .getRisk(upper)
+                     }
+                     if(loRad>0){
+                        args.Ie$loRad <- lower; args.Ie$loRisk <- .getRisk(lower)
+                     }
+                     leastFavR <- try(
+                         uniroot(fct.Ie, lower = lower, upper = upper,
+                         tol = .Machine$double.eps^0.25)$root, silent = TRUE)
+                     isE <- is(leastFavR, "try-error")
+                     if(isE) print(conditionMessage(attr(leastFavR,"condition")))
+                   }
+                   if(isE)
+                      stop("Problem: Zero search in getIneffDiff did not converge.")
+                   else warning(paste("Had to modify radius bounds to [", lower,
+                        upper, "] after", warnRund, "iterations."))
+                }
+
+                options(ow)
                 cat("current radius:\t", r, "\tinefficiency:\t", ineff, "\n")
                 return(ineff)
             }
-            leastFavR <- optimize(leastFavFct, lower = 1e-4, upper = upRad, 
+            leastFavR <- optimize(leastFavFct, lower = 1e-4, upper = upRad,
                             tol = .Machine$double.eps^0.25, maximum = TRUE,
                             L2Fam = L2Fam, neighbor = neighbor, risk = risk,
-                            rho = rho, upper.b = upper, MaxIter = maxiter, 
+                            rho = rho, upper.b = upper, MaxIter = maxiter,
                             eps = tol, warn = warn)
 
-            res <- list(rho = rho, leastFavorableRadius = leastFavR$maximum, 
+            res <- list(rho = rho, leastFavorableRadius = leastFavR$maximum,
                         ineff = leastFavR$objective)
             names(res)[3] <- paste(class(risk)[1], "-inefficiency", sep="")
             return(res)
@@ -116,12 +144,14 @@
                         L2derivDistrSymm <- new("DistrSymmList", L2)
                     }
                 }
-   
-                std <- if(is(normtype,"QFNorm")) 
+
+                std <- if(is(normtype,"QFNorm"))
                        QuadForm(normtype) else diag(nrow(trafo))
-   
-                leastFavFct <- function(r, L2Fam, neighbor, risk, rho, 
+
+
+                leastFavFct <- function(r, L2Fam, neighbor, risk, rho,
                                         z.start, A.start, upper.b, MaxIter, eps, warn){
+
                     loRad <- r*rho
                     upRad <- r/rho
                     lower <- ifelse(identical(all.equal(loRad, 0), TRUE), 1e-4, loRad)
@@ -129,75 +159,109 @@
                     ow <- options("warn")
                     on.exit(options(ow))
                     options(warn = -1)
+
+                    .getRisk <- function(rad, fac = 1){
+                        neighbor at radius <- rad
+                        res <- getInfRobIC(L2deriv = L2deriv, neighbor = neighbor,
+                                        risk = risk, 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, OptOrIter = OptOrIter,
+                                        maxiter = MaxIter, tol = eps, warn = warn,
+                                        verbose = verbose)
+                        risk.0 <- risk
+                        normtype(risk.0) <- .Norm <- res$normtype
+                        .Risk <- getAsRisk(risk = risk.0, L2deriv = L2deriv,
+                                          neighbor = neighbor, biastype = biastype,
+                                          normtype = normtype, clip = res$b,
+                                          cent = res$a, stand = res$A,
+                                          trafo = trafo)[[1]]
+                        return(list(Risk=.Risk,Norm=.Norm))
+                    }
+
                     if(identical(all.equal(loRad, 0), TRUE)){
                         loRad <- 0
                         loRisk <- sum(diag(std%*%FI0))
-                        loNorm <- normtype                    
+                        loNorm <- normtype
                     }else{
-                        neighbor at radius <- loRad
-                        resLo <- getInfRobIC(L2deriv = L2deriv, neighbor = neighbor, risk = risk, 
-                                    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,
-                                    OptOrIter = OptOrIter, maxiter = MaxIter,
-                                    tol = eps, warn = warn, verbose = verbose)
-                        riskLo <- risk
-                        normtype(riskLo) <- resLo$normtype
-                        loRisk <- getAsRisk(risk = riskLo, L2deriv = L2deriv, neighbor = neighbor, 
-                                            biastype = biastype, normtype = normtype,
-                                            clip = resLo$b, cent = resLo$a, 
-                                            stand = resLo$A, trafo = trafo)[[1]]
-                        loNorm <- resLo$normtype                    
+                        rL <- .getRisk(loRad)
+                        loRisk <- rL$Risk
+                        loNorm <- rL$Norm
                     }
 
                     if(upRad == Inf){
-                        biasR <- getAsRisk(risk = asBias(biastype = biastype(risk), 
-                                      normtype = normtype), L2deriv = L2deriv, 
-                                      neighbor = neighbor, biastype = biastype, 
+                        biasR <- getAsRisk(risk = asBias(biastype = biastype(risk),
+                                      normtype = normtype), L2deriv = L2deriv,
+                                      neighbor = neighbor, biastype = biastype,
                                       normtype = normtype,
-                                      Distr = L2Fam at distribution, 
-                                      DistrSymm = L2Fam at distrSymm, 
-                                      L2derivSymm = L2derivSymm, 
-                                      L2derivDistrSymm= L2derivDistrSymm,                                       
+                                      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,
-                                maxiter = maxiter, tol = tol, 
+                                maxiter = maxiter, tol = tol,
                                 warn = warn, verbose = verbose)
                         bmin <- biasR$asBias
                         upRisk <- bmin^2
                         upNorm <- biasR$normtype
                     }else{
-                        neighbor at radius <- upRad
-                        resUp <- getInfRobIC(L2deriv = L2deriv, neighbor = neighbor, risk = risk, 
-                                    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,
-                                    OptOrIter = OptOrIter, maxiter = maxiter,
-                                    tol = tol, warn = warn, verbose = verbose)
-                         riskUp <- risk
-                         normtype(riskUp) <- resUp$normtype
-                         upRisk <- getAsRisk(risk = riskUp, L2deriv = L2deriv, neighbor = neighbor, 
-                                        biastype = biastype, normtype = normtype,
-                                        clip = resUp$b, cent = resUp$a, stand = resUp$A, trafo = trafo)[[1]]
-                         upNorm <- resUp$normtype                    
+                        rL <- .getRisk(upRad)
+                        upRisk <- rL$Risk
+                        upNorm <- rL$Norm
                     }
-                    leastFavR <- uniroot(getIneffDiff, lower = lower, upper = upper, 
-                                    tol = .Machine$double.eps^0.25, L2Fam = L2Fam,
+                    args.Ie <- list(L2Fam = L2Fam,
                                     neighbor = neighbor, z.start = z.start,
                                     A.start = A.start, upper.b = upper.b,
-                                    risk = risk, 
+                                    risk = risk,
                                     loRad = loRad, upRad = upRad,
                                     loRisk = loRisk, upRisk = upRisk,
                                     eps = eps, OptOrIter = OptOrIter,
                                     MaxIter = MaxIter, warn = warn,
-                                    loNorm = loNorm, upNorm = upNorm)$root
+                                    loNorm = loNorm, upNorm = upNorm,
+                                    withRetIneff = TRUE)
+                    ineff <- NULL
+                    fct.Ie <- function(x){
+                      args.Ie$radius  <- x
+                      res <- do.call(getIneffDiff,args.Ie)
+                      ineff <<- res["up"]
+                      res["up"] - res["lo"]
+                    }
+                    leastFavR <- try(uniroot(fct.Ie, lower = lower, upper = upper,
+                                 tol = .Machine$double.eps^0.25)$root, silent =TRUE)
+                    if(is(leastFavR, "try-error")){
+                       warnRund <- 1; isE <- TRUE
+                       while(warnRund < 7 && isE ){
+                         warnRund <- warnRund + 1
+                         lower <- lower * 2;  upper <- upper / 2
+                         if( warnRund == 4 ) min(upper, 1.5)
+                         if(is.finite(upRad)){
+                            args.Ie$upRad <- upper; rL <- .getRisk(upper)
+                            args.Ie$upRisk <- rL$Risk; args.Ie$upNorm <- rL$Norm
+                         }
+                         if(loRad>0){
+                            args.Ie$upRad <- upper; rL <- .getRisk(upper)
+                            args.Ie$upRisk <- rL$Risk; args.Ie$upNorm <- rL$Norm
+                         }
+                         leastFavR <- try(
+                             uniroot(fct.Ie, lower = lower, upper = upper,
+                             tol = .Machine$double.eps^0.25)$root, silent = TRUE)
+                         isE <- is(leastFavR, "try-error")
+                         if(isE) print(conditionMessage(attr(leastFavR,"condition")))
+                       }
+                       if(isE)
+                          stop("Problem: Zero search in getIneffDiff did not converge.")
+                       else warning(paste("Had to modify radius bounds to [", lower,
+                            upper, "] after", warnRund, "iterations."))
+                    }
                     options(ow)
 
                     if(verbose)
                        cat(paste(rep("-",75), sep = "", collapse = ""),"\n")
-                    cat("current radius:   ", round(radius,4),
+                    cat("current radius:   ", round(r,4),
                         "\tinefficiency:   ", round(ineff,4))
                     if(verbose)
                        cat(paste("\n",paste(rep("-",75), sep = "",
@@ -208,13 +272,13 @@
                 }
                 if(is.null(z.start)) z.start <- numeric(L2derivDim)
                 if(is.null(A.start)) A.start <- trafo
-                leastFavR <- optimize(leastFavFct, lower = 1e-4, upper = upRad, 
+                leastFavR <- optimize(leastFavFct, lower = 1e-4, upper = upRad,
                                 tol = .Machine$double.eps^0.25, maximum = TRUE,
                                 L2Fam = L2Fam, neighbor = neighbor, risk = risk,
-                                rho = rho, z.start = z.start, A.start = A.start, 
+                                rho = rho, z.start = z.start, A.start = A.start,
                                 upper.b = upper, MaxIter = maxiter, eps = tol, warn = warn)
 
-                res <- list(rho = rho, leastFavorableRadius = leastFavR$maximum, 
+                res <- list(rho = rho, leastFavorableRadius = leastFavR$maximum,
                             ineff = leastFavR$objective)
                 names(res)[3] <- paste(class(risk)[1], "-inefficiency", sep="")
                 return(res)

Modified: branches/robast-0.9/pkg/ROptEst/R/radiusMinimaxIC.R
===================================================================
--- branches/robast-0.9/pkg/ROptEst/R/radiusMinimaxIC.R	2013-01-21 18:55:57 UTC (rev 547)
+++ branches/robast-0.9/pkg/ROptEst/R/radiusMinimaxIC.R	2013-01-23 01:46:37 UTC (rev 548)
@@ -27,84 +27,70 @@
         if(is(normtype(risk),"SelfNorm")||is(normtype(risk),"InfoNorm"))
            upRad <- min(upRad,10) 
 
+        ineff <- NULL
+        args.IC <- list(L2deriv = NULL, neighbor = NULL,
+                       risk = risk, symm = NULL,
+                       Finfo = L2Fam at FisherInfo, upper = upper, lower = lower,
+                       trafo = trafo, maxiter = maxiter, tol = tol,
+                       warn = warn, verbose = verbose)
+
+        args.R  <- list(risk = risk, L2deriv = NULL,
+                         neighbor = neighbor, biastype = biastype,
+                         normtype = normtype(risk), trafo = trafo)
+
+        args.Ie <- list(radius = NULL, L2Fam = L2Fam, neighbor,
+                        risk = risk, upper.b = upper, lower.b = lower,
+                        loRad = loRad, upRad = upRad, eps = tol,
+                        MaxIter = maxiter, warn = warn,
+                        loNorm = NormType(), upNorm = NormType(),
+                        verbose=verbose, withRetIneff = TRUE)
+        fct.Ie <- function(x){
+                 args.Ie$radius  <- x
+#                 print(with(args.Ie, list(loRisk,upRisk,loRad,upRad)))
+                 res <- do.call(getIneffDiff,args.Ie)
+#                 print(res)
+                 ineff <<- res["up"]
+                 res["up"] - res["lo"]
+        }
+
+        .getRisk <- function(rad, fac = 1){
+               neighbor at radius <- rad
+               args.IC$maxiter <- maxiter * fac
+               args.IC$neighbor <- args.R$neighbor <- neighbor
+               res <- do.call(getInfRobIC, args.IC)
+               args.R$clip <- res$b; args.R$stand <- res$A; args.R$cent <- res$a
+               Norm <- NormType()
+               if(L2derivDim > 1){
+                   risk.0 <- risk; normtype(risk.0) <- Norm <- res$normtype
+                   args.R$risk <- risk.0
+               }
+               return(list(Risk=do.call(getAsRisk, args.R)[[1]], Norm=Norm))
+        }
+
         if(L2derivDim == 1){
             options(warn = -1)
-            upper.b <- upper
-            lower.b <- lower
-            lower <- max(loRad, loRad0)
-            upper <- if(upRad == Inf) max(lower+2, 4) else upRad
-            if(is(neighbor,"TotalVarNeighborhood")) {upper <- upper/2}
-            if(identical(all.equal(loRad, 0), TRUE)){
-                loRad <- 0
-                loRisk <- 1/as.vector(L2Fam at FisherInfo)
-            }else{
-                neighbor at radius <- loRad
-                resLo <- getInfRobIC(L2deriv = L2Fam at L2derivDistr[[1]], neighbor = neighbor, 
-                            risk = risk, symm = L2Fam at L2derivDistrSymm[[1]],
-                            Finfo = L2Fam at FisherInfo, upper = upper.b, lower = lower.b,
-                            trafo = trafo, maxiter = maxiter*6, tol = tol,
-                            warn = warn, verbose = verbose)
-                loRisk <- getAsRisk(risk = risk, L2deriv = L2Fam at L2derivDistr[[1]], 
-                                    neighbor = neighbor, biastype = biastype,
-                                    normtype = normtype(risk),
-                                    clip = resLo$b, cent = resLo$a, 
-                                    stand = resLo$A, trafo = trafo)[[1]]
-            }
+            args.R$L2deriv <- args.IC$L2deriv <- L2Fam at L2derivDistr[[1]]
+            args.IC$symm <-  L2Fam at L2derivDistrSymm[[1]]
 
-            if(upRad == Inf){
-                bmin <- getAsRisk(risk = asBias(biastype = biastype), 
-                                  L2deriv = L2Fam at L2derivDistr[[1]], 
-                                  neighbor = neighbor, biastype = biastype, 
-                                  normtype = normtype(risk),
-                                  trafo = trafo)$asBias
-                upRisk <- bmin^2
-            }else{
-                neighbor at radius <- upRad
-                resUp <- getInfRobIC(L2deriv = L2Fam at L2derivDistr[[1]], neighbor = neighbor, 
-                            risk = risk, symm = L2Fam at L2derivDistrSymm[[1]],
-                            Finfo = L2Fam at FisherInfo, upper = upper.b, lower=lower.b,
-                            trafo = trafo, maxiter = maxiter, tol = tol,
-                            warn = warn, verbose = verbose)
-                upRisk <- getAsRisk(risk = risk, L2deriv = L2Fam at L2derivDistr[[1]], 
-                                    neighbor = neighbor, biastype = biastype, 
-                                    normtype = normtype(risk),
-                                    clip = resUp$b, cent = resUp$a, 
-                                    stand = resUp$A, trafo = trafo)[[1]]
-            }
+            if(is(neighbor,"TotalVarNeighborhood")) upper <- upper/2
 
-            loNorm<- upNorm <- NormType()
-##            print(c(loRad=loRad,loRisk=loRisk,lower=lower,upRad=upRad,upRisk=upRisk,upper=upper))
-            leastFavR <- uniroot(getIneffDiff, lower = lower, upper = upper,
-                            tol = .Machine$double.eps^0.25, L2Fam = L2Fam, neighbor = neighbor, 
-                            upper.b = upper.b, lower.b = lower.b, risk = risk, loRad = loRad, upRad = upRad,
-                            loRisk = loRisk, upRisk = upRisk, eps = tol, 
-                            MaxIter = maxiter, warn = warn, 
-                            loNorm = loNorm, upNorm = upNorm, verbose=verbose)$root
-            neighbor at radius <- leastFavR
-            res <- getInfRobIC(L2deriv = L2Fam at L2derivDistr[[1]], neighbor = neighbor, 
-                        risk = risk, symm = L2Fam at L2derivSymm[[1]],
-                        Finfo = L2Fam at FisherInfo, upper = upper.b, lower = lower.b,
-                        trafo = trafo, maxiter = maxiter, tol = tol, 
-                        warn = warn, verbose = verbose)
-            options(ow)
-            res$info <- c("radiusMinimaxIC", paste("radius minimax IC for radius interval [", 
-                            round(loRad, 3), ", ", round(upRad, 3), "]", sep=""))
-            res$info <- rbind(res$info, c("radiusMinimaxIC", 
-                            paste("least favorable radius: ", round(leastFavR, 3), sep="")))
-            res$info <- rbind(res$info, c("radiusMinimaxIC", 
-                            paste("maximum ", sQuote(class(risk)[1]), "-inefficiency: ",
-                            round(ineff, 3), sep="")))
-            res <- c(res, modifyIC = getModifyIC(L2FamIC = L2Fam, 
-                                                 neighbor = neighbor, 
-                                                 risk = risk))
-            return(generateIC(neighbor, L2Fam, res))
+            args.Ie$loRisk <- if(identical(all.equal(loRad, 0), TRUE))
+                1/as.vector(L2Fam at FisherInfo) else .getRisk(loRad, 6)$Risk
+
+            if(upRad == Inf){
+                args.lR <- args.R
+                args.lR$risk <- asBias(biastype = biastype)
+                args.Ie$upRisk <- (do.call(getAsRisk, args.lR)$asBias)^2
+            }else args.Ie$upRisk <- .getRisk(upRad)$Risk
+#            print(c(rlo=loRad, Rlo=args.Ie$loRisk, rup=upRad,Rup=args.Ie$upRisk))
         }else{
             if(is(L2Fam at distribution, "UnivariateDistribution")){
-                if((length(L2Fam at L2deriv) == 1) & is(L2Fam at L2deriv[[1]], "RealRandVariable")){
+               if((length(L2Fam at L2deriv) == 1) &
+                     is(L2Fam at L2deriv[[1]], "RealRandVariable")){
                     L2deriv <- L2Fam at L2deriv[[1]]
                     L2derivSymm <- L2Fam at L2derivSymm
                     L2derivDistrSymm <- L2Fam at L2derivDistrSymm
-                }else{
+               }else{
                     L2deriv <- diag(dimension(L2Fam at L2deriv)) %*% L2Fam at L2deriv
                     L2deriv <- RealRandVariable(Map = L2deriv at Map, Domain = L2deriv at Domain)
                     nrvalues <- numberOfMaps(L2deriv)
@@ -118,110 +104,105 @@
                         L2derivSymm <- new("FunSymmList", L1)
                         L2derivDistrSymm <- new("DistrSymmList", L2)
                     }
-                }
-                normtype <- normtype(risk)
+               }
+               li.0 <- list(z.start = z.start, A.start = A.start,
+                            OptOrIter = OptOrIter)
+               li.1 <- list(Distr = L2Fam at distribution,
+                            DistrSymm = L2Fam at distrSymm,
+                            L2derivSymm = L2derivSymm,
+                            L2derivDistrSymm = L2derivDistrSymm)
 
-                Finfo <- L2Fam at FisherInfo
+               args.R$L2deriv <- args.IC$L2deriv <- L2deriv
+               args.IC <- c(args.IC, li.1, li.0)
+               args.Ie <- c(args.Ie, li.0)
 
[TRUNCATED]

To get the complete diff run:
    svnlook diff /svnroot/robast -r 548


More information about the Robast-commits mailing list