[RsR] Huber rho function evaluation

[Martin Maechler]

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


------- 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

    > Note that the Huber psi function evaluation (deriv=0) works well

    >> Mpsi(seq(-5,5,by=0.1),psi="huber",cc=1.35,deriv=0)
    >   [1] -1.35 -1.35 -1.35 -1.35 -1.35 -1.35 -1.35 -1.35 -1.35 -1.35 -1.35 
    > -1.35 -1.35 -1.35 -1.35 -1.35 -1.35 -1.35 -1.35 -1.35 -1.35 -1.35 -1.35 
    > -1.35
    >  [25] -1.35 -1.35 -1.35 -1.35 -1.35 -1.35 -1.35 -1.35 -1.35 -1.35 -1.35 
    > -1.35 -1.35 -1.30 -1.20 -1.10 -1.00 -0.90 -0.80 -0.70 -0.60 -0.50 -0.40 
    > -0.30
    >  [49] -0.20 -0.10  0.00  0.10  0.20  0.30  0.40  0.50 0.60  0.70  0.80  
    > 0.90  1.00  1.10  1.20  1.30  1.35  1.35  1.35 1.35  1.35  1.35  1.35  1.35
    >  [73]  1.35  1.35  1.35  1.35  1.35  1.35  1.35  1.35 1.35  1.35  1.35  
    > 1.35  1.35  1.35  1.35  1.35  1.35  1.35  1.35 1.35  1.35  1.35  1.35  1.35
    >  [97]  1.35  1.35  1.35  1.35  1.35

    > and also the rho evaluation for other choices of functions (e.g. Tukey's 
    > biweight)

    >> Mpsi(seq(-5,5,by=0.1),psi="biweight",cc=1.35,deriv=-1)
    >   [1] 0.303750000 0.303750000 0.303750000 0.303750000 0.303750000 
    > 0.303750000 0.303750000 0.303750000 0.303750000 0.303750000 0.303750000 
    > 0.303750000
    >  [13] 0.303750000 0.303750000 0.303750000 0.303750000 0.303750000 
    > 0.303750000 0.303750000 0.303750000 0.303750000 0.303750000 0.303750000 
    > 0.303750000
    >  [25] 0.303750000 0.303750000 0.303750000 0.303750000 0.303750000 
    > 0.303750000 0.303750000 0.303750000 0.303750000 0.303750000 0.303750000 
    > 0.303750000
    >  [37] 0.303750000 0.303633276 0.300941930 0.292219930 0.275829615 
    > 0.251666667 0.220780758 0.185032340 0.146785551 0.108637255 0.073182210 
    > 0.042814358
    >  [49] 0.019564254 0.004972615 0.000000000 0.004972615 0.019564254 
    > 0.042814358 0.073182210 0.108637255 0.146785551 0.185032340 0.220780758 
    > 0.251666667
    >  [61] 0.275829615 0.292219930 0.300941930 0.303633276 0.303750000 
    > 0.303750000 0.303750000 0.303750000 0.303750000 0.303750000 0.303750000 
    > 0.303750000
    >  [73] 0.303750000 0.303750000 0.303750000 0.303750000 0.303750000 
    > 0.303750000 0.303750000 0.303750000 0.303750000 0.303750000 0.303750000 
    > 0.303750000
    >  [85] 0.303750000 0.303750000 0.303750000 0.303750000 0.303750000 
    > 0.303750000 0.303750000 0.303750000 0.303750000 0.303750000 0.303750000 
    > 0.303750000
    >  [97] 0.303750000 0.303750000 0.303750000 0.303750000 0.303750000


    > Thank you for your help,
    > Eva

    > -- 
    > Prof. Eva Cantoni
    > Research Center for Statistics and
    > Geneva School of Economics and Management
    > University of Geneva, Bd du Pont d'Arve 40, CH-1211 Genève 4
    > gsem.unige.ch/rcs/cantoni

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