[Robast-commits] r283 - in pkg/RobLox: . R inst man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Apr 18 11:06:33 CEST 2009
Author: stamats
Date: 2009-04-18 11:06:32 +0200 (Sat, 18 Apr 2009)
New Revision: 283
Added:
pkg/RobLox/R/finiteSampleCorrection.R
pkg/RobLox/man/0RobLox-package.Rd
pkg/RobLox/man/finiteSampleCorrection.Rd
Modified:
pkg/RobLox/DESCRIPTION
pkg/RobLox/NAMESPACE
pkg/RobLox/R/colRoblox.R
pkg/RobLox/R/roblox.R
pkg/RobLox/R/rowRoblox.R
pkg/RobLox/R/sysdata.rda
pkg/RobLox/inst/CITATION
pkg/RobLox/inst/NEWS
pkg/RobLox/man/roblox.Rd
pkg/RobLox/man/rowRoblox.Rd
Log:
+ introduction of finite-sample correction
+ handle marginal cases n = 1 or 2 as well as eps = 0 and mad = 0
+ package Rd-file added
Modified: pkg/RobLox/DESCRIPTION
===================================================================
--- pkg/RobLox/DESCRIPTION 2009-04-14 18:20:40 UTC (rev 282)
+++ pkg/RobLox/DESCRIPTION 2009-04-18 09:06:32 UTC (rev 283)
@@ -1,9 +1,9 @@
Package: RobLox
Version: 0.7
-Date: 2009-04-14
-Title: Optimally robust influence curves for location and scale
+Date: 2009-04-18
+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: pkg/RobLox/NAMESPACE
===================================================================
--- pkg/RobLox/NAMESPACE 2009-04-14 18:20:40 UTC (rev 282)
+++ pkg/RobLox/NAMESPACE 2009-04-18 09:06:32 UTC (rev 283)
@@ -1,4 +1,5 @@
-export(rlsOptIC.AL,
+export(finiteSampleCorrection,
+ rlsOptIC.AL,
rlsOptIC.M,
rlsOptIC.BM,
rlsOptIC.Hu1,
Modified: pkg/RobLox/R/colRoblox.R
===================================================================
--- pkg/RobLox/R/colRoblox.R 2009-04-14 18:20:40 UTC (rev 282)
+++ pkg/RobLox/R/colRoblox.R 2009-04-18 09:06:32 UTC (rev 283)
@@ -2,7 +2,7 @@
## Evaluate roblox on columns of a matrix
###############################################################################
colRoblox <- function(x, mean, sd, eps, eps.lower, eps.upper, initial.est,
- k = 1, finiteSampleCorrection = TRUE){
+ k = 1L, fsCor = TRUE, mad0 = 1e-4){
call.est <- match.call()
if(missing(x))
stop("'x' is missing with no default")
@@ -16,7 +16,7 @@
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,
- finiteSampleCorrection = finiteSampleCorrection)
+ fsCor = fsCor, mad0 = mad0)
res at estimate.call <- call.est
return(res)
}
Added: pkg/RobLox/R/finiteSampleCorrection.R
===================================================================
--- pkg/RobLox/R/finiteSampleCorrection.R (rev 0)
+++ pkg/RobLox/R/finiteSampleCorrection.R 2009-04-18 09:06:32 UTC (rev 283)
@@ -0,0 +1,26 @@
+###############################################################################
+## Function for finite-sample correction of the neighborhood radius
+###############################################################################
+finiteSampleCorrection <- function(r, n, model = "locsc"){
+ if(r >= 1.74) 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: pkg/RobLox/R/roblox.R
===================================================================
--- pkg/RobLox/R/roblox.R 2009-04-14 18:20:40 UTC (rev 282)
+++ pkg/RobLox/R/roblox.R 2009-04-18 09:06:32 UTC (rev 283)
@@ -160,27 +160,13 @@
return(list(est = est, A1 = A1, A2 = A2, a = a, b = b, asvar = asVar))
}
-.finiteSampleCorrection.locsc <- function(r, n){
- if(r >= 1.74) return(r)
- 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))
- }
- return(approx(x = epss, y = .finiteSampleRadius.locsc[,ind], xout = eps, rule = 2)$y)
-}
-
###############################################################################
## optimally robust estimator for normal location and/or scale
###############################################################################
-roblox <- function(x, mean, sd, eps, eps.lower, eps.upper, initial.est, k = 1,
- finiteSampleCorrection = TRUE, 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")
@@ -196,6 +182,44 @@
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.")
+ 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
@@ -214,10 +238,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)
@@ -230,13 +292,17 @@
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)
@@ -245,7 +311,7 @@
if(!missing(eps)){
r <- sqrt(length(x))*eps
- if(finiteSampleCorrection) r <- .finiteSampleCorrection.locsc(r = r, n = length(x))
+ if(fsCor) r <- finiteSampleCorrection(r = r, n = length(x), model = "locsc")
if(r > 10){
b <- sd*1.618128043
const <- 1.263094656
@@ -262,10 +328,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
@@ -313,6 +386,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,
@@ -320,9 +394,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,
@@ -342,7 +416,10 @@
r <- uniroot(.getlsInterval, lower = rlo+1e-8, upper = rup,
tol = .Machine$double.eps^0.25, rlo = rlo, rup = rup)$root
}
- if(finiteSampleCorrection) r <- .finiteSampleCorrection.locsc(r = r, n = length(x))
+ if(fsCor){
+ r.as <- r
+ r <- finiteSampleCorrection(r = r, n = length(x), model = "locsc")
+ }
if(r > 10){
b <- sd*1.618128043
const <- 1.263094656
@@ -358,24 +435,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
@@ -423,6 +509,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,
@@ -430,15 +517,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,
@@ -465,6 +546,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)
@@ -474,10 +556,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
@@ -505,6 +594,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),
@@ -524,6 +614,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)
@@ -543,12 +637,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
@@ -576,12 +679,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),
@@ -597,10 +695,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")
@@ -611,6 +712,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)
@@ -622,10 +724,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
@@ -673,6 +782,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),
@@ -692,6 +802,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)
@@ -713,12 +827,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
@@ -766,12 +889,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$est,
samplesize = length(x), asvar = as.matrix(robEst$A-r^2*robEst$b^2),
Modified: pkg/RobLox/R/rowRoblox.R
===================================================================
--- pkg/RobLox/R/rowRoblox.R 2009-04-14 18:20:40 UTC (rev 282)
+++ pkg/RobLox/R/rowRoblox.R 2009-04-18 09:06:32 UTC (rev 283)
@@ -72,7 +72,7 @@
## Evaluate roblox on rows of a matrix
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/robast -r 283
More information about the Robast-commits
mailing list