[RsR] [R] M-estimator R function question

[Martin Maechler]

>>>>> "Valentin" == Valentin Todorov <valentin.to using gmail.com>
>>>>>     on Tue, 20 Dec 2005 13:20:33 +0100 writes:

    Valentin> Dear Kjell ,
    Valentin> beware of a minor error in the pairwise code that you posted: the
    Valentin> weights in gk.sigmamu() will never become 0, even if abs(x) > c1.

Indeed!

    Valentin> The correction should be:

>>      gk.sigmamu <- function(x, c1 = 4.5, c2 = 3.0, mu.too = FALSE, ...)
>>      {
>> 	n <- length(x)
>> 	medx <- median(x)
>> 	sigma0 <- median(abs(x - medx))
>> 
>> 	# VT::19.12.2005 - should give 0 weights if abs(x - medx) / sigma0 > c1
>> 	#
>> 	#  w <- (x - medx) / sigma0
>> 	#  w <- (1.0 - (w / c1)^2)^2
>> 	#  w[w < 0.0] <- 0.0
>> 
>> 	w <- abs(x - medx) / sigma0
>> 	w <- ifelse(w<=c1,(1.0 - (w / c1)^2)^2,0)
>> 
>> 	# VT::19.12.2005 - END
>> 
>> 	mu <- sum(x * w) / sum(w)
>> 	...... ... ... ..

    Valentin> best regard
    Valentin> valentin

Note however that Kjell was trying to be smart, namely *not*
using ifelse() in such a case -- for efficiency reasons. But he
forgot to do it properly: Here,  pmax()  is much more efficient
than ifelse().  Note that gk.sigmamu() is really the function that is called
many times, so it makes sense to write it as efficiently as possible
--- still keeping it well readable.

I had also applied a few minor cosmetic improvements to Kjell's
original code --- but overlooked the obvious buglet ! ---  
which now leads to 

gk.sigmamu <- function(x, c1 = 4.5, c2 = 3.0, mu.too = FALSE, ...)
{
    n <- length(x)
    medx <- median(x)
    x. <- abs(x - medx)
    c.sigma0 <- c1 * median(x.)
    w <- pmax(0, 1 - (x. / c.sigma0)^2)^2
    mu <- sum(x * w) / sum(w)

    x <- (x - mu) / sigma0
    rho <- x^2
    rho[rho > c2^2] <- c2^2
    ## sigma2 <- sigma0^2 / n * sum(rho)

    ## return
    c(if(mu.too) mu,
      ## sqrt(sigma2) == sqrt( sigma0^2 / n * sum(rho) ) :
      sigma0 * sqrt(sum(rho)/n))
}


Martin Maechler, ETH Zurich


    Valentin> On 12/5/05, Kjell Konis <konis using stats.ox.ac.uk> wrote:
    >> Here is an implementation of the OGK estimator completely in R.  I
    >> haven't touched it for a while and I forget how thoroughly I tested
    >> it so use it with a bit of caution.
    >> 
    >> http://www.stats.ox.ac.uk/~konis/pairwise.q
    >> 
    >> Ricardo, can you be more specific about what you don't like in the
    >> implementation of OGK in S-Plus?
    >> 
    >> Kjell
    >> 
    >> 
    >> 
    >> On Dec 3, 2005, at 9:55 PM, Ricardo Maronna wrote:
    >> 
    >> > One observation about the following.
    >> >
    >> >> One thing we'd be very interested is the OGK estimator of
    >> >> Maronna and Zamar (JASA, 2002).  Unfortunately, there, too, some
    >> >> of the authors sold their code exclusively to Insightful (the
    >> >> S-plus company).
    >> >
    >> >     Insightful has implemented the procedure in some way I ignore
    >> > (actually,
    >> > I don't like the way it works), but unfortunately I got not a
    >> > single dollar
    >> > from that!. Anybody is free to implement the method (which is
    >> > extremely
    >> > simple) in whatever language.
    >> >             Ricardo
    >> >
    >> > _______________________________________________
    >> > R-SIG-Robust using r-project.org mailing list
    >> > https://stat.ethz.ch/mailman/listinfo/r-sig-robust
    >> 
    >> _______________________________________________
    >> R-SIG-Robust using r-project.org mailing list
    >> https://stat.ethz.ch/mailman/listinfo/r-sig-robust
    >> 

    Valentin> _______________________________________________
    Valentin> R-SIG-Robust using r-project.org mailing list
    Valentin> https://stat.ethz.ch/mailman/listinfo/r-sig-robust