R versions 2.1.0 to present.

Examples shown were computed under Windows R-devel, current SVN, but ix86

Linux shows similar behaviour (sometimes NaN or -Inf rather than Inf,

depending on the compiler and optimization level used).

The replacement pgamma algorithm used from R 2.1.0 has an inadequate

design and no supporting documentation whatsoever. There is no reference

given to support the algorithm, and it seems very desirable to use only

algorithms with a published (and preferably refereed) analysis, or at

least of impeccable provenance.

The following errors were found by investigating an example in the

d-p-q-r-tests.R regression tests that gave NaN on a real system.

These errors were not present in R 2.0.0, which used a normal

approximation in that region. We could fix this by reverting where needed

to a normal approximation, but that leaves the problem that we have no

proof of the validity or accuracy of the rest of the algorithm (if indeed

it is accurate).

?pgamma says

As from R 2.1.0 'pgamma()' uses a new algorithm (mainly by Morten

Welinder) which should be uniformly as accurate as AS 239.

Well, it 'should be' but it is not, and we should not be making statements

like that. Those in the email quoted in pgamma.c exhibit hubris.

There are also at least two examples of sloppy coding that lead to numeric

overflow and complete loss of accuracy.

Consider

> pgamma(seq(0.75, 1.25, by=0.05)*1e100, shape = 1e100, log=T)

[1] -3.768207e+98 NaN NaN NaN NaN

[6] -6.931472e-01 NaN NaN NaN NaN

[11] 0.000000e+00

Warning message:

NaNs produced in: pgamma(q, shape, scale, lower.tail, log.p)

> pgamma(seq(0.75, 1.25, by=0.05)*1e100, shape = 1e100, log=T, lower=F)

[1] 0.000000e+00 NaN NaN NaN NaN

[6] -6.931472e-01 Inf Inf Inf Inf

[11] -2.685645e+98

Warning message:

NaNs produced in: pgamma(q, shape, scale, lower.tail, log.p)

> pgamma(c(1-1e-10, 1+1e-10)*1e100, shape = 1e100)

[1] NaN NaN

(shape=1e25 is enough to cause a breakdown in the first of these, and

1e60 in the rest.)

The code has four branches

1) x <= 1

2) x <= alph - 1 && x < 0.8 * (alph + 50)). This has the comment

/* incl. large alph */, but that is false.

3) if (alph - 1 < x && alph < 0.8 * (x + 50)). This has the comment

/* incl. large x */, but again false.

4) The rest, which uses an asymptotic expansion in

pt_ = -x * (log(1 + (lambda-x)/x) - (lambda-x)/x)

= -x * log((lambda-x)/x) - (lambda-x)

and naively assumes that this is small enough to use a power series

expansion in 1/x with coefficients as powers of pt_. To make matters

worse, consider

> pgamma(0.9*1e25, 1e25, log=T)

pgamma_raw(x=9e+024, alph=1e+025, low=1, log=1)

using ppois_asymp()

pt_ = 5.3605156578263e+022

pp*_asymp(): f=-2.0803930924076e-013 np=-5.3605156578263e+022

nd=-5.3605156578263e+022 f*nd=11151979746.284

[1] -Inf

Hmm, how did that manage to lose *all* the accuracy here? Hubris again

appears in the comments.

Here np and nd are on log scale and if they are large they will be almost

equal (and negative), and f is not large (and if it were we could have

computed log f). So we can compute the log of np+f*nd accurately as

log(np*(1+f*nd/np)) = lnp + log(1+f*nd/np) = lnp + log1p(f*exp(lnd-lnp))

Almost all the mass of gamma(shape=1e100) is concentrated at the nearest

representable value to 1e100:

> qgamma(c(1e-16, 1-1e-16), 1e100)-1e100

[1] 0 0

(if it can be believed, but this can be verified independently). So being

accurate in the middle of the range is pretty academic, but one can at

least avoid returning the nonsense of non-monotone cdfs and NaN/infinite

probabilities.

--

Brian D. Ripley,

[hidden email]
Professor of Applied Statistics,

http://www.stats.ox.ac.uk/~ripley/University of Oxford, Tel: +44 1865 272861 (self)

1 South Parks Road, +44 1865 272866 (PA)

Oxford OX1 3TG, UK Fax: +44 1865 272595

______________________________________________

[hidden email] mailing list

https://stat.ethz.ch/mailman/listinfo/r-devel