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

[Ben Bolker]

  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, <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?
> 
> Serguei.
> 
>>   Solution:
>> work with log-probabilities instead.
>>
>> 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
>>
> 
>