R: Safe (generalized) Huber M-Estimator of Location
huberM
R Documentation
Safe (generalized) Huber M-Estimator of Location
Description
(Generalized) Huber M-estimator of location with MAD scale, being
sensible also when the scale is zero where huber()
returns an error.
Usage
huberM(x, k = 1.5, weights = NULL, tol = 1e-06,
mu = if(is.null(weights)) median(x) else wgt.himedian(x, weights),
s = if(is.null(weights)) mad(x, center=mu)
else wgt.himedian(abs(x - mu), weights),
se = FALSE,
warn0scale = getOption("verbose"))
Arguments
x
numeric vector.
k
positive factor; the algorithm winsorizes at k
standard deviations.
weights
numeric vector of non-negative weights of same length
as x, or NULL.
tol
convergence tolerance.
mu
initial location estimator.
s
scale estimator held constant through the iterations.
se
logical indicating if the standard error should be computed
and returned (as SE component). Currently only available
when weights is NULL.
warn0scale
logical; if true, and s is 0 and
length(x) > 1, this will be warned about.
Details
Note that currently, when non-NULLweights are
specified, the default for initial location mu and scale
s is wgt.himedian, where strictly speaking a
weighted “non-hi” median should be used for consistency.
Since s is not updated, the results slightly differ, see the
examples below.
When se = TRUE, the standard error is computed using the
τ correction factor but no finite sample correction.
Value
list of location and scale parameters, and number of iterations used.
mu
location estimate
s
the s argument, typically the mad.
it
the number of “Huber iterations” used.
Author(s)
Martin Maechler, building on the MASS code mentioned.
References
Huber, P. J. (1981)
Robust Statistics.
Wiley.
See Also
hubers (and huber) in package MASS;
mad.
Examples
huberM(c(1:9, 1000))
mad (c(1:9, 1000))
mad (rep(9, 100))
huberM(rep(9, 100))
## When you have "binned" aka replicated observations:
set.seed(7)
x <- c(round(rnorm(1000),1), round(rnorm(50, m=10, sd = 10)))
t.x <- table(x) # -> unique values and multiplicities
x.uniq <- as.numeric(names(t.x)) ## == sort(unique(x))
x.mult <- unname(t.x)
str(Hx <- huberM(x.uniq, weights = x.mult), digits = 7)
str(Hx. <- huberM(x, s = Hx$s, se=TRUE), digits = 7) ## should be ~= Hx
stopifnot(all.equal(Hx[-4], Hx.[-4]))
str(Hx2 <- huberM(x, se=TRUE), digits = 7)## somewhat different, since 's' differs
## Confirm correctness of std.error :
system.time(
SS <- replicate(10000, vapply(huberM(rnorm(400), se=TRUE), as.double, 1.))
) # ~ 12.2 seconds
rbind(mean(SS["SE",]), sd(SS["mu",]))# both ~ 0.0508
stopifnot(all.equal(mean(SS["SE",]),
sd ( SS["mu",]), tolerance= 0.002))