From bbolker @ending from gm@il@com  Tue Dec  4 17:23:09 2018
From: bbolker @ending from gm@il@com (Ben Bolker)
Date: Tue, 4 Dec 2018 11:23:09 -0500
Subject: [Rd] Bug report: Function ppois(0:20,
 lambda=0.9) does not generate a non-decreasing result.
In-Reply-To: <2787db5d-a009-40f5-bea6-8ddb10b10a50@insa-toulouse.fr>
References: <f49f51f.a07fd.16774649817.Coremail.qweytr1@mail.ustc.edu.cn>
 <CALEXWq1WFLzQBgm8oF1myry25DZEyEA_jT3kFEnLYP55-06QKg@mail.gmail.com>
 <2787db5d-a009-40f5-bea6-8ddb10b10a50@insa-toulouse.fr>
Message-ID: <8c485115-a156-23f1-12e3-d57876299b02@gmail.com>


  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 at 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 at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-devel
>>
> 
>



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