From: <terra_at_gnome.org>

Date: Wed 27 Oct 2004 - 16:15:38 EST

}

}

#endif

R-devel@stat.math.ethz.ch mailing list

https://stat.ethz.ch/mailman/listinfo/r-devel Received on Wed Oct 27 16:23:54 2004

Date: Wed 27 Oct 2004 - 16:15:38 EST

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);}

R-devel@stat.math.ethz.ch mailing list

https://stat.ethz.ch/mailman/listinfo/r-devel Received on Wed Oct 27 16:23:54 2004

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