[Rd] Bug report: Function ppois(0:20, lambda=0.9) does not generate a non-decreasing result.

[Martin Maechler]

>>>>> Serguei Sokol 
>>>>>     on Tue, 4 Dec 2018 11:46:32 +0100 writes:

    > Le 04/12/2018 à 11:27, Iñaki Ucar a écrit :
    >> On Tue, 4 Dec 2018 at 11:12, <qweytr1 using mail.ustc.edu.cn> wrote:
    >>> function ppois is a function calculate the CDF of Poisson distribution, it should generate a non-decreasing result, but what I got is:
    >>> 
    >>>> any(diff(ppois(0:19,lambda=0.9))<0)
    >>> [1] TRUE
    >>> 
    >>> Actually,
    >>> 
    >>>> ppois(19,lambda=0.9)<ppois(18,lambda=0.9)
    >>> [1] TRUE
    >>> 
    >>> Which could not be TRUE.
    >> This is just another manifestation of
    >> 
    >> 0.1 * 3 > 0.3
    >> #> [1] TRUE
    >> 
    >> This discussion returns to this list from time to time. TLDR; this is
    >> not an R issue, but an unavoidable floating point issue.
    > Well, here the request may be interpreted not as "do it without round 
    > error" which is indeed unavoidable but rather "please cope with rounding 
    > errors in a way that return consistent result for ppois()". You have 
    > indicated one way to do so (I have just added exp() in the row):

    > any(diff(exp(ppois(0:19, lambda=0.9, log.p=TRUE))) < 0)
    > #[1] FALSE

    > But may be there is another, more economic way?

Well, log probabilites *are* very economic for many such p*()
functions.

OTOH, I'm a bit surprised that nobody mentioned the
'lower.tail=FALSE' option which here makes so very much sense,
and is I think slightly more intuitive than the log-probabilities: 
It gives much much much more accurate results for such outermost
right tail probabilities where p*() ~= 1 :

> ppois(15:19, lambda=0.9)
[1] 1 1 1 1 1
> ppois(15:19, lambda=0.9, lower.tail=FALSE)
[1] 3.801404e-15 2.006332e-16 1.000417e-17 4.727147e-19 2.122484e-20
>

... and if you compare with

> ppois(15:19, lambda=0.9, log.p=TRUE)

and notice how similar the numbers are, you may remember that indeed

log(1-p)  ~=  -p   when  |p| ≪ 1

Martin

    > Serguei.

    >> Solution:
    >> work with log-probabilities instead.

or 

    >> 
    >> any(diff(ppois(0:40, lambda=0.9, log.p=TRUE))<0)
    >> #> [1] FALSE
    >> 
    >> Iñaki
    >> 
    >> ______________________________________________
    >> R-devel using r-project.org mailing list
    >> https://stat.ethz.ch/mailman/listinfo/r-devel
    >> 


    > -- 
    > Serguei Sokol
    > Ingenieur de recherche INRA

    > Cellule mathématiques
    > LISBP, INSA/INRA UMR 792, INSA/CNRS UMR 5504
    > 135 Avenue de Rangueil
    > 31077 Toulouse Cedex 04

    > tel: +33 5 62 25 01 27
    > email: sokol using insa-toulouse.fr
    > http://www.lisbp.fr

    > ______________________________________________
    > R-devel using r-project.org mailing list
    > https://stat.ethz.ch/mailman/listinfo/r-devel