[Rd] pgamma discontinuity (PR#7307)

[maechler at stat.math.ethz.ch]

>>>>> "Morten" == Morten Welinder <terra at gnome.org>
>>>>>     on Mon, 25 Oct 2004 12:04:08 -0400 (EDT) writes:

    Morten> A little code study, formula study and experimentation reveals that the
    Morten> situation is mostly fixable:

    Morten> 1. Get rid of the explicit alpha limit.  (A form of
    Morten> it is implicit in (2) and (3) below.)

    Morten> 2. Use the series formula when

    Morten> (x < alph + 10 && x < 0.99 * (alph + 1000))

    Morten> This guarantees that the sum converges reasonably
    Morten> fast.  (For extremely large x and alpha that will
    Morten> take about -53/log2(0.99) iterations for 53
    Morten> significant bits, i.e., about 3700 iterations.)

    Morten> 3. Use the continued fraction formula when

    Morten> (alph < x && alph - 100 < 0.99 * x)

    Morten> Aka, you don't want to use the formula either near
    Morten> the critical point where alpha/x ~ 1 unless the
    Morten> numbers are small.

    Morten> 4a. Go to a library and figure out how Temme does it
    Morten> for alpha near x, both large.  In this case the 0.99
    Morten> from above could probably be lowered a lot for
    Morten> faster convergence.

    Morten> or

    Morten> 4b. Use the pnorm approximation.  It seems to do a
    Morten> whole lot better for alpha near x than it did for
    Morten> the 10:1 case I quoted.

    Morten> Comments, please.

Hi Morten,
thanks a lot for your investigation.

I have spent quite a few working days exploring pgamma() and
also some alternatives.  The discontinuity is "well know".  I
vaguely remember my conclusions were a bit different - at least
over all problems (not only the one mentioned), it needed more
fixing. I think I had included Temme's paper (and others) in my study.

But really, I'm just talking from the top of my head; I will
take the time to look into my old testing scripts and
alternative C code; but just not this week (release of
bioconductor upcoming; other urgencies).

I'll be glad to also communicate with you offline on this topic
{and about pbeta() !}.  But "just not now".

Martin Maechler, ETH Zurich