Re: [R] Exactness of ppois
> "MM" == Martin Maechler <[EMAIL PROTECTED]>
> on Fri, 16 Jan 2004 14:27:03 +0100 writes:
> "Matthias" == Matthias Kohl <[EMAIL PROTECTED]>
> on Thu, 15 Jan 2004 13:55:22 + 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 <[EMAIL PROTECTED]> http://stat.ethz.ch/~maechler/
Seminar fuer Statistik, ETH-Zentrum LEO C16Leonhardstr. 27
ETH (Federal Inst. Technology) 8092 Zurich SWITZERLAND
phone: x-41-1-632-3408 fax: ...-1228 <><
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Exactness of ppois
> "Matthias" == Matthias Kohl <[EMAIL PROTECTED]>
> on Thu, 15 Jan 2004 13:55:22 + 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.
correct. Namely, because
ppois(x, lambda, lower_tail, log_p) :=
pgamma(lambda, x + 1, 1., !lower_tail, log_p)
and pgamma(x, alph, scale) uses AS 239, currently.
So this thread is really about the precision of R's current pgamma().
In your example, (x = 977.8, alph = 1002, scale=1) and
in pgamma.c,
alphlimit = 1000;
and later
/* use a normal approximation if alph > alphlimit */
if (alph > alphlimit) {
pn1 = sqrt(alph) * 3. * (pow(x/alph, 1./3.) + 1. / (9. * alph) - 1.);
return pnorm(pn1, 0., 1., lower_tail, log_p);
}
So, we could conceivably
improve the situation by increasing `alphlimit'.
Though, I don't see a real need for this (and it will cost CPU
cycles in these cases).
Matthias> Do you think this could cause other problems apart
Matthias> from that admittedly extreme example?
no, I don't think. Look at
> lam <- 977.8
> (p1 <- ppois(1001, lam))
[1] 0.77643705
> (p2 <- sum(dpois(0:1001, lam)))
[1] 0.77643187
Can you imagine a situation where this difference matters?
Matthias> Thanks for your attention!
You're welcome.
Martin Maechler <[EMAIL PROTECTED]> http://stat.ethz.ch/~maechler/
Seminar fuer Statistik, ETH-Zentrum LEO C16Leonhardstr. 27
ETH (Federal Inst. Technology) 8092 Zurich SWITZERLAND
phone: x-41-1-632-3408 fax: ...-1228 <><
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Exactness of ppois
Hello,
by checking the precision of a convolution algorithm, we found the
following "inexactness":
We work with R Version 1.8.1 (2003-11-21) on Windows systems (NT, 2000,
XP).
Try the code:
## Kolmogorov distance between two methods to
## determine P(Poisson(lambda)<=x)
Kolm.dist <- function(lam, eps){
x <- seq(0,qpois(1-eps, lambda=lam), by=1)
max(abs(ppois(x, lambda=lam)-cumsum(dpois(x, lambda=lam
}
erg<-optimize(Kolm.dist, lower=900, upper=1000, maximum=TRUE, eps=1e-15)
erg
Kolm1.dist <- function(lam, eps){
x <- seq(0,qpois(1-eps, lambda=lam), by=1)
which.max(abs(ppois(x, lambda=lam)-cumsum(dpois(x, lambda=lam
}
Kolm1.dist(lam=erg$max, eps=1e-15)
So for lambda=977.8 and x=1001 we get a distance of about 5.2e-06.
(This inexactness seems to hold for all lambda values greater than about
900.)
BUT, summing about 1000 terms of exactness around 1e-16,
we would expect an error of order 1e-13.
We suspect algorithm AS 239 to cause that flaw.
Do you think this could cause other problems apart from
that admittedly extreme example?
Thanks for your attention!
Matthias
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
