[RsR] Huber rho function evaluation (plus plot)

[Martin Maechler]

>>>>> Martin Maechler 
>>>>>     on Fri, 15 Jun 2018 10:16:24 +0200 writes:

>>>>> Eva Cantoni 
>>>>>     on Thu, 14 Jun 2018 14:56:14 +0200 writes:

    >> Dear all, I would like to use the function Mpsi() to
    >> compute Huber's rho function (deriv=-1), but I
    >> consistently get "Inf" values (and even NaN). Am I
    >> missing something?

    >> -----
    >>> Mpsi(seq(-5,5,by=0.1),psi="huber",cc=1.35,deriv=-1)
    >>   [1] Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf
    >> Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf
    >> Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf  [38] Inf Inf Inf
    >> Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf NaN Inf Inf Inf
    >> Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf
    >> Inf Inf Inf Inf Inf Inf  [75] Inf Inf Inf Inf Inf Inf Inf
    >> Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf
    >> Inf Inf Inf Inf Inf Inf
    >> ----

    > Dear Eva,

    > The above is a bug.

    > I'm giving you (and more importantly everyone else reading
    > this) some context, and then a workaround (= "solution for
    > now").  As an appetizer for the latter, use

    >     plot(huberPsi)

    > to get the nice plot appended, showing that robustbase
    > *can* compute rho() {and much more} nicely.

and I forgot the plot.  Appended now (hopefully).

-------------- next part --------------
A non-text attachment was scrubbed...
Name: huberPsi.png
Type: image/png
Size: 23005 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-robust/attachments/20180615/697a7654/attachment.png>

-------------- next part --------------



    > ------- Context --------- As you know, in the newer
    > robustness literature, the definition of rho() is
    > unfortunately available in two different versions, one of
    > them, for redescending psi() only -- scaled to have
    > rho~(Inf) = 1.

    > In robustbase, we call that version chi() or rho~()
    > ["tilde"], and use rho() for the "principal function" of
    > psi(), i.e., psi() = rho'().

    > The Mpsi(), Mwgt(), Mchi() functions were implemented in
    > C, notably to be used by lmrob(), but of course also with
    > an R level interface that you are using.

    > Unfortunately (for this case), the C code was heavily
    > geared toward redescending psi() functions, and so,
    > internally, rho() is defined as "un-scaling" the rho~() =
    > chi() :

    >    rho(x) := rho~(x) * rho(Inf)

    > and of course this is nonsense for Huber (and other
    > non-redescenders) and gives the Inf's you have seen [apart
    > from rho(0) <- rho~(x) * rho(Inf) = 0 * Inf = NaN ]

    > ------- Workaround ---------

    > 1. Use Mchi(*, "huber", deriv=0 ) instead of Mpsi(*,
    > "huber", deriv=-1)

    >           (Uses the fast internal C code) :

    >> > Mchi(seq(-3,3, by=1/4), psi="huber", cc=1.35)
    >  [1] 3.13875 2.80125 2.46375 2.12625 1.78875 1.45125
    > 1.11375 0.78125 [9] 0.50000 0.28125 0.12500 0.03125
    > 0.00000 0.03125 0.12500 0.28125 [17] 0.50000 0.78125
    > 1.11375 1.45125 1.78875 2.12625 2.46375 2.80125 [25]
    > 3.13875
    >> 
   
    > 2. Use the psiFunc() based 'huberPsi' object, {mentioned
    > above that `plot(huberPsi)` produced the nice plot in the
    > appendix},
   

    >    huberPsi using rho(x, seq(-3,3, by=1/4)) # is for the default
    > cc = 1.345

    > or for another k than the default :

    >> Hub_1.5 <- chgDefaults(huberPsi, k = 1.5)
    >> Hub_1.5 using rho(seq(-3,3, by=1/4))
    >  [1] 3.37500 3.00000 2.62500 2.25000 1.87500 1.50000
    > 1.12500 0.78125 [9] 0.50000 0.28125 0.12500 0.03125
    > 0.00000 0.03125 0.12500 0.28125 [17] 0.50000 0.78125
    > 1.12500 1.50000 1.87500 2.25000 2.62500 3.00000 [25]
    > 3.37500
    >> 

    > --------------

    > Of course, I will fix the problem so that Mpsi(x,
    > psi="huber", 1.35, deriv=-1)

    > will also work in the next version of robustbase.

    > Best regards, Martin