From: <maechler_at_stat.math.ethz.ch>

Date: Tue 26 Oct 2004 - 04:19:18 EST

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

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.99Morten> 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

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
*