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

classic Classic list List threaded Threaded
6 messages Options
Reply | Threaded
Open this post in threaded view
|

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

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


Code is tested in both R 3.5.1 and Microsoft R Open 3.5.1.


               _                          
platform       x86_64-w64-mingw32        
arch           x86_64                    
os             mingw32                    
system         x86_64, mingw32            
status                                    
major          3                          
minor          5.1                        
year           2018                      
month          07                        
day            02                        
svn rev        74947                      
language       R                          
version.string R version 3.5.1 (2018-07-02)
nickname       Feather Spray      







PS. if set lambda=1, such bug will not appear.








        [[alternative HTML version deleted]]

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

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

Iñaki Ucar
On Tue, 4 Dec 2018 at 11:12, <[hidden email]> 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. Solution:
work with log-probabilities instead.

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

Iñaki

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

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

Serguei Sokol
Le 04/12/2018 à 11:27, Iñaki Ucar a écrit :

> On Tue, 4 Dec 2018 at 11:12, <[hidden email]> 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?

Serguei.

>   Solution:
> work with log-probabilities instead.
>
> any(diff(ppois(0:40, lambda=0.9, log.p=TRUE))<0)
> #> [1] FALSE
>
> Iñaki
>
> ______________________________________________
> [hidden email] 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: [hidden email]
http://www.lisbp.fr

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

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

bbolker

  I do think it's plausible to expect that we could get *non-decreasing*
results.

I get

any(diff(exp(ppois(0:19, lambda=0.9, log.p=TRUE)))<0)

as FALSE.

But I do get diff(ppois(18:19, lambda=0.9)) < 0.

Looking at the code of ppois, it's just (within C code) calling pgamma
with pgamma(lambda, shape=1+x, scale=1, lower.tail=FALSE):


identical(
  ppois(18:19,0.9),
  pgamma(0.9,shape=1+(18:19),scale=1,lower.tail=FALSE)
)

is TRUE.  So the problem (if we choose to define it as such) would be in
pgamma (upper tail should be a non-decreasing function of the shape
parameter) ... the code is here
https://github.com/wch/r-source/blob/5a156a0865362bb8381dcd69ac335f5174a4f60c/src/nmath/pgamma.c
 if anyone wants to dig in ...




On 2018-12-04 5:46 a.m., Serguei Sokol wrote:

> Le 04/12/2018 à 11:27, Iñaki Ucar a écrit :
>> On Tue, 4 Dec 2018 at 11:12, <[hidden email]> 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?
>
> Serguei.
>
>>   Solution:
>> work with log-probabilities instead.
>>
>> any(diff(ppois(0:40, lambda=0.9, log.p=TRUE))<0)
>> #> [1] FALSE
>>
>> Iñaki
>>
>> ______________________________________________
>> [hidden email] mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-devel
>>
>
>

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

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

Martin Maechler
In reply to this post by Serguei Sokol
>>>>> 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, <[hidden email]> 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
    >>
    >> ______________________________________________
    >> [hidden email] 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: [hidden email]
    > http://www.lisbp.fr

    > ______________________________________________
    > [hidden email] mailing list
    > https://stat.ethz.ch/mailman/listinfo/r-devel

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

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

Serguei Sokol
Le 04/12/2018 à 19:12, Martin Maechler a écrit :

>>>>>> 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, <[hidden email]> 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
> ...
>      > 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.
I have not doubt about it. My "economic way" was related to get ppois()
*non decreasing*, at least, more economic than exp-log.p trick.

Serguei.

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel