[R] Exactness of ppois

[Martin Maechler]

>>>>> "MM" == Martin Maechler <maechler at stat.math.ethz.ch>
>>>>>     on Fri, 16 Jan 2004 14:27:03 +0100 writes:

>>>>> "Matthias" == Matthias Kohl <Matthias.Kohl at uni-bayreuth.de>
>>>>>     on Thu, 15 Jan 2004 13:55:22 +0000 writes:

    Matthias> Hello, by checking the precision of a convolution
    Matthias> algorithm, we found the following "inexactness":
    Matthias> We work with R Version 1.8.1 (2003-11-21) on
    Matthias> Windows systems (NT, 2000, XP).

    Matthias> Try the code:

    Matthias> So for lambda=977.8 and x=1001 we get a distance
    Matthias> of about 5.2e-06.  (This inexactness seems to hold
    Matthias> for all lambda values greater than about 900.)

    Matthias> BUT, summing about 1000 terms of exactness around 1e-16,
    Matthias> we would expect an error of order 1e-13.

    Matthias> We suspect algorithm AS 239 to cause that flaw.

    MM> correct.   Namely, because

    MM> ppois(x, lambda, lower_tail, log_p) :=  
    MM> pgamma(lambda, x + 1, 1., !lower_tail, log_p)

    MM> and pgamma(x, alph, scale) uses AS 239, currently. 
    MM> So this thread is really about the precision of R's current pgamma().

    MM> In your example, (x = 977.8, alph = 1002, scale=1) and 
    MM> in pgamma.c,
    MM>   alphlimit = 1000;
    MM> and later

    MM> /* use a normal approximation if alph > alphlimit */
    MM>  if (alph > alphlimit) {
    MM>    pn1 = sqrt(alph) * 3. * (pow(x/alph, 1./3.) + 1. / (9. * alph) - 1.);
    MM>    return pnorm(pn1, 0., 1., lower_tail, log_p);
    MM>  }

    MM> So, we could conceivably 
    MM> improve the situation by increasing `alphlimit'.
   
    MM> Though, I don't see a real need for this (and it will cost CPU
    MM> cycles in these cases).

I've now investigated this in some detail.
As a consequence, I will substantially increas 'alphlimit' for
R-devel aka 1.9.0-to-bet,
from 1000. to probably 100'000.

    Matthias> Do you think this could cause other problems apart
    Matthias> from that admittedly extreme example?

    MM> no, I don't think.  Look at

    >> lam <- 977.8
    >> (p1 <- ppois(1001, lam))
    MM> [1] 0.77643705
    >> (p2 <- sum(dpois(0:1001, lam)))
    MM> [1] 0.77643187

    MM> Can you imagine a situation where this difference matters?

Martin Maechler <maechler at stat.math.ethz.ch>	http://stat.ethz.ch/~maechler/
Seminar fuer Statistik, ETH-Zentrum  LEO C16	Leonhardstr. 27
ETH (Federal Inst. Technology)	8092 Zurich	SWITZERLAND
phone: x-41-1-632-3408		fax: ...-1228			<><