pgamma - inadequate algorithm design and poor coding (PR#8528)

classic Classic list List threaded Threaded
1 message Options
Reply | Threaded
Open this post in threaded view
|

pgamma - inadequate algorithm design and poor coding (PR#8528)

Prof Brian Ripley
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