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, <qwey...@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@r-project.org mailing list >> https://stat.ethz.ch/mailman/listinfo/r-devel >> > > ______________________________________________ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
