For the record, this is what I am going to use in Gnumeric for now.
Morten
/*
* Problems fixed.
*
* 1. x = -inf, scale = +inf: result now 0, was nan.
*
* 2. No discontinuity for very small alph.
* Old pgamma(1,1e-8,1,FALSE,TRUE) --> -19.9376 [ok]
* Old pgamma(1,0.9999999999e-8,1,FALSE,TRUE) --> -5556027.755 [bad]
* That's more than <dr_evil>two million</dr_evil> orders of magnitude caused
* by the pnorm approximation being used in an area it was not intended for.
*
* 3. No discontinuity at alph=100000.
* Old pgamma(1e6,100000,1,FALSE,TRUE) --> -669750.36332851 [ok]
* Old pgamma(1e6,100001,1,FALSE,TRUE) --> -599731.36192347 [bad]
*
* 4. Improved precision in calculating log of series sum using log1p.
*
* 5. Avoid using series and continued fraction too close to the critical
* line.
*
* 6. Improved precision of argument in pnorm approximation.
*
* 7. Use a power of 2 for rescaling in continued fraction so rescaling does
* not introduce any rounding errors.
*/
/*
* Abramowitz and Stegun 6.5.29
*/
static double
pgamma_series_sum (double x, double alph)
{
double sum = 0, c = 1, a = alph;
do {
a += 1.;
c *= x / a;
sum += c;
} while (c > DBL_EPSILON * sum);
return sum;
}
/*
* Abramowitz and Stegun 6.5.31
*/
static double
pgamma_cont_frac (double x, double alph)
{
double a, b, n, an, pn1, pn2, pn3, pn4, pn5, pn6, sum, osum;
double xlarge = pow (2, 128);
a = 1. - alph;
b = a + x + 1.;
pn1 = 1.;
pn2 = x;
pn3 = x + 1.;
pn4 = x * b;
sum = pn3 / pn4;
for (n = 1; ; n++) {
a += 1.;/* = n+1 -alph */
b += 2.;/* = 2(n+1)-alph+x */
an = a * n;
pn5 = b * pn3 - an * pn1;
pn6 = b * pn4 - an * pn2;
if (fabs(pn6) > 0.) {
osum = sum;
sum = pn5 / pn6;
if (fabs(osum - sum) <= DBL_EPSILON * fmin2(1., sum))
break;
}
pn1 = pn3;
pn2 = pn4;
pn3 = pn5;
pn4 = pn6;
if (fabs(pn5) >= xlarge) {
/* re-scale the terms in continued fraction if they are large */
#ifdef DEBUG_p
REprintf(" [r] ");
#endif
pn1 /= xlarge;
pn2 /= xlarge;
pn3 /= xlarge;
pn4 /= xlarge;
}
}
return sum;
}
double pgamma(double x, double alph, double scale, int lower_tail, int log_p)
{
double arg;
#ifdef IEEE_754
if (ISNAN(x) || ISNAN(alph) || ISNAN(scale))
return x + alph + scale;
#endif
if(alph <= 0. || scale <= 0.)
ML_ERR_return_NAN;
if (x <= 0.)
return R_DT_0;
x /= scale;
#ifdef IEEE_754
if (ISNAN(x)) /* eg. original x = scale = +Inf */
return x;
#endif
if (x < 0.99 * (alph + 1000) && x < 10 + alph && x < 10 * alph) {
/*
* The first condition makes sure that terms in the sum get at
* least 1% smaller, no later than after 1000 terms.
*
* The second condition makes sure that the terms start getting
* smaller (even if just a tiny bit) after no more than 10 terms.
*
* The third condition makes sure that the terms don't get huge
* before they start getting smaller.
*/
double C1mls2p = /* 1 - log (sqrt (2 * M_PI)) */
0.081061466795327258219670263594382360138602526362216587182846;
double logsum = log1p (pgamma_series_sum (x, alph));
/*
* The first two terms here would cause cancellation for extremely
* large and near-equal x and alph. The first condition above
* prevents that.
*/
arg = (alph - x) + alph * log (x / (alph + 1)) + C1mls2p -
log1p (alph) / 2 - stirlerr (alph + 1) + logsum;
} else if (alph < x && alph - 100 < 0.99 * x) {
/*
* The first condition guarantees that we will start making
* a tiny bit of progress right now.
*
* The second condition guarantees that we will make decent
* progress after no more than 100 loops.
*/
double lfrac = log (pgamma_cont_frac (x, alph));
arg = lfrac + alph * log(x) - x - lgammafn(alph);
lower_tail = !lower_tail;
} else {
/*
* Approximation using pnorm. We only get here for fairly large
* x and alph that are within 1% of each other.
*/
double r3m1 = expm1 (log (x / alph) / 3);
return pnorm (sqrt (alph) * 3 * (r3m1 + 1 / (9 * alph)),
0, 1, lower_tail, log_p);
}
if (lower_tail)
return log_p ? arg : exp(arg);
else
return log_p ? R_Log1_Exp(arg) : -expm1(arg);
}
______________________________________________
[EMAIL PROTECTED] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
