[Robast-commits] r1039 - in pkg/ROptEstOld: . R inst
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Jul 23 22:48:43 CEST 2018
Author: ruckdeschel
Date: 2018-07-23 22:48:42 +0200 (Mon, 23 Jul 2018)
New Revision: 1039
Modified:
pkg/ROptEstOld/DESCRIPTION
pkg/ROptEstOld/NAMESPACE
pkg/ROptEstOld/R/AllInitialize.R
pkg/ROptEstOld/R/AllPlot.R
pkg/ROptEstOld/R/ContIC.R
pkg/ROptEstOld/R/getAsRisk.R
pkg/ROptEstOld/R/getIneffDiff.R
pkg/ROptEstOld/R/getInfCent.R
pkg/ROptEstOld/R/getInfRobIC_asBias.R
pkg/ROptEstOld/R/getInfRobIC_asCov.R
pkg/ROptEstOld/R/getInfRobIC_asHampel.R
pkg/ROptEstOld/R/infoPlot.R
pkg/ROptEstOld/R/ksEstimator.R
pkg/ROptEstOld/R/leastFavorableRadius.R
pkg/ROptEstOld/R/lowerCaseRadius.R
pkg/ROptEstOld/R/radiusMinimaxIC.R
pkg/ROptEstOld/inst/NEWS
pkg/ROptEstOld/inst/TOBEDONE
Log:
[ROptEstOld] merged branch 1.1 to trunk
Modified: pkg/ROptEstOld/DESCRIPTION
===================================================================
--- pkg/ROptEstOld/DESCRIPTION 2018-07-23 20:47:19 UTC (rev 1038)
+++ pkg/ROptEstOld/DESCRIPTION 2018-07-23 20:48:42 UTC (rev 1039)
@@ -1,17 +1,15 @@
Package: ROptEstOld
-Version: 0.9.2
-Date: 2013-09-12
-Title: Optimally robust estimation - old version
+Version: 1.1.0
+Date: 2018-07-17
+Title: Optimally Robust Estimation - Old Version
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.4), RandVar(>= 0.9.2), evd
-Author: Matthias Kohl
-Maintainer: Matthias Kohl <Matthias.Kohl at stamats.de>
-LazyLoad: yes
+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("aut", "cre", "cph"), email="Matthias.Kohl at stamats.de")
ByteCompile: yes
License: LGPL-3
URL: http://robast.r-forge.r-project.org/
Encoding: latin1
LastChangedDate: {$LastChangedDate$}
LastChangedRevision: {$LastChangedRevision$}
-SVNRevision: 696
+VCS/SVNRevision: 940
Modified: pkg/ROptEstOld/NAMESPACE
===================================================================
--- pkg/ROptEstOld/NAMESPACE 2018-07-23 20:47:19 UTC (rev 1038)
+++ pkg/ROptEstOld/NAMESPACE 2018-07-23 20:48:42 UTC (rev 1039)
@@ -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: pkg/ROptEstOld/R/AllInitialize.R
===================================================================
--- pkg/ROptEstOld/R/AllInitialize.R 2018-07-23 20:47:19 UTC (rev 1038)
+++ pkg/ROptEstOld/R/AllInitialize.R 2018-07-23 20:48:42 UTC (rev 1039)
@@ -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: pkg/ROptEstOld/R/AllPlot.R
===================================================================
--- pkg/ROptEstOld/R/AllPlot.R 2018-07-23 20:47:19 UTC (rev 1038)
+++ pkg/ROptEstOld/R/AllPlot.R 2018-07-23 20:48:42 UTC (rev 1039)
@@ -13,8 +13,8 @@
plot(e1)
if(is(e1, "AbscontDistribution")){
- lower <- ifelse(is.finite(q(e1)(0)), q(e1)(0), q(e1)(getdistrOption("TruncQuantile")))
- upper <- ifelse(is.finite(q(e1)(1)), q(e1)(1), q(e1)(1 - getdistrOption("TruncQuantile")))
+ lower <- ifelse(is.finite(q.l(e1)(0)), q.l(e1)(0), q.l(e1)(getdistrOption("TruncQuantile")))
+ upper <- ifelse(is.finite(q.l(e1)(1)), q.l(e1)(1), q.l(e1)(1 - getdistrOption("TruncQuantile")))
h <- upper - lower
x.vec <- seq(from = lower - 0.1*h, to = upper + 0.1*h, length = 1000)
plty <- "l"
@@ -69,8 +69,8 @@
if(!is(e1, "UnivariateDistribution")) stop("not yet implemented")
if(is(e1, "AbscontDistribution")){
- lower <- ifelse(is.finite(q(e1)(0)), q(e1)(0), q(e1)(getdistrOption("TruncQuantile")))
- upper <- ifelse(is.finite(q(e1)(1)), q(e1)(1), q(e1)(1 - getdistrOption("TruncQuantile")))
+ lower <- ifelse(is.finite(q.l(e1)(0)), q.l(e1)(0), q.l(e1)(getdistrOption("TruncQuantile")))
+ upper <- ifelse(is.finite(q.l(e1)(1)), q.l(e1)(1), q.l(e1)(1 - getdistrOption("TruncQuantile")))
h <- upper - lower
x.vec <- seq(from = lower - 0.1*h, to = upper + 0.1*h, length = 1000)
plty <- "l"
Modified: pkg/ROptEstOld/R/ContIC.R
===================================================================
--- pkg/ROptEstOld/R/ContIC.R 2018-07-23 20:47:19 UTC (rev 1038)
+++ pkg/ROptEstOld/R/ContIC.R 2018-07-23 20:48:42 UTC (rev 1039)
@@ -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: pkg/ROptEstOld/R/getAsRisk.R
===================================================================
--- pkg/ROptEstOld/R/getAsRisk.R 2018-07-23 20:47:19 UTC (rev 1038)
+++ pkg/ROptEstOld/R/getAsRisk.R 2018-07-23 20:48:42 UTC (rev 1039)
@@ -29,7 +29,7 @@
L2deriv = "UnivariateDistribution",
neighbor = "ContNeighborhood"),
function(risk, L2deriv, neighbor, trafo){
- z <- q(L2deriv)(0.5)
+ z <- q.l(L2deriv)(0.5)
bias <- abs(as.vector(trafo))/E(L2deriv, function(x, z){abs(x - z)},
useApply = FALSE, z = z)
@@ -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: pkg/ROptEstOld/R/getIneffDiff.R
===================================================================
--- pkg/ROptEstOld/R/getIneffDiff.R 2018-07-23 20:47:19 UTC (rev 1038)
+++ pkg/ROptEstOld/R/getIneffDiff.R 2018-07-23 20:48:42 UTC (rev 1039)
@@ -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: pkg/ROptEstOld/R/getInfCent.R
===================================================================
--- pkg/ROptEstOld/R/getInfCent.R 2018-07-23 20:47:19 UTC (rev 1038)
+++ pkg/ROptEstOld/R/getInfCent.R 2018-07-23 20:48:42 UTC (rev 1039)
@@ -9,8 +9,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 <- q.l(L2deriv)(getdistrOption("TruncQuantile"))
+ upper <- q.l(L2deriv)(1-getdistrOption("TruncQuantile"))
return(uniroot(z.fct, lower = lower, upper = upper, tol = tol.z,
c0=clip, D1=L2deriv)$root)
@@ -24,8 +24,8 @@
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))
}
- lower <- q(L2deriv)(getdistrOption("TruncQuantile"))
- upper <- q(L2deriv)(1-getdistrOption("TruncQuantile"))
+ lower <- q.l(L2deriv)(getdistrOption("TruncQuantile"))
+ upper <- q.l(L2deriv)(1-getdistrOption("TruncQuantile"))
return(uniroot(g.fct, lower = lower, upper = upper, tol = tol.z,
c0 = clip, D1 = D1)$root)
Modified: pkg/ROptEstOld/R/getInfRobIC_asBias.R
===================================================================
--- pkg/ROptEstOld/R/getInfRobIC_asBias.R 2018-07-23 20:47:19 UTC (rev 1038)
+++ pkg/ROptEstOld/R/getInfRobIC_asBias.R 2018-07-23 20:48:42 UTC (rev 1039)
@@ -8,7 +8,7 @@
upper, maxiter, tol, warn){
zi <- sign(as.vector(trafo))
A <- as.matrix(zi)
- z <- q(L2deriv)(0.5)
+ z <- q.l(L2deriv)(0.5)
b <- zi*as.vector(trafo)/E(L2deriv, function(x, z){abs(x - z)}, z = z)
if(is(L2deriv, "AbscontDistribution"))
Modified: pkg/ROptEstOld/R/getInfRobIC_asCov.R
===================================================================
--- pkg/ROptEstOld/R/getInfRobIC_asCov.R 2018-07-23 20:47:19 UTC (rev 1038)
+++ pkg/ROptEstOld/R/getInfRobIC_asCov.R 2018-07-23 20:48:42 UTC (rev 1039)
@@ -7,7 +7,7 @@
function(L2deriv, risk, neighbor, Finfo, trafo){
info <- c("optimal IC in sense of Cramer-Rao bound")
A <- trafo %*% solve(Finfo)
- b <- abs(as.vector(A))*max(abs(q(L2deriv)(1)),abs(q(L2deriv)(0)))
+ b <- abs(as.vector(A))*max(abs(q.l(L2deriv)(1)),abs(q.l(L2deriv)(0)))
Risk <- list(asCov = A %*% t(trafo), asBias = b)
return(list(A = A, a = 0, b = b, d = NULL, risk = Risk, info = info))
@@ -18,7 +18,7 @@
function(L2deriv, risk, neighbor, Finfo, trafo){
info <- c("optimal IC in sense of Cramer-Rao bound")
A <- trafo %*% solve(Finfo)
- b <- abs(as.vector(A))*(q(L2deriv)(1)-q(L2deriv)(0))
+ b <- abs(as.vector(A))*(q.l(L2deriv)(1)-q.l(L2deriv)(0))
Risk <- list(asCov = A %*% t(trafo), asBias = b)
return(list(A = A, a = -b/2, b = b, d = NULL, risk = Risk, info = info))
@@ -31,8 +31,8 @@
A <- trafo %*% solve(Finfo)
IC <- A %*% L2deriv
if(is(Distr, "UnivariateDistribution")){
- lower <- ifelse(is.finite(q(Distr)(0)), q(Distr)(1e-8), q(Distr)(0))
- upper <- ifelse(is.finite(q(Distr)(1)), q(Distr)(1-1e-8), q(Distr)(1))
+ lower <- ifelse(is.finite(q.l(Distr)(0)), q.l(Distr)(1e-8), q.l(Distr)(0))
+ upper <- ifelse(is.finite(q.l(Distr)(1)), q.l(Distr)(1-1e-8), q.l(Distr)(1))
x <- seq(from = lower, to = upper, length = 1e5)
x <- x[x!=0] # problems with NaN=log(0)!
b <- evalRandVar(IC, as.matrix(x))^2
Modified: pkg/ROptEstOld/R/getInfRobIC_asHampel.R
===================================================================
--- pkg/ROptEstOld/R/getInfRobIC_asHampel.R 2018-07-23 20:47:19 UTC (rev 1038)
+++ pkg/ROptEstOld/R/getInfRobIC_asHampel.R 2018-07-23 20:48:42 UTC (rev 1039)
@@ -8,7 +8,7 @@
upper, maxiter, tol, warn){
A <- trafo / E(L2deriv, function(x){x^2})
b <- risk at bound
- bmax <- abs(as.vector(A))*max(abs(q(L2deriv)(0)), q(L2deriv)(1))
+ bmax <- abs(as.vector(A))*max(abs(q.l(L2deriv)(0)), q.l(L2deriv)(1))
if(b >= bmax){
if(warn) cat("'b >= maximum asymptotic bias' => (classical) optimal IC\n",
"in sense of Cramer-Rao bound is returned\n")
@@ -70,8 +70,8 @@
if(is.null(A.start)) A.start <- trafo
ClassIC <- trafo %*% solve(Finfo) %*% L2deriv
- lower <- q(Distr)(getdistrOption("TruncQuantile"))
- upper <- q(Distr)(1-getdistrOption("TruncQuantile"))
+ lower <- q.l(Distr)(getdistrOption("TruncQuantile"))
+ upper <- q.l(Distr)(1-getdistrOption("TruncQuantile"))
x <- seq(from = lower, to = upper, by = 0.01)
bmax <- evalRandVar(ClassIC, as.matrix(x))^2
bmax <- sqrt(max(colSums(bmax)))
Modified: pkg/ROptEstOld/R/infoPlot.R
===================================================================
--- pkg/ROptEstOld/R/infoPlot.R 2018-07-23 20:47:19 UTC (rev 1038)
+++ pkg/ROptEstOld/R/infoPlot.R 2018-07-23 20:48:42 UTC (rev 1039)
@@ -7,8 +7,8 @@
if(is(e1, "UnivariateDistribution")){
if(is(e1, "AbscontDistribution")){
- ifelse(is.finite(q(e1)(0)), lower <- q(e1)(0), lower <- q(e1)(getdistrOption("TruncQuantile")))
- ifelse(is.finite(q(e1)(1)), upper <- q(e1)(1), upper <- q(e1)(1 - getdistrOption("TruncQuantile")))
+ ifelse(is.finite(q.l(e1)(0)), lower <- q.l(e1)(0), lower <- q.l(e1)(getdistrOption("TruncQuantile")))
+ ifelse(is.finite(q.l(e1)(1)), upper <- q.l(e1)(1), upper <- q.l(e1)(1 - getdistrOption("TruncQuantile")))
h <- upper - lower
x.vec <- seq(from = lower - 0.1*h, to = upper + 0.1*h, length = 1000)
plty <- "l"
Modified: pkg/ROptEstOld/R/ksEstimator.R
===================================================================
--- pkg/ROptEstOld/R/ksEstimator.R 2018-07-23 20:47:19 UTC (rev 1038)
+++ pkg/ROptEstOld/R/ksEstimator.R 2018-07-23 20:48:42 UTC (rev 1039)
@@ -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: pkg/ROptEstOld/R/leastFavorableRadius.R
===================================================================
--- pkg/ROptEstOld/R/leastFavorableRadius.R 2018-07-23 20:47:19 UTC (rev 1038)
+++ pkg/ROptEstOld/R/leastFavorableRadius.R 2018-07-23 20:48:42 UTC (rev 1039)
@@ -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: pkg/ROptEstOld/R/lowerCaseRadius.R
===================================================================
--- pkg/ROptEstOld/R/lowerCaseRadius.R 2018-07-23 20:47:19 UTC (rev 1038)
+++ pkg/ROptEstOld/R/lowerCaseRadius.R 2018-07-23 20:48:42 UTC (rev 1039)
@@ -13,7 +13,7 @@
w0 <- options("warn")
options(warn = -1)
L2deriv <- L2Fam at L2derivDistr[[1]]
- m <- q(L2deriv)(0.5)
+ m <- q.l(L2deriv)(0.5)
wsm <- d(L2deriv)(m)
supp <- support(L2deriv)
Modified: pkg/ROptEstOld/R/radiusMinimaxIC.R
===================================================================
--- pkg/ROptEstOld/R/radiusMinimaxIC.R 2018-07-23 20:47:19 UTC (rev 1038)
+++ pkg/ROptEstOld/R/radiusMinimaxIC.R 2018-07-23 20:48:42 UTC (rev 1039)
@@ -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,
Modified: pkg/ROptEstOld/inst/NEWS
===================================================================
--- pkg/ROptEstOld/inst/NEWS 2018-07-23 20:47:19 UTC (rev 1038)
+++ pkg/ROptEstOld/inst/NEWS 2018-07-23 20:48:42 UTC (rev 1039)
@@ -8,6 +8,24 @@
information)
#######################################
+version 1.1
+#######################################
+
+user-visible CHANGES:
++ DESCRIPTION tag SVNRevision changed to VCS/SVNRevision
+
+under the hood:
++ wherever possible also use q.l internally instead of q to
+ provide functionality in IRKernel
+
+#######################################
+version 1.0
+#######################################
+
+user-visible CHANGES:
++ title changed to title style / capitalization
+
+#######################################
version 0.9
#######################################
Modified: pkg/ROptEstOld/inst/TOBEDONE
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/robast -r 1039
More information about the Robast-commits
mailing list