[Robast-commits] r314 - in branches/robast-0.7/pkg/RobLox: . R chm inst man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Jun 30 15:05:57 CEST 2009
Author: stamats
Date: 2009-06-30 15:05:57 +0200 (Tue, 30 Jun 2009)
New Revision: 314
Added:
branches/robast-0.7/pkg/RobLox/R/finiteSampleCorrection.R
branches/robast-0.7/pkg/RobLox/inst/CITATION
branches/robast-0.7/pkg/RobLox/inst/NEWS
branches/robast-0.7/pkg/RobLox/man/0RobLox-package.Rd
branches/robast-0.7/pkg/RobLox/man/finiteSampleCorrection.Rd
Modified:
branches/robast-0.7/pkg/RobLox/DESCRIPTION
branches/robast-0.7/pkg/RobLox/NAMESPACE
branches/robast-0.7/pkg/RobLox/R/colRoblox.R
branches/robast-0.7/pkg/RobLox/R/rlsOptIC_TuMad.R
branches/robast-0.7/pkg/RobLox/R/roblox.R
branches/robast-0.7/pkg/RobLox/R/rowRoblox.R
branches/robast-0.7/pkg/RobLox/R/sysdata.rda
branches/robast-0.7/pkg/RobLox/chm/00Index.html
branches/robast-0.7/pkg/RobLox/chm/RobLox.chm
branches/robast-0.7/pkg/RobLox/chm/RobLox.toc
branches/robast-0.7/pkg/RobLox/chm/rlOptIC.html
branches/robast-0.7/pkg/RobLox/chm/rlsOptIC.AL.html
branches/robast-0.7/pkg/RobLox/chm/rlsOptIC.An1.html
branches/robast-0.7/pkg/RobLox/chm/rlsOptIC.An2.html
branches/robast-0.7/pkg/RobLox/chm/rlsOptIC.AnMad.html
branches/robast-0.7/pkg/RobLox/chm/rlsOptIC.BM.html
branches/robast-0.7/pkg/RobLox/chm/rlsOptIC.Ha3.html
branches/robast-0.7/pkg/RobLox/chm/rlsOptIC.Ha4.html
branches/robast-0.7/pkg/RobLox/chm/rlsOptIC.HaMad.html
branches/robast-0.7/pkg/RobLox/chm/rlsOptIC.Hu1.html
branches/robast-0.7/pkg/RobLox/chm/rlsOptIC.Hu2.html
branches/robast-0.7/pkg/RobLox/chm/rlsOptIC.Hu2a.html
branches/robast-0.7/pkg/RobLox/chm/rlsOptIC.Hu3.html
branches/robast-0.7/pkg/RobLox/chm/rlsOptIC.HuMad.html
branches/robast-0.7/pkg/RobLox/chm/rlsOptIC.M.html
branches/robast-0.7/pkg/RobLox/chm/rlsOptIC.MM2.html
branches/robast-0.7/pkg/RobLox/chm/rlsOptIC.Tu1.html
branches/robast-0.7/pkg/RobLox/chm/rlsOptIC.Tu2.html
branches/robast-0.7/pkg/RobLox/chm/rlsOptIC.TuMad.html
branches/robast-0.7/pkg/RobLox/chm/roblox.html
branches/robast-0.7/pkg/RobLox/chm/rowRoblox.html
branches/robast-0.7/pkg/RobLox/chm/rsOptIC.html
branches/robast-0.7/pkg/RobLox/man/roblox.Rd
branches/robast-0.7/pkg/RobLox/man/rowRoblox.Rd
Log:
merged branch and trunk
Modified: branches/robast-0.7/pkg/RobLox/DESCRIPTION
===================================================================
--- branches/robast-0.7/pkg/RobLox/DESCRIPTION 2009-06-30 13:02:54 UTC (rev 313)
+++ branches/robast-0.7/pkg/RobLox/DESCRIPTION 2009-06-30 13:05:57 UTC (rev 314)
@@ -1,10 +1,9 @@
Package: RobLox
Version: 0.7
-Date: 2008-11-27
-Title: Optimally robust influence curves for location and scale
+Date: 2009-04-21
+Title: Optimally robust influence curves and estimators for location and scale
Description: Functions for the determination of optimally robust
- influence curves and estimators in case of normal location with
- unknown scale
+ influence curves and estimators in case of normal location and/or scale
Depends: R(>= 2.7.0), stats, distrMod(>= 2.0.1), RobAStBase(>= 0.1.1)
Suggests: Biobase
Author: Matthias Kohl <Matthias.Kohl at stamats.de>
Modified: branches/robast-0.7/pkg/RobLox/NAMESPACE
===================================================================
--- branches/robast-0.7/pkg/RobLox/NAMESPACE 2009-06-30 13:02:54 UTC (rev 313)
+++ branches/robast-0.7/pkg/RobLox/NAMESPACE 2009-06-30 13:05:57 UTC (rev 314)
@@ -1,4 +1,5 @@
-export(rlsOptIC.AL,
+export(finiteSampleCorrection,
+ rlsOptIC.AL,
rlsOptIC.M,
rlsOptIC.BM,
rlsOptIC.Hu1,
Modified: branches/robast-0.7/pkg/RobLox/R/colRoblox.R
===================================================================
--- branches/robast-0.7/pkg/RobLox/R/colRoblox.R 2009-06-30 13:02:54 UTC (rev 313)
+++ branches/robast-0.7/pkg/RobLox/R/colRoblox.R 2009-06-30 13:05:57 UTC (rev 314)
@@ -1,7 +1,8 @@
###############################################################################
## Evaluate roblox on columns of a matrix
###############################################################################
-colRoblox <- function(x, mean, sd, eps, eps.lower, eps.upper, initial.est, k = 1){
+colRoblox <- function(x, mean, sd, eps, eps.lower, eps.upper, initial.est,
+ k = 1L, fsCor = TRUE, mad0 = 1e-4){
call.est <- match.call()
if(missing(x))
stop("'x' is missing with no default")
@@ -14,7 +15,8 @@
or 'data.matrix'")
res <- rowRoblox(x = t(x), mean = mean, sd = sd, eps = eps, eps.lower = eps.lower,
- eps.upper = eps.upper, initial.est = initial.est, k = k)
+ eps.upper = eps.upper, initial.est = initial.est, k = k,
+ fsCor = fsCor, mad0 = mad0)
res at estimate.call <- call.est
return(res)
}
Copied: branches/robast-0.7/pkg/RobLox/R/finiteSampleCorrection.R (from rev 308, pkg/RobLox/R/finiteSampleCorrection.R)
===================================================================
--- branches/robast-0.7/pkg/RobLox/R/finiteSampleCorrection.R (rev 0)
+++ branches/robast-0.7/pkg/RobLox/R/finiteSampleCorrection.R 2009-06-30 13:05:57 UTC (rev 314)
@@ -0,0 +1,27 @@
+###############################################################################
+## Function for finite-sample correction of the neighborhood radius
+###############################################################################
+finiteSampleCorrection <- function(r, n, model = "locsc"){
+ if(model == "locsc" & r >= 1.74) return(r)
+ if(model %in% c("loc", "sc") & r >= 3.0) return(r)
+ if(n == 1) return(Inf)
+ if(n == 2) return(Inf)
+
+ eps <- r/sqrt(n)
+ ns <- c(3:50, seq(55, 100, by = 5), seq(110, 200, by = 10),
+ seq(250, 500, by = 50))
+ epss <- c(seq(0.001, 0.01, by = 0.001), seq(0.02, to = 0.5, by = 0.01))
+ if(n %in% ns){
+ ind <- ns == n
+ }else{
+ ind <- which.min(abs(ns-n))
+ }
+ if(model == "locsc")
+ return(max(r, approx(x = epss, y = .finiteSampleRadius.locsc[,ind], xout = eps, rule = 2)$y))
+ if(model == "loc")
+ return(max(r, approx(x = epss, y = .finiteSampleRadius.loc[,ind], xout = eps, rule = 2)$y))
+ if(model == "sc")
+ return(max(r, approx(x = epss, y = .finiteSampleRadius.sc[,ind], xout = eps, rule = 2)$y))
+ else
+ stop("argument 'model' has to be 'locsc', 'loc' or 'sc'")
+}
Modified: branches/robast-0.7/pkg/RobLox/R/rlsOptIC_TuMad.R
===================================================================
--- branches/robast-0.7/pkg/RobLox/R/rlsOptIC_TuMad.R 2009-06-30 13:02:54 UTC (rev 313)
+++ branches/robast-0.7/pkg/RobLox/R/rlsOptIC_TuMad.R 2009-06-30 13:05:57 UTC (rev 314)
@@ -19,7 +19,7 @@
a.mad <- qnorm(0.75)
b.mad <- 1/(4*a.mad*dnorm(a.mad))
- b.2 <- sqrt(A.loc^2*(16/25/sqrt(5)*a^5)^2 + b.mad^2)
+ b.2 <- sqrt(A.loc^2*(16/25/sqrt(5)*a^5)^2 + b.mad^2)^2
return(.TuMadrlsGetvar(a = a) + r^2*b.2)
}
Modified: branches/robast-0.7/pkg/RobLox/R/roblox.R
===================================================================
--- branches/robast-0.7/pkg/RobLox/R/roblox.R 2009-06-30 13:02:54 UTC (rev 313)
+++ branches/robast-0.7/pkg/RobLox/R/roblox.R 2009-06-30 13:05:57 UTC (rev 314)
@@ -165,8 +165,8 @@
###############################################################################
## optimally robust estimator for normal location and/or scale
###############################################################################
-roblox <- function(x, mean, sd, eps, eps.lower, eps.upper, initial.est, k = 1,
- returnIC = FALSE){
+roblox <- function(x, mean, sd, eps, eps.lower, eps.upper, initial.est, k = 1L,
+ fsCor = TRUE, returnIC = FALSE, mad0 = 1e-4){
es.call <- match.call()
if(missing(x))
stop("'x' is missing with no default")
@@ -182,6 +182,46 @@
stop("number of rows and columns/variables > 1. Please, do use 'rowRoblox'
resp. 'colRoblox'.")
}
+ if(length(x) <= 2){
+ if(missing(mean) && missing(sd)){
+ warning("Sample size <= 2! => Median and MAD are used for estimation.")
+ robEst <- c(median(x, na.rm = TRUE), mad(x, na.rm = TRUE))
+ names(robEst) <- c("mean", "sd")
+ Info.matrix <- matrix(c("roblox",
+ paste("median and MAD")),
+ ncol = 2, dimnames = list(NULL, c("method", "message")))
+ return(new("ALEstimate", name = "Median and MAD",
+ estimate.call = es.call, estimate = robEst,
+ samplesize = length(x), asvar = NULL,
+ asbias = NULL, pIC = NULL, Infos = Info.matrix))
+ }
+ if(missing(mean)){
+ warning("Sample size <= 2! => Median is used for estimation.")
+ robEst <- median(x, na.rm = TRUE)
+ names(robEst) <- "mean"
+ Info.matrix <- matrix(c("roblox",
+ paste("median")),
+ ncol = 2, dimnames = list(NULL, c("method", "message")))
+ return(new("ALEstimate", name = "Median",
+ estimate.call = es.call, estimate = robEst,
+ samplesize = length(x), asvar = NULL,
+ asbias = NULL, pIC = NULL, Infos = Info.matrix))
+ }
+ if(missing(sd)){
+ warning("Sample size <= 2! => MAD is used for estimation.")
+ if(length(mean) != 1)
+ stop("mean has length != 1")
+ robEst <- mad(x, center = mean, na.rm = TRUE)
+ names(robEst) <- "sd"
+ Info.matrix <- matrix(c("roblox",
+ paste("MAD")),
+ ncol = 2, dimnames = list(NULL, c("method", "message")))
+ return(new("ALEstimate", name = "MAD",
+ estimate.call = es.call, estimate = robEst,
+ samplesize = length(x), asvar = NULL,
+ asbias = NULL, pIC = NULL, Infos = Info.matrix))
+ }
+ }
if(missing(eps) && missing(eps.lower) && missing(eps.upper)){
eps.lower <- 0
eps.upper <- 0.5
@@ -200,10 +240,48 @@
}else{
if(length(eps) != 1)
stop("'eps' has to be of length 1")
- if(eps == 0)
- stop("'eps = 0'! => use functions 'mean' and 'sd' for estimation")
if((eps < 0) || (eps > 0.5))
stop("'eps' has to be in (0, 0.5]")
+ if(eps == 0){
+ if(missing(mean) && missing(sd)){
+ warning("eps = 0! => Mean and sd are used for estimation.")
+ n <- sum(!is.na(x))
+ robEst <- c(mean(x, na.rm = TRUE), sqrt((n-1)/n)*sd(x, na.rm = TRUE))
+ names(robEst) <- c("mean", "sd")
+ Info.matrix <- matrix(c("roblox",
+ paste("mean and sd")),
+ ncol = 2, dimnames = list(NULL, c("method", "message")))
+ return(new("ALEstimate", name = "Mean and sd",
+ estimate.call = es.call, estimate = robEst,
+ samplesize = n, asvar = NULL,
+ asbias = NULL, pIC = NULL, Infos = Info.matrix))
+ }
+ if(missing(mean)){
+ warning("eps = 0! => Mean is used for estimation.")
+ robEst <- mean(x, na.rm = TRUE)
+ names(robEst) <- "mean"
+ Info.matrix <- matrix(c("roblox",
+ paste("mean")),
+ ncol = 2, dimnames = list(NULL, c("method", "message")))
+ return(new("ALEstimate", name = "Mean",
+ estimate.call = es.call, estimate = robEst,
+ samplesize = length(x), asvar = NULL,
+ asbias = NULL, pIC = NULL, Infos = Info.matrix))
+ }
+ if(missing(sd)){
+ warning("eps = 0! => sd is used for estimation.")
+ n <- sum(!is.na(x))
+ robEst <- sqrt((n-1)/n)*sd(x, na.rm = TRUE)
+ names(robEst) <- "sd"
+ Info.matrix <- matrix(c("roblox",
+ paste("sd")),
+ ncol = 2, dimnames = list(NULL, c("method", "message")))
+ return(new("ALEstimate", name = "sd",
+ estimate.call = es.call, estimate = robEst,
+ samplesize = n, asvar = NULL,
+ asbias = NULL, pIC = NULL, Infos = Info.matrix))
+ }
+ }
}
if(!is.integer(k))
k <- as.integer(k)
@@ -213,17 +291,20 @@
if(length(k) != 1){
stop("'k' has to be of length 1")
}
-
if(missing(mean) && missing(sd)){
if(missing(initial.est)){
mean <- median(x, na.rm = TRUE)
- sd <- mad(x, na.rm = TRUE)
- if(sd == 0)
- stop("'mad(x, na.rm = TRUE) == 0' => cannot compute a valid initial estimate,
- please specify one via 'initial.est'")
+ sd <- mad(x, center = mean, na.rm = TRUE)
+ if(sd == 0){
+ warning("'mad(x, na.rm = TRUE) = 0' => cannot compute a valid initial estimate.
+ To avoid division by zero 'mad0' is used. You could also specify
+ a valid scale estimate via 'initial.est'. Note that you have to provide
+ a location and scale estimate.")
+ sd <- mad0
+ }
}else{
if(!is.numeric(initial.est) || length(initial.est) != 2)
- stop("'initial.est' needs to be a numeric vector of length 2 or missing")
+ stop("'initial.est' needs to be a numeric vector of length 2 or missing")
mean <- initial.est[1]
sd <- initial.est[2]
if(sd <= 0)
@@ -232,6 +313,7 @@
if(!missing(eps)){
r <- sqrt(length(x))*eps
+ if(fsCor) r <- finiteSampleCorrection(r = r, n = length(x), model = "locsc")
if(r > 10){
b <- sd*1.618128043
const <- 1.263094656
@@ -248,10 +330,17 @@
}
robEst <- .kstep.locsc(x = x, initial.est = c(mean, sd), A1 = A1, A2 = A2, a = a, b = b, k = k)
names(robEst$est) <- c("mean", "sd")
- Info.matrix <- matrix(c("roblox",
- paste("optimally robust estimate for contamination 'eps' =", round(eps, 3),
- "and 'asMSE'")),
- ncol = 2, dimnames = list(NULL, c("method", "message")))
+ if(fsCor){
+ Info.matrix <- matrix(c("roblox",
+ paste("finite-sample corrected optimally robust estimate for contamination 'eps' =",
+ round(eps, 3), "and 'asMSE'")),
+ ncol = 2, dimnames = list(NULL, c("method", "message")))
+ }else{
+ Info.matrix <- matrix(c("roblox",
+ paste("optimally robust estimate for contamination 'eps' =", round(eps, 3),
+ "and 'asMSE'")),
+ ncol = 2, dimnames = list(NULL, c("method", "message")))
+ }
if(returnIC){
w <- new("HampelWeight")
clip(w) <- robEst$b
@@ -299,6 +388,7 @@
}
L2Fam <- substitute(NormLocationScaleFamily(mean = m1, sd = s1),
list(m1 = robEst$est[1], s1 = robEst$est[2]))
+ info <- c("roblox", "optimally robust IC for AL estimators and 'asMSE'")
IC1 <- generateIC(neighbor = ContNeighborhood(radius = r),
L2Fam = eval(L2Fam),
res = list(A = diag(c(robEst$A1, robEst$A2)), a = robEst$a,
@@ -306,9 +396,9 @@
risk = list(asMSE = mse, asBias = robEst$b,
trAsCov = mse - r^2*robEst$b^2,
asCov = robEst$asVar),
- info = c("roblox", "optimally robust IC for AL estimators and 'asMSE'"),
- w = w, biastype = symmetricBias(), normtype = NormType(),
- modifyIC = modIC))
+ info = info, w = w, biastype = symmetricBias(),
+ normtype = NormType(), modifyIC = modIC))
+ Infos(IC1) <- Info.matrix
return(new("kStepEstimate", name = "Optimally robust estimate",
estimate.call = es.call, estimate = robEst$est,
samplesize = length(x), asvar = robEst$asvar,
@@ -328,6 +418,10 @@
r <- uniroot(.getlsInterval, lower = rlo+1e-8, upper = rup,
tol = .Machine$double.eps^0.25, rlo = rlo, rup = rup)$root
}
+ if(fsCor){
+ r.as <- r
+ r <- finiteSampleCorrection(r = r, n = length(x), model = "locsc")
+ }
if(r > 10){
b <- sd*1.618128043
const <- 1.263094656
@@ -343,24 +437,33 @@
mse <- A1 + A2
}
if(rlo == 0){
- ineff <- (A1 + A2 - b^2*r^2)/(1.5*sd^2)
+ ineff <- (A1 + A2 - b^2*r.as^2)/(1.5*sd^2)
}else{
if(rlo > 10){
ineff <- 1
}else{
A1lo <- sd^2*.getA1.locsc(rlo)
A2lo <- sd^2*.getA2.locsc(rlo)
- ineff <- (A1 + A2 - b^2*(r^2 - rlo^2))/(A1lo + A2lo)
+ ineff <- (A1 + A2 - b^2*(r.as^2 - rlo^2))/(A1lo + A2lo)
}
}
robEst <- .kstep.locsc(x = x, initial.est = c(mean, sd), A1 = A1, A2 = A2, a = a, b = b, k = k)
names(robEst$est) <- c("mean", "sd")
- Info.matrix <- matrix(c(rep("roblox", 3),
- paste("radius-minimax estimate for contamination interval [",
- round(eps.lower, 3), ", ", round(eps.upper, 3), "]", sep = ""),
- paste("least favorable contamination: ", round(r/sqrtn, 3), sep = ""),
- paste("maximum MSE-inefficiency: ", round(ineff, 3), sep = "")),
- ncol = 2, dimnames = list(NULL, c("method", "message")))
+ if(fsCor){
+ Info.matrix <- matrix(c(rep("roblox", 3),
+ paste("finite-sample corrected radius-minimax estimate for contamination interval [",
+ round(eps.lower, 3), ", ", round(eps.upper, 3), "]", sep = ""),
+ paste("least favorable (uncorrected) contamination: ", round(r.as/sqrtn, 3), sep = ""),
+ paste("maximum asymptotic MSE-inefficiency: ", round(ineff, 3), sep = "")),
+ ncol = 2, dimnames = list(NULL, c("method", "message")))
+ }else{
+ Info.matrix <- matrix(c(rep("roblox", 3),
+ paste("radius-minimax estimate for contamination interval [",
+ round(eps.lower, 3), ", ", round(eps.upper, 3), "]", sep = ""),
+ paste("least favorable contamination: ", round(r/sqrtn, 3), sep = ""),
+ paste("maximum asymptotic MSE-inefficiency: ", round(ineff, 3), sep = "")),
+ ncol = 2, dimnames = list(NULL, c("method", "message")))
+ }
if(returnIC){
w <- new("HampelWeight")
clip(w) <- robEst$b
@@ -408,6 +511,7 @@
}
L2Fam <- substitute(NormLocationScaleFamily(mean = m1, sd = s1),
list(m1 = robEst$est[1], s1 = robEst$est[2]))
+ info <- c("roblox", "optimally robust IC for AL estimators and 'asMSE'")
IC1 <- generateIC(neighbor = ContNeighborhood(radius = r),
L2Fam = eval(L2Fam),
res = list(A = diag(c(robEst$A1, robEst$A2)), a = robEst$a,
@@ -415,15 +519,9 @@
risk = list(asMSE = mse, asBias = robEst$b,
trAsCov = mse - r^2*robEst$b^2,
asCov = robEst$asvar),
- info = c("roblox", "optimally robust IC for AL estimators and 'asMSE'"),
- w = w, biastype = symmetricBias(), normtype = NormType(),
- modifyIC = modIC))
- Infos(IC1) <- matrix(c(rep("roblox", 3),
- paste("radius-minimax IC for contamination interval [",
- round(eps.lower, 3), ", ", round(eps.upper, 3), "]", sep = ""),
- paste("least favorable contamination: ", round(r/sqrtn, 3), sep = ""),
- paste("maximum MSE-inefficiency: ", round(ineff, 3), sep = "")),
- ncol = 2, dimnames = list(NULL, c("method", "message")))
+ info = info, w = w, biastype = symmetricBias(),
+ normtype = NormType(), modifyIC = modIC))
+ Infos(IC1) <- Info.matrix
return(new("kStepEstimate", name = "Optimally robust estimate",
estimate.call = es.call, estimate = robEst$est,
samplesize = length(x), asvar = robEst$asvar,
@@ -450,6 +548,7 @@
if(!missing(eps)){
r <- sqrt(length(x))*eps
+ if(fsCor) r <- finiteSampleCorrection(r = r, n = length(x), model = "loc")
if(r > 10){
b <- sd*sqrt(pi/2)
A <- b^2*(1+r^2)
@@ -459,10 +558,17 @@
}
robEst <- .kstep.loc(x = x, initial.est = mean, A = A, b = b, sd = sd, k = k)
names(robEst) <- "mean"
- Info.matrix <- matrix(c("roblox",
- paste("optimally robust estimate for contamination 'eps' =", round(eps, 3),
- "and 'asMSE'")),
- ncol = 2, dimnames = list(NULL, c("method", "message")))
+ if(fsCor){
+ Info.matrix <- matrix(c("roblox",
+ paste("finite-sample corrected optimally robust estimate for contamination 'eps' =", round(eps, 3),
+ "and 'asMSE'")),
+ ncol = 2, dimnames = list(NULL, c("method", "message")))
+ }else{
+ Info.matrix <- matrix(c("roblox",
+ paste("optimally robust estimate for contamination 'eps' =", round(eps, 3),
+ "and 'asMSE'")),
+ ncol = 2, dimnames = list(NULL, c("method", "message")))
+ }
if(returnIC){
w <- new("HampelWeight")
clip(w) <- b
@@ -490,6 +596,7 @@
info = c("roblox", "optimally robust IC for AL estimators and 'asMSE'"),
w = w, biastype = symmetricBias(), normtype = NormType(),
modifyIC = modIC))
+ Infos(IC1) <- Info.matrix
return(new("kStepEstimate", name = "Optimally robust estimate",
estimate.call = es.call, estimate = robEst,
samplesize = length(x), asvar = as.matrix(A-r^2*b^2),
@@ -509,6 +616,10 @@
r <- uniroot(.getlInterval, lower = rlo+1e-8, upper = rup,
tol = .Machine$double.eps^0.25, rlo = rlo, rup = rup)$root
}
+ if(fsCor){
+ r.as <- r
+ r <- finiteSampleCorrection(r = r, n = length(x), model = "loc")
+ }
if(r > 10){
b <- sd*sqrt(pi/2)
A <- b^2*(1+r^2)
@@ -528,12 +639,21 @@
}
robEst <- .kstep.loc(x = x, initial.est = mean, A = A, b = b, sd = sd, k = k)
names(robEst) <- "mean"
- Info.matrix <- matrix(c(rep("roblox", 3),
- paste("radius-minimax estimate for contamination interval [",
- round(eps.lower, 3), ", ", round(eps.upper, 3), "]", sep = ""),
- paste("least favorable contamination: ", round(r/sqrtn, 3), sep = ""),
- paste("maximum MSE-inefficiency: ", round(ineff, 3), sep = "")),
- ncol = 2, dimnames = list(NULL, c("method", "message")))
+ if(fsCor){
+ Info.matrix <- matrix(c(rep("roblox", 3),
+ paste("finite-sample corrected radius-minimax estimate for contamination interval [",
+ round(eps.lower, 3), ", ", round(eps.upper, 3), "]", sep = ""),
+ paste("least favorable (uncorrected) contamination: ", round(r.as/sqrtn, 3), sep = ""),
+ paste("maximum asymptotic MSE-inefficiency: ", round(ineff, 3), sep = "")),
+ ncol = 2, dimnames = list(NULL, c("method", "message")))
+ }else{
+ Info.matrix <- matrix(c(rep("roblox", 3),
+ paste("radius-minimax estimate for contamination interval [",
+ round(eps.lower, 3), ", ", round(eps.upper, 3), "]", sep = ""),
+ paste("least favorable contamination: ", round(r/sqrtn, 3), sep = ""),
+ paste("maximum MSE-inefficiency: ", round(ineff, 3), sep = "")),
+ ncol = 2, dimnames = list(NULL, c("method", "message")))
+ }
if(returnIC){
w <- new("HampelWeight")
clip(w) <- b
@@ -561,12 +681,7 @@
info = c("roblox", "optimally robust IC for AL estimators and 'asMSE'"),
w = w, biastype = symmetricBias(), normtype = NormType(),
modifyIC = modIC))
- Infos(IC1) <- matrix(c(rep("roblox", 3),
- paste("radius-minimax IC for contamination interval [",
- round(eps.lower, 3), ", ", round(eps.upper, 3), "]", sep = ""),
- paste("least favorable contamination: ", round(r/sqrtn, 3), sep = ""),
- paste("maximum MSE-inefficiency: ", round(ineff, 3), sep = "")),
- ncol = 2, dimnames = list(NULL, c("method", "message")))
+ Infos(IC1) <- Info.matrix
return(new("kStepEstimate", name = "Optimally robust estimate",
estimate.call = es.call, estimate = robEst,
samplesize = length(x), asvar = as.matrix(A-r^2*b^2),
@@ -582,10 +697,13 @@
if(length(mean) != 1)
stop("mean has length != 1")
if(missing(initial.est)){
- sd <- mad(x, na.rm = TRUE)
- if(sd == 0)
- stop("'mad(x, na.rm = TRUE) == 0' => cannot compute a valid initial estimate,
- please specify one via 'initial.est'")
+ sd <- mad(x, center = mean, na.rm = TRUE)
+ if(sd == 0){
+ warning("'mad(x, na.rm = TRUE) = 0' => cannot compute a valid initial estimate.
+ To avoid division by zero 'mad0' is used. You could also specify
+ a valid scale estimate via 'initial.est'.")
+ sd <- mad0
+ }
}else{
if(!is.numeric(initial.est) || length(initial.est) != 1)
stop("'initial.est' needs to be a numeric vector of length 1 or missing")
@@ -596,6 +714,7 @@
if(!missing(eps)){
r <- sqrt(length(x))*eps
+ if(fsCor) r <- finiteSampleCorrection(r = r, n = length(x), model = "sc")
if(r > 10){
b <- sd/(4*qnorm(0.75)*dnorm(qnorm(0.75)))
A <- b^2*(1+r^2)
@@ -607,10 +726,17 @@
}
robEst <- .kstep.sc(x = x, initial.est = sd, A = A, a = a, b = b, mean = mean, k = k)
names(robEst$est) <- "sd"
- Info.matrix <- matrix(c("roblox",
- paste("optimally robust estimate for contamination 'eps' =", round(eps, 3),
- "and 'asMSE'")),
- ncol = 2, dimnames = list(NULL, c("method", "message")))
+ if(fsCor){
+ Info.matrix <- matrix(c("roblox",
+ paste("finite-sample corrected optimally robust estimate for contamination 'eps' =", round(eps, 3),
+ "and 'asMSE'")),
+ ncol = 2, dimnames = list(NULL, c("method", "message")))
+ }else{
+ Info.matrix <- matrix(c("roblox",
+ paste("optimally robust estimate for contamination 'eps' =", round(eps, 3),
+ "and 'asMSE'")),
+ ncol = 2, dimnames = list(NULL, c("method", "message")))
+ }
if(returnIC){
w <- new("HampelWeight")
clip(w) <- robEst$b
@@ -658,6 +784,7 @@
info = c("roblox", "optimally robust IC for AL estimators and 'asMSE'"),
w = w, biastype = symmetricBias(), normtype = NormType(),
modifyIC = modIC))
+ Infos(IC1) <- Info.matrix
return(new("kStepEstimate", name = "Optimally robust estimate",
estimate.call = es.call, estimate = robEst$est,
samplesize = length(x), asvar = as.matrix(robEst$A-r^2*robEst$b^2),
@@ -677,6 +804,10 @@
r <- uniroot(.getsInterval, lower = rlo+1e-8, upper = rup,
tol = .Machine$double.eps^0.25, rlo = rlo, rup = rup)$root
}
+ if(fsCor){
+ r.as <- r
+ r <- finiteSampleCorrection(r = r, n = length(x), model = "sc")
+ }
if(r > 10){
b <- sd/(4*qnorm(0.75)*dnorm(qnorm(0.75)))
A <- b^2*(1+r^2)
@@ -698,12 +829,21 @@
}
robEst <- .kstep.sc(x = x, initial.est = sd, A = A, a = a, b = b, mean = mean, k = k)
names(robEst$est) <- "sd"
- Info.matrix <- matrix(c(rep("roblox", 3),
- paste("radius-minimax estimate for contamination interval [",
- round(eps.lower, 3), ", ", round(eps.upper, 3), "]", sep = ""),
- paste("least favorable contamination: ", round(r/sqrtn, 3), sep = ""),
- paste("maximum MSE-inefficiency: ", round(ineff, 3), sep = "")),
- ncol = 2, dimnames = list(NULL, c("method", "message")))
+ if(fsCor){
+ Info.matrix <- matrix(c(rep("roblox", 3),
+ paste("finite-sample corrected radius-minimax estimate for contamination interval [",
+ round(eps.lower, 3), ", ", round(eps.upper, 3), "]", sep = ""),
+ paste("least favorable (uncorrected) contamination: ", round(r.as/sqrtn, 3), sep = ""),
+ paste("maximum asymptotic MSE-inefficiency: ", round(ineff, 3), sep = "")),
+ ncol = 2, dimnames = list(NULL, c("method", "message")))
+ }else{
+ Info.matrix <- matrix(c(rep("roblox", 3),
+ paste("radius-minimax estimate for contamination interval [",
+ round(eps.lower, 3), ", ", round(eps.upper, 3), "]", sep = ""),
+ paste("least favorable contamination: ", round(r/sqrtn, 3), sep = ""),
+ paste("maximum MSE-inefficiency: ", round(ineff, 3), sep = "")),
+ ncol = 2, dimnames = list(NULL, c("method", "message")))
+ }
if(returnIC){
w <- new("HampelWeight")
clip(w) <- robEst$b
@@ -751,12 +891,7 @@
info = c("roblox", "optimally robust IC for AL estimators and 'asMSE'"),
w = w, biastype = symmetricBias(), normtype = NormType(),
modifyIC = modIC))
- Infos(IC1) <- matrix(c(rep("roblox", 3),
- paste("radius-minimax IC for contamination interval [",
- round(eps.lower, 3), ", ", round(eps.upper, 3), "]", sep = ""),
- paste("least favorable contamination: ", round(r/sqrtn, 3), sep = ""),
- paste("maximum MSE-inefficiency: ", round(ineff, 3), sep = "")),
- ncol = 2, dimnames = list(NULL, c("method", "message")))
+ Infos(IC1) <- Info.matrix
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/robast -r 314
More information about the Robast-commits
mailing list