[Robast-commits] r961 - in branches/robast-1.1/pkg/ROptEstOld: . R

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Tue Jul 17 14:26:22 CEST 2018


Author: ruckdeschel
Date: 2018-07-17 14:26:22 +0200 (Tue, 17 Jul 2018)
New Revision: 961

Modified:
   branches/robast-1.1/pkg/ROptEstOld/DESCRIPTION
   branches/robast-1.1/pkg/ROptEstOld/NAMESPACE
   branches/robast-1.1/pkg/ROptEstOld/R/AllInitialize.R
   branches/robast-1.1/pkg/ROptEstOld/R/ContIC.R
   branches/robast-1.1/pkg/ROptEstOld/R/getAsRisk.R
   branches/robast-1.1/pkg/ROptEstOld/R/getIneffDiff.R
   branches/robast-1.1/pkg/ROptEstOld/R/ksEstimator.R
   branches/robast-1.1/pkg/ROptEstOld/R/leastFavorableRadius.R
   branches/robast-1.1/pkg/ROptEstOld/R/radiusMinimaxIC.R
Log:
[ROptEst] branch 1.1 some NAMESPACE items added, changed shakey call to sys.frame(...) 

Modified: branches/robast-1.1/pkg/ROptEstOld/DESCRIPTION
===================================================================
--- branches/robast-1.1/pkg/ROptEstOld/DESCRIPTION	2018-07-17 11:20:51 UTC (rev 960)
+++ branches/robast-1.1/pkg/ROptEstOld/DESCRIPTION	2018-07-17 12:26:22 UTC (rev 961)
@@ -5,7 +5,7 @@
 Description: Optimally robust estimation using S4 classes and methods. Old version still needed
         for current versions of ROptRegTS and RobRex.
 Depends: R(>= 2.14.0), methods, distr(>= 2.5.2), distrEx(>= 2.2), RandVar(>= 0.9.2), evd
-Authors at R: person("Matthias", "Kohl", role=c("cre", "cph"), email="Matthias.Kohl at stamats.de")
+Authors at R: person("Matthias", "Kohl", role=c("aut", "cre", "cph"), email="Matthias.Kohl at stamats.de")
 ByteCompile: yes
 License: LGPL-3
 URL: http://robast.r-forge.r-project.org/

Modified: branches/robast-1.1/pkg/ROptEstOld/NAMESPACE
===================================================================
--- branches/robast-1.1/pkg/ROptEstOld/NAMESPACE	2018-07-17 11:20:51 UTC (rev 960)
+++ branches/robast-1.1/pkg/ROptEstOld/NAMESPACE	2018-07-17 12:26:22 UTC (rev 961)
@@ -3,6 +3,11 @@
 import("distrEx")
 import("RandVar")
 import("evd")
