From: Martin Maechler <maechler_at_stat.math.ethz.ch>

Date: Thu 22 Dec 2005 - 21:12:02 EST

* >>
*

* >>
*

>> That's Cauchy, BTW, a very extreme case.

R-help@stat.math.ethz.ch mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Thu Dec 22 21:16:39 2005

Date: Thu 22 Dec 2005 - 21:12:02 EST

>>>>> "Murray" == Murray Jorgensen <maj@stats.waikato.ac.nz> >>>>> on Thu, 22 Dec 2005 22:13:45 +1300 writes:

Murray> Prof Brian Ripley wrote:

>> On Thu, 22 Dec 2005, Murray Jorgensen wrote:

* >>
*

>>> We have a choice when calculating the Huber location estimate: >>> > set.seed(221205) >>> > y <- 7 + 3*rt(30,1)

>> That's Cauchy, BTW, a very extreme case.

Murray> Sure, the sort of situation where one might want a robust estimator.

Murray> [...]

>> Note that the huber() scale estimate is the initial MAD, whereas rlm

* >> iterates. Because during iteration the 'center' for MAD is known to be
** >> zero, the results differ. For symmetric distributions there is little
** >> difference, but your sample is not close to symmetric.
*

Murray> [Blush] yes I knew that and somehow I forgot. But leave rlm() alone for Murray> a while and do IRLS with fixed scale:

Murray> th <- median(y) Murray> s <- mad(y) Murray> # paste this in a few times: Murray> w <- ifelse((y-th< 1.345*s & y-th>-1.345*s), 1, 1.345*s/abs(y-th)) Murray> th <- weighted.mean(y,w) Murray> th

Murray> We converge to

* >> th
*

Murray> [1] 5.9203

Murray> close to the answer given by rlm() different from

>> huber(y)$mu

Murray> [1] 5.9117

Murray> So the variable scale does not account for the difference.

No, the main difference is the different default:

huber() has k=1.5

and rlm() has k=1.345 :

Try this

set.seed(221205)

y <- 7 + 3*rt(30,1)

str(huber(y, k = 1.345), digits = 5)

## List of 2 ## $ mu: num 5.9203 ## $ s : num 4.0967

str(rlm(y ~ 1)[c("coefficients", "s")], digits = 5) #

## (edited to) ## $ coefficients: num 5.9204 ## $ s : num 3.7463

which gives 'mu' very close, even for the iterated vs. non-iterated scales.

Martin Maechler, ETH Zurich

R-help@stat.math.ethz.ch mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Thu Dec 22 21:16:39 2005

*
This archive was generated by hypermail 2.1.8
: Fri 03 Mar 2006 - 03:41:40 EST
*