Re: [Rd] pgamma discontinuity (PR#7307)

From: <maechler_at_stat.math.ethz.ch>
Date: Tue 26 Oct 2004 - 04:19:18 EST

>>>>> "Morten" == Morten Welinder <terra@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.

```

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

R-devel@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel Received on Tue Oct 26 04:31:38 2004

This archive was generated by hypermail 2.1.8 : Wed 03 Nov 2004 - 22:45:25 EST