DescTools (version 0.99.8.1)

HuberM: 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.

Value

  • list of location and scale parameters, and number of iterations used.
  • mulocation estimate
  • sthe s argument, typically the mad.
  • itthe number of Huber iterations used.

concept

robust location

Details

Note that currently, when non-NULL weights 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 $\tau$ correction factor but no finite sample correction.

References

Huber, P. J. (1981) Robust Statistics. Wiley.

See Also

hubers (and huber) in package MASS; mad.

Examples

Run this code
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",]), tol= 0.002))

Run the code above in your browser using DataLab