>>>>> "Morten" == Morten Welinder <[EMAIL PROTECTED]>
>>>>> 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
______________________________________________
[EMAIL PROTECTED] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