+importFrom("grDevices", "grey")
+importFrom("graphics", "legend", "lines", "par", "title")
+importFrom("stats", "approxfun", "dbinom", "ecdf", "fft", "ks.test",
+             "optim", "optimize", "pbinom", "pnorm", "ppois", "qpois",
+             "uniroot")
 
 exportClasses("FunctionSymmetry",
               "NonSymmetric",

Modified: branches/robast-1.1/pkg/ROptEstOld/R/AllInitialize.R
===================================================================
--- branches/robast-1.1/pkg/ROptEstOld/R/AllInitialize.R	2018-07-17 11:20:51 UTC (rev 960)
+++ branches/robast-1.1/pkg/ROptEstOld/R/AllInitialize.R	2018-07-17 12:26:22 UTC (rev 961)
@@ -14,6 +14,7 @@
         body(.Object at p) <- substitute({p1 <- pgumbel(q, loc = loc1, scale = scale1, lower.tail = lower.tail) 
                                        return(if(log.p) log(p1) else p1)},
                                      list(loc1 = loc, scale1 = scale))
+        loc1 <- loc; scale1 <- scale
         .Object at q <- function(p, loc = loc1, scale = scale1, lower.tail = TRUE, log.p = FALSE){}
             body(.Object at q) <- substitute({
                         ## P.R.: changed to vectorized form 

Modified: branches/robast-1.1/pkg/ROptEstOld/R/ContIC.R
===================================================================
--- branches/robast-1.1/pkg/ROptEstOld/R/ContIC.R	2018-07-17 11:20:51 UTC (rev 960)
+++ branches/robast-1.1/pkg/ROptEstOld/R/ContIC.R	2018-07-17 12:26:22 UTC (rev 961)
@@ -56,17 +56,17 @@
         Y <- as(A %*% L2Fam at L2deriv - a, "EuclRandVariable")
         if(nrvalues == 1){
             if(!is.null(d)){
-                ICfct[[1]] <- function(x){ 
-                                    ind <- (Y(x) != 0) 
-                                    b*(ind*Y(x)/(ind*absY(x) + (1-ind)) + zi*(1-ind)*d)
-                              }
+                ICfct[[1]] <- function(x){}#
+                                    #ind <- (Y(x) != 0)
+                                   # b*(ind*Y(x)/(ind*absY(x) + (1-ind)) + zi*(1-ind)*d)
+                              #}
                 body(ICfct[[1]]) <- substitute(
                                         { ind <- (Y(x) != 0) 
                                           b*(ind*Y(x)/(ind*absY(x) + (1-ind)) + zi*(1-ind)*d) },
                                         list(Y = Y at Map[[1]], absY = abs(Y)@Map[[1]], b = b, d = d, 
                                              zi = sign(L2Fam at param@trafo)))
             }else{
-                ICfct[[1]] <- function(x){ Y(x)*pmin(1, b/absY(x)) }
+                ICfct[[1]] <- function(x){}# Y(x)*pmin(1, b/absY(x)) }
                 body(ICfct[[1]]) <- substitute({ Y(x)*pmin(1, b/absY(x)) },
                                                  list(Y = Y at Map[[1]], absY = abs(Y)@Map[[1]], b = b))
             }
@@ -74,13 +74,13 @@
             absY <- sqrt(Y %*% Y)
             if(!is.null(d))
                 for(i in 1:nrvalues){
-                    ICfct[[i]] <- function(x){ ind <- (Yi(x) != 0) ; ind*b*Yi(x)/absY(x) + (1-ind)*d }
+                    ICfct[[i]] <- function(x){}# ind <- (Yi(x) != 0) ; ind*b*Yi(x)/absY(x) + (1-ind)*d }
                     body(ICfct[[i]]) <- substitute({ ind <- (Yi(x) != 0) ; ind*b*Yi(x)/absY(x) + (1-ind)*d },
                                                  list(Yi = Y at Map[[i]], absY = absY at Map[[1]], b = b, d = d[i]))
                 }
             else
                 for(i in 1:nrvalues){
-                    ICfct[[i]] <- function(x){ Yi(x)*pmin(1, b/absY(x)) }
+                    ICfct[[i]] <- function(x){}# Yi(x)*pmin(1, b/absY(x)) }
                     body(ICfct[[i]]) <- substitute({ Yi(x)*pmin(1, b/absY(x)) },
                                                  list(Yi = Y at Map[[i]], absY = absY at Map[[1]], b = b))
                 }

Modified: branches/robast-1.1/pkg/ROptEstOld/R/getAsRisk.R
===================================================================
--- branches/robast-1.1/pkg/ROptEstOld/R/getAsRisk.R	2018-07-17 11:20:51 UTC (rev 960)
+++ branches/robast-1.1/pkg/ROptEstOld/R/getAsRisk.R	2018-07-17 12:26:22 UTC (rev 961)
@@ -120,7 +120,7 @@
         nrvalues <- nrow(stand)
         ICfct <- vector(mode = "list", length = nrvalues)
         for(i in 1:nrvalues){
-            ICfct[[i]] <- function(x){ Yi(x)*pmin(1, b/absY(x)) }
+            ICfct[[i]] <- function(x){}# Yi(x)*pmin(1, b/absY(x)) }
             body(ICfct[[i]]) <- substitute({ Yi(x)*pmin(1, b/absY(x)) },
                                     list(Yi = Y at Map[[i]], absY = absY at Map[[1]], b = clip))
         }

Modified: branches/robast-1.1/pkg/ROptEstOld/R/getIneffDiff.R
===================================================================
--- branches/robast-1.1/pkg/ROptEstOld/R/getIneffDiff.R	2018-07-17 11:20:51 UTC (rev 960)
+++ branches/robast-1.1/pkg/ROptEstOld/R/getIneffDiff.R	2018-07-17 12:26:22 UTC (rev 961)
@@ -20,10 +20,11 @@
                 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))
+## changed: shakey...            assign("ineff", ineffUp, envir = sys.frame(which = -4))
+#            return(ineffUp - ineffLo)
+        return(c(ineff=ineffUp,ineffDiff=ineffUp - ineffLo))
+        }else{
 
-            return(ineffUp - ineffLo)
-        }else{
             if(is(L2Fam at distribution, "UnivariateDistribution")){
                 if((length(L2Fam at L2deriv) == 1) & is(L2Fam at L2deriv[[1]], "RealRandVariable")){
                     L2deriv <- L2Fam at L2deriv[[1]]
@@ -57,10 +58,11 @@
                     ineffUp <- res$b^2/upRisk
                 else
                     ineffUp <- (sum(diag(res$A%*%t(trafo))) - res$b^2*(radius^2-upRad^2))/upRisk
-                assign("ineff", ineffUp, envir = sys.frame(which = -4))
+     ## changed: shakey  assign("ineff", ineffUp, envir = sys.frame(which = -4))
                 cat("current radius:\t", radius, "\tMSE-inefficiency difference:\t", ineffUp - ineffLo, "\n")
 
-                return(ineffUp - ineffLo)
+     ##           return(ineffUp - ineffLo)
+                return(c(ineff=ineffUp,ineffDiff=ineffUp - ineffLo))
             }else{
                 stop("not yet implemented")
             }

Modified: branches/robast-1.1/pkg/ROptEstOld/R/ksEstimator.R
===================================================================
--- branches/robast-1.1/pkg/ROptEstOld/R/ksEstimator.R	2018-07-17 11:20:51 UTC (rev 960)
+++ branches/robast-1.1/pkg/ROptEstOld/R/ksEstimator.R	2018-07-17 12:26:22 UTC (rev 961)
@@ -22,23 +22,23 @@
         if(param == "prob"){
             if(max(x) > size(distribution))
                 stop("maximum of 'x' > 'size' of distribution")
-            KSdist <- function(prob, x, size){
+            KSdist.p <- function(prob, x, size){
                 supp <- 0:size
                 edf <- ecdf(x)
                 return(max(abs(edf(supp)-pbinom(supp, size = size, prob = prob))))
             }
-            res <- optimize(f = KSdist, interval = c(0, 1), tol = eps, 
+            res <- optimize(f = KSdist.p, interval = c(0, 1), tol = eps,
                      x = x, size = size(distribution))$minimum
             return(list(size = size(distribution), prob = res))
         }
         if(param == "size"){
-            KSdist <- function(size, x, prob){
+            KSdist.s <- function(size, x, prob){
                 supp <- 0:size
                 edf <- ecdf(x)
                 return(max(abs(edf(supp)-pbinom(supp, size = size, prob = prob))))
             }
             size <- max(x):(max(x) + 100)
-            ind <- which.min(sapply(size, KSdist, x = x, prob = prob(distribution)))
+            ind <- which.min(sapply(size, KSdist.s, x = x, prob = prob(distribution)))
 
             return(list(size = size[ind], prob = prob(distribution)))
         }
@@ -70,24 +70,24 @@
                 if(param[2] <= 0) return(Inf)
                 return(ks.test(x, "pnorm", mean = param[1], sd = param[2])$statistic)
             }
-            res <- optim(c(mean(distribution), sd(distribution)), f = KSdist, 
+            res <- optim(c(mean(distribution), sd(distribution)), fn = KSdist,
                         method = "Nelder-Mead", control=list(reltol=eps), 
                         x = x)$par
             return(list(mean = res[1], sd = res[2]))
         }
         if(param == "mean"){
-            KSdist <- function(mean, x, sd){
+            KSdist.m <- function(mean, x, sd){
                 return(ks.test(x, "pnorm", mean = mean, sd = sd)$statistic)
             }
-            res <- optimize(f = KSdist, interval = c(min(x), max(x)), 
+            res <- optimize(f = KSdist.m, interval = c(min(x), max(x)),
                         tol = eps, x = x, sd = sd(distribution))$minimum
             return(list(mean = res, sd = sd(distribution)))
         }
         if(param == "sd"){
-            KSdist <- function(sd, x, mean){
+            KSdist.s <- function(sd, x, mean){
                 return(ks.test(x, "pnorm", mean = mean, sd = sd)$statistic)
             }
-            res <- optimize(f = KSdist, 
+            res <- optimize(f = KSdist.s,
                         interval = c(.Machine$double.eps^0.5, max(x)-min(x)), 
                         tol = eps, x = x, mean = mean(distribution))$minimum
             return(list(mean = mean(distribution), sd = res))
@@ -105,24 +105,24 @@
                 if(param[2] <= 0) return(Inf)
                 return(ks.test(x, "plnorm", meanlog = param[1], sdlog = param[2])$statistic)
             }
-            res <- optim(c(meanlog(distribution), sdlog(distribution)), f = KSdist, 
+            res <- optim(c(meanlog(distribution), sdlog(distribution)), fn = KSdist,
                         method = "Nelder-Mead", control=list(reltol=eps), 
                         x = x)$par
             return(list(meanlog = res[1], sdlog = res[2]))
         }
         if(param == "meanlog"){
-            KSdist <- function(meanlog, x, sdlog){
+            KSdist.m <- function(meanlog, x, sdlog){
                 return(ks.test(x, "plnorm", meanlog = meanlog, sdlog = sdlog)$statistic)
             }
-            res <- optimize(f = KSdist, interval = c(min(x), max(x)), 
+            res <- optimize(f = KSdist.m, interval = c(min(x), max(x)),
                         tol = eps, x = x, sdlog = sdlog(distribution))$minimum
             return(list(meanlog = res, sdlog = sdlog(distribution)))
         }
         if(param == "sdlog"){
-            KSdist <- function(sdlog, x, meanlog){
+            KSdist.s <- function(sdlog, x, meanlog){
                 return(ks.test(x, "plnorm", meanlog = meanlog, sdlog = sdlog)$statistic)
             }
-            res <- optimize(f = KSdist, 
+            res <- optimize(f = KSdist.s,
                         interval = c(.Machine$double.eps^0.5, max(x)-min(x)), 
                         tol = eps, x = x, meanlog = meanlog(distribution))$minimum
             return(list(meanlog = meanlog(distribution), sdlog = res))
@@ -140,24 +140,24 @@
                 if(param[2] <= 0) return(Inf)
                 return(ks.test(x, "pgumbel", loc = param[1], scale = param[2])$statistic)
             }
-            res <- optim(c(mean(distribution), sd(distribution)), f = KSdist, 
+            res <- optim(c(mean(distribution), sd(distribution)), fn = KSdist,
                         method = "Nelder-Mead", control=list(reltol=eps), 
                         x = x)$par
             return(list(loc = res[1], scale = res[2]))
         }
         if(param == "loc"){
-            KSdist <- function(loc, x, scale){
+            KSdist.l <- function(loc, x, scale){
                 return(ks.test(x, "pgumbel", loc = loc, scale = scale)$statistic)
             }
-            res <- optimize(f = KSdist, interval = c(min(x), max(x)), 
+            res <- optimize(f = KSdist.l, interval = c(min(x), max(x)),
                         tol = eps, x = x, scale = scale(distribution))$minimum
             return(list(loc = res, scale = scale(distribution)))
         }
         if(param == "scale"){
-            KSdist <- function(scale, x, loc){
+            KSdist.s <- function(scale, x, loc){
                 return(ks.test(x, "pgumbel", loc = loc, scale = scale)$statistic)
             }
-            res <- optimize(f = KSdist, 
+            res <- optimize(f = KSdist.s,
                         interval = c(.Machine$double.eps^0.5, max(x)-min(x)), 
                         tol = eps, x = x, loc = loc(distribution))$minimum
             return(list(loc = loc(distribution), scale = res))
@@ -187,23 +187,23 @@
                 if((param[1] <= 0) || (param[2] <= 0)) return(Inf)
                 return(ks.test(x, "pgamma", scale = param[1], shape = param[2])$statistic)
             }
-            res <- optim(c(scale(distribution), shape(distribution)), f = KSdist, method = "Nelder-Mead", 
+            res <- optim(c(scale(distribution), shape(distribution)), fn = KSdist, method = "Nelder-Mead",
                         control=list(reltol=eps), x = x)$par
             return(list(scale = res[1], shape = res[2]))
         }
         if(param == "scale"){
-            KSdist <- function(scale, x, shape){
+            KSdist.sc <- function(scale, x, shape){
                 return(ks.test(x, "pgamma", scale = scale, shape = shape)$statistic)
             }
-            res <- optimize(f = KSdist, interval = c(min(x), max(x)), 
+            res <- optimize(f = KSdist.sc, interval = c(min(x), max(x)),
                         tol = eps, x = x, shape = shape(distribution))$minimum
             return(list(scale = res, shape = shape(distribution)))
         }
         if(param == "shape"){
-            KSdist <- function(shape, x, scale){
+            KSdist.sh <- function(shape, x, scale){
                 return(ks.test(x, "pgamma", scale = scale, shape = shape)$statistic)
             }
-            res <- optimize(f = KSdist, 
+            res <- optimize(f = KSdist.sh,
                         interval = c(.Machine$double.eps^0.5, max(x)-min(x)), 
                         tol = eps, x = x, scale = scale(distribution))$minimum
             return(list(scale = scale(distribution), shape = res))

Modified: branches/robast-1.1/pkg/ROptEstOld/R/leastFavorableRadius.R
===================================================================
--- branches/robast-1.1/pkg/ROptEstOld/R/leastFavorableRadius.R	2018-07-17 11:20:51 UTC (rev 960)
+++ branches/robast-1.1/pkg/ROptEstOld/R/leastFavorableRadius.R	2018-07-17 12:26:22 UTC (rev 961)
@@ -15,7 +15,7 @@
 
         L2derivDim <- numberOfMaps(L2Fam at L2deriv)
         if(L2derivDim == 1){
-            leastFavFct <- function(r, L2Fam, neighbor, risk, rho, 
+            leastFavFct.1 <- function(r, L2Fam, neighbor, risk, rho,
                                     upper.b, MaxIter, eps, warn){
                 loRad <- r*rho
                 upRad <- r/rho
@@ -51,16 +51,23 @@
                                         neighbor = neighbor, clip = resUp$b, cent = resUp$a, 
                                         stand = resUp$A, trafo = L2Fam at param@trafo)[[1]]
                 }
-                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)$root
+
+                ineff <- NULL
+                getIneffDiff.1 <- function(x){
+                            res <- getIneffDiff(x, L2Fam = L2Fam, neighbor = neighbor,
+                              upper.b = upper.b, risk = risk, loRad = loRad, upRad = upRad,
+                              loRisk = loRisk, upRisk = upRisk, eps = .Machine$double.eps^0.25,
+                              MaxIter = MaxIter, warn = warn)
+                            ineff <<- res["ineff"]
+                            return(res["ineffDiff"])
+                }
+                leastFavR <- uniroot(getIneffDiff.1, lower = lower, upper = upper,
+                                tol = .Machine$double.eps^0.25)$root
                 options(ow)
                 cat("current radius:\t", r, "\tinefficiency:\t", ineff, "\n")
                 return(ineff)
             }
-            leastFavR <- optimize(leastFavFct, lower = 1e-4, upper = upRad, 
+            leastFavR <- optimize(leastFavFct.1, 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, 
@@ -91,7 +98,7 @@
                         L2derivDistrSymm <- new("DistrSymmList", L2)
                     }
                 }
-                leastFavFct <- function(r, L2Fam, neighbor, risk, rho, 
+                leastFavFct.p <- function(r, L2Fam, neighbor, risk, rho,
                                         z.start, A.start, upper.b, MaxIter, eps, warn){
                     loRad <- r*rho
                     upRad <- r/rho
@@ -132,18 +139,24 @@
                          upRisk <- getAsRisk(risk = risk, L2deriv = L2deriv, neighbor = neighbor, 
                                         clip = resUp$b, cent = resUp$a, stand = resUp$A, trafo = trafo)[[1]]
                     }
-                    leastFavR <- uniroot(getIneffDiff, lower = lower, upper = upper, 
-                                    tol = .Machine$double.eps^0.25, L2Fam = L2Fam, neighbor = neighbor, 
-                                    z.start = z.start, A.start = A.start, upper.b = upper.b, risk = risk, 
-                                    loRad = loRad, upRad = upRad, loRisk = loRisk, upRisk = upRisk,
-                                    eps = eps, MaxIter = MaxIter, warn = warn)$root
+                    ineff <- NULL
+                    getIneffDiff.p <- function(x){
+                            res <- getIneffDiff(x, L2Fam = L2Fam, neighbor = neighbor,
+                                z.start = z.start, A.start = A.start, upper.b = upper.b, risk = risk,
+                                loRad = loRad, upRad = upRad, loRisk = loRisk, upRisk = upRisk,
+                                eps = .Machine$double.eps^0.25, MaxIter = MaxIter, warn = warn)
+                            ineff <<- res["ineff"]
+                            return(res["ineffDiff"])
+                            }
+                    leastFavR <- uniroot(getIneffDiff.p, lower = lower, upper = upper,
+                            tol = .Machine$double.eps^0.25)$root
                     options(ow)
                     cat("current radius:\t", r, "\tinefficiency:\t", ineff, "\n")
                     return(ineff)
                 }
                 if(is.null(z.start)) z.start <- numeric(L2derivDim)
                 if(is.null(A.start)) A.start <- L2Fam at param@trafo
-                leastFavR <- optimize(leastFavFct, lower = 1e-4, upper = upRad, 
+                leastFavR <- optimize(leastFavFct.p, 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, 

Modified: branches/robast-1.1/pkg/ROptEstOld/R/radiusMinimaxIC.R
===================================================================
--- branches/robast-1.1/pkg/ROptEstOld/R/radiusMinimaxIC.R	2018-07-17 11:20:51 UTC (rev 960)
+++ branches/robast-1.1/pkg/ROptEstOld/R/radiusMinimaxIC.R	2018-07-17 12:26:22 UTC (rev 961)
@@ -50,11 +50,17 @@
                                     stand = resUp$A, trafo = L2Fam at param@trafo)[[1]]
             }
 
-            leastFavR <- uniroot(getIneffDiff, lower = lower, upper = upper, 
-                            tol = .Machine$double.eps^0.25, L2Fam = L2Fam, neighbor = neighbor, 
-                            upper.b = upper.b, risk = risk, loRad = loRad, upRad = upRad, 
-                            loRisk = loRisk, upRisk = upRisk, eps = tol, 
-                            MaxIter = maxiter, warn = warn)$root
+            ineff <- NULL
+            getIneffDiff.1 <- function(x){
+                           res <- getIneffDiff(x, L2Fam = L2Fam, neighbor = neighbor,
+                            upper.b = upper.b, risk = risk, loRad = loRad, upRad = upRad,
+                            loRisk = loRisk, upRisk = upRisk, eps = .Machine$double.eps^0.25,
+                            MaxIter = maxiter, warn = warn)
+                            ineff <<- res["ineff"]
+                            return(res["ineffDiff"])
+            }
+            leastFavR <- uniroot(getIneffDiff.1, lower = lower, upper = upper,
+                            tol = .Machine$double.eps^0.25)$root
             neighbor at radius <- leastFavR
             res <- getInfRobIC(L2deriv = L2Fam at L2derivDistr[[1]], neighbor = neighbor, 
                         risk = risk, symm = L2Fam at L2derivSymm[[1]],
@@ -130,11 +136,17 @@
                                 clip = resUp$b, cent = resUp$a, stand = resUp$A, trafo = trafo)[[1]]
                 }
 
-                leastFavR <- uniroot(getIneffDiff, lower = lower, upper = upper, 
-                                tol = .Machine$double.eps^0.25, L2Fam = L2Fam, neighbor = neighbor, 
-                                z.start = z.start, A.start = A.start, upper.b = upper.b, risk = risk, 
-                                loRad = loRad, upRad = upRad, loRisk = loRisk, upRisk = upRisk, 
-                                eps = tol, MaxIter = maxiter, warn = warn)$root
+                ineff <- NULL
+                getIneffDiff.p <- function(x){
+                            res <- getIneffDiff(x, L2Fam = L2Fam, neighbor = neighbor,
+                                z.start = z.start, A.start = A.start, upper.b = upper.b, risk = risk,
+                                loRad = loRad, upRad = upRad, loRisk = loRisk, upRisk = upRisk,
+                                eps = .Machine$double.eps^0.25, MaxIter = maxiter, warn = warn)
+                            ineff <<- res["ineff"]
+                            return(res["ineffDiff"])
+                            }
+                leastFavR <- uniroot(getIneffDiff.p, lower = lower, upper = upper,
+                            tol = .Machine$double.eps^0.25)$root
                 neighbor at radius <- leastFavR
                 res <- getInfRobIC(L2deriv = L2deriv, neighbor = neighbor, risk = risk, 
                             Distr = L2Fam at distribution, DistrSymm = L2Fam at distrSymm, 



More information about the Robast-commits mailing list