1 /*
2 * Mathlib : A C Library of Special Functions
3 * Copyright (C) 2005-6 Morten Welinder <terra@gnome.org>
4 * Copyright (C) 2005-10 The R Foundation
5 * Copyright (C) 2006-10 The R Core Team
6 *
7 * This program is free software; you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation; either version 2 of the License, or
10 * (at your option) any later version.
11 *
12 * This program is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with this program; if not, a copy is available at
19 * http://www.r-project.org/Licenses/
20 *
21 * SYNOPSIS
22 *
23 * #include <Rmath.h>
24 *
25 * double pgamma (double x, double alph, double scale,
26 * int lower_tail, int log_p)
27 *
28 * double log1pmx (double x)
29 * double lgamma1p (double a)
30 *
31 * double logspace_add (double logx, double logy)
32 * double logspace_sub (double logx, double logy)
33 *
34 *
35 * DESCRIPTION
36 *
37 * This function computes the distribution function for the
38 * gamma distribution with shape parameter alph and scale parameter
39 * scale. This is also known as the incomplete gamma function.
40 * See Abramowitz and Stegun (6.5.1) for example.
41 *
42 * NOTES
43 *
44 * Complete redesign by Morten Welinder, originally for Gnumeric.
45 * Improvements (e.g. "while NEEDED_SCALE") by Martin Maechler
46 *
47 * REFERENCES
48 *
49 */
50
51 #include "nmath.h"
52 #include "dpq.h"
53 /*----------- DEBUGGING -------------
54 * make CFLAGS='-DDEBUG_p -g -I/usr/local/include -I../include'
55 * (cd ~/R/D/r-devel/Linux-inst/src/nmath; gcc -std=gnu99 -I. -I../../src/include -I../../../R/src/include -I/usr/local/include -DDEBUG_p -g -O2 -c ../../../R/src/nmath/pgamma.c -o pgamma.o)
56 */
57
58 /* Scalefactor:= (2^32)^8 = 2^256 = 1.157921e+77 */
59 #define SQR(x) ((x)*(x))
60 static const double scalefactor = SQR(SQR(SQR(4294967296.0)));
61 #undef SQR
62
63 /* If |x| > |k| * M_cutoff, then log[ exp(-x) * k^x ] =~= -x */
64 static const double M_cutoff = M_LN2 * DBL_MAX_EXP / DBL_EPSILON;/*=3.196577e18*/
65
66 /* Continued fraction for calculation of
67 * 1/i + x/(i+d) + x^2/(i+2*d) + x^3/(i+3*d) + ... = sum_{k=0}^Inf x^k/(i+k*d)
68 *
69 * auxilary in log1pmx() and lgamma1p()
70 */
71 static double
logcf(double x,double i,double d,double eps)72 logcf (double x, double i, double d,
73 double eps /* ~ relative tolerance */)
74 {
75 double c1 = 2 * d;
76 double c2 = i + d;
77 double c4 = c2 + d;
78 double a1 = c2;
79 double b1 = i * (c2 - i * x);
80 double b2 = d * d * x;
81 double a2 = c4 * c2 - b2;
82
83 #if 0
84 assert (i > 0);
85 assert (d >= 0);
86 #endif
87
88 b2 = c4 * b1 - i * b2;
89
90 while (fabs(a2 * b1 - a1 * b2) > fabs(eps * b1 * b2)) {
91 double c3 = c2*c2*x;
92 c2 += d;
93 c4 += d;
94 a1 = c4 * a2 - c3 * a1;
95 b1 = c4 * b2 - c3 * b1;
96
97 c3 = c1 * c1 * x;
98 c1 += d;
99 c4 += d;
100 a2 = c4 * a1 - c3 * a2;
101 b2 = c4 * b1 - c3 * b2;
102
103 if (fabs (b2) > scalefactor) {
104 a1 /= scalefactor;
105 b1 /= scalefactor;
106 a2 /= scalefactor;
107 b2 /= scalefactor;
108 } else if (fabs (b2) < 1 / scalefactor) {
109 a1 *= scalefactor;
110 b1 *= scalefactor;
111 a2 *= scalefactor;
112 b2 *= scalefactor;
113 }
114 }
115
116 return a2 / b2;
117 }
118
119 /* Accurate calculation of log(1+x)-x, particularly for small x. */
log1pmx(double x)120 double log1pmx (double x)
121 {
122 static const double minLog1Value = -0.79149064;
123
124 if (x > 1 || x < minLog1Value)
125 return log1p(x) - x;
126 else { /* -.791 <= x <= 1 -- expand in [x/(2+x)]^2 =: y :
127 * log(1+x) - x = x/(2+x) * [ 2 * y * S(y) - x], with
128 * ---------------------------------------------
129 * S(y) = 1/3 + y/5 + y^2/7 + ... = \sum_{k=0}^\infty y^k / (2k + 3)
130 */
131 double r = x / (2 + x), y = r * r;
132 if (fabs(x) < 1e-2) {
133 static const double two = 2;
134 return r * ((((two / 9 * y + two / 7) * y + two / 5) * y +
135 two / 3) * y - x);
136 } else {
137 static const double tol_logcf = 1e-14;
138 return r * (2 * y * logcf (y, 3, 2, tol_logcf) - x);
139 }
140 }
141 }
142
143
144 /* Compute log(gamma(a+1)) accurately also for small a (0 < a < 0.5). */
lgamma1p(double a)145 double lgamma1p (double a)
146 {
147 const double eulers_const = 0.5772156649015328606065120900824024;
148
149 /* coeffs[i] holds (zeta(i+2)-1)/(i+2) , i = 0:(N-1), N = 40 : */
150 const int N = 40;
151 static const double coeffs[40] = {
152 0.3224670334241132182362075833230126e-0,/* = (zeta(2)-1)/2 */
153 0.6735230105319809513324605383715000e-1,/* = (zeta(3)-1)/3 */
154 0.2058080842778454787900092413529198e-1,
155 0.7385551028673985266273097291406834e-2,
156 0.2890510330741523285752988298486755e-2,
157 0.1192753911703260977113935692828109e-2,
158 0.5096695247430424223356548135815582e-3,
159 0.2231547584535793797614188036013401e-3,
160 0.9945751278180853371459589003190170e-4,
161 0.4492623673813314170020750240635786e-4,
162 0.2050721277567069155316650397830591e-4,
163 0.9439488275268395903987425104415055e-5,
164 0.4374866789907487804181793223952411e-5,
165 0.2039215753801366236781900709670839e-5,
166 0.9551412130407419832857179772951265e-6,
167 0.4492469198764566043294290331193655e-6,
168 0.2120718480555466586923135901077628e-6,
169 0.1004322482396809960872083050053344e-6,
170 0.4769810169363980565760193417246730e-7,
171 0.2271109460894316491031998116062124e-7,
172 0.1083865921489695409107491757968159e-7,
173 0.5183475041970046655121248647057669e-8,
174 0.2483674543802478317185008663991718e-8,
175 0.1192140140586091207442548202774640e-8,
176 0.5731367241678862013330194857961011e-9,
177 0.2759522885124233145178149692816341e-9,
178 0.1330476437424448948149715720858008e-9,
179 0.6422964563838100022082448087644648e-10,
180 0.3104424774732227276239215783404066e-10,
181 0.1502138408075414217093301048780668e-10,
182 0.7275974480239079662504549924814047e-11,
183 0.3527742476575915083615072228655483e-11,
184 0.1711991790559617908601084114443031e-11,
185 0.8315385841420284819798357793954418e-12,
186 0.4042200525289440065536008957032895e-12,
187 0.1966475631096616490411045679010286e-12,
188 0.9573630387838555763782200936508615e-13,
189 0.4664076026428374224576492565974577e-13,
190 0.2273736960065972320633279596737272e-13,
191 0.1109139947083452201658320007192334e-13/* = (zeta(40+1)-1)/(40+1) */
192 };
193
194 const double c = 0.2273736845824652515226821577978691e-12;/* zeta(N+2)-1 */
195 const double tol_logcf = 1e-14;
196 double lgam;
197 int i;
198
199 if (fabs (a) >= 0.5)
200 return lgammafn (a + 1);
201
202 /* Abramowitz & Stegun 6.1.33 : for |x| < 2,
203 * <==> log(gamma(1+x)) = -(log(1+x) - x) - gamma*x + x^2 * \sum_{n=0}^\infty c_n (-x)^n
204 * where c_n := (Zeta(n+2) - 1)/(n+2) = coeffs[n]
205 *
206 * Here, another convergence acceleration trick is used to compute
207 * lgam(x) := sum_{n=0..Inf} c_n (-x)^n
208 */
209 lgam = c * logcf(-a / 2, N + 2, 1, tol_logcf);
210 for (i = N - 1; i >= 0; i--)
211 lgam = coeffs[i] - a * lgam;
212
213 return (a * lgam - eulers_const) * a - log1pmx (a);
214 } /* lgamma1p */
215
216
217
218 /*
219 * Compute the log of a sum from logs of terms, i.e.,
220 *
221 * log (exp (logx) + exp (logy))
222 *
223 * without causing overflows and without throwing away large handfuls
224 * of accuracy.
225 */
logspace_add(double logx,double logy)226 double logspace_add (double logx, double logy)
227 {
228 return fmax2 (logx, logy) + log1p (exp (-fabs (logx - logy)));
229 }
230
231
232 /*
233 * Compute the log of a difference from logs of terms, i.e.,
234 *
235 * log (exp (logx) - exp (logy))
236 *
237 * without causing overflows and without throwing away large handfuls
238 * of accuracy.
239 */
logspace_sub(double logx,double logy)240 double logspace_sub (double logx, double logy)
241 {
242 return logx + R_Log1_Exp(logy - logx);
243 }
244
245
246 /* dpois_wrap (x_P_1, lambda, g_log) ==
247 * dpois (x_P_1 - 1, lambda, g_log) := exp(-L) L^k / gamma(k+1) , k := x_P_1 - 1
248 */
249 static double
dpois_wrap(double x_plus_1,double lambda,int give_log)250 dpois_wrap (double x_plus_1, double lambda, int give_log)
251 {
252 #ifdef DEBUG_p
253 REprintf (" dpois_wrap(x+1=%.14g, lambda=%.14g, log=%d)\n",
254 x_plus_1, lambda, give_log);
255 #endif
256 if (!R_FINITE(lambda))
257 return R_D__0;
258 if (x_plus_1 > 1)
259 return dpois_raw (x_plus_1 - 1, lambda, give_log);
260 if (lambda > fabs(x_plus_1 - 1) * M_cutoff)
261 return R_D_exp(-lambda - lgammafn(x_plus_1));
262 else {
263 double d = dpois_raw (x_plus_1, lambda, give_log);
264 #ifdef DEBUG_p
265 REprintf (" -> d=dpois_raw(..)=%.14g\n", d);
266 #endif
267 return give_log
268 ? d + log (x_plus_1 / lambda)
269 : d * (x_plus_1 / lambda);
270 }
271 }
272
273 /*
274 * Abramowitz and Stegun 6.5.29 [right]
275 */
276 static double
pgamma_smallx(double x,double alph,int lower_tail,int log_p)277 pgamma_smallx (double x, double alph, int lower_tail, int log_p)
278 {
279 double sum = 0, c = alph, n = 0, term;
280
281 #ifdef DEBUG_p
282 REprintf (" pg_smallx(x=%.12g, alph=%.12g): ", x, alph);
283 #endif
284
285 /*
286 * Relative to 6.5.29 all terms have been multiplied by alph
287 * and the first, thus being 1, is omitted.
288 */
289
290 do {
291 n++;
292 c *= -x / n;
293 term = c / (alph + n);
294 sum += term;
295 } while (fabs (term) > DBL_EPSILON * fabs (sum));
296
297 #ifdef DEBUG_p
298 REprintf (" %d terms --> conv.sum=%g;", n, sum);
299 #endif
300 if (lower_tail) {
301 double f1 = log_p ? log1p (sum) : 1 + sum;
302 double f2;
303 if (alph > 1) {
304 f2 = dpois_raw (alph, x, log_p);
305 f2 = log_p ? f2 + x : f2 * exp (x);
306 } else if (log_p)
307 f2 = alph * log (x) - lgamma1p (alph);
308 else
309 f2 = pow (x, alph) / exp (lgamma1p (alph));
310 #ifdef DEBUG_p
311 REprintf (" (f1,f2)= (%g,%g)\n", f1,f2);
312 #endif
313 return log_p ? f1 + f2 : f1 * f2;
314 } else {
315 double lf2 = alph * log (x) - lgamma1p (alph);
316 #ifdef DEBUG_p
317 REprintf (" 1:%.14g 2:%.14g\n", alph * log (x), lgamma1p (alph));
318 REprintf (" sum=%.14g log(1+sum)=%.14g lf2=%.14g\n",
319 sum, log1p (sum), lf2);
320 #endif
321 if (log_p)
322 return R_Log1_Exp (log1p (sum) + lf2);
323 else {
324 double f1m1 = sum;
325 double f2m1 = expm1 (lf2);
326 return -(f1m1 + f2m1 + f1m1 * f2m1);
327 }
328 }
329 } /* pgamma_smallx() */
330
331 static double
pd_upper_series(double x,double y,int log_p)332 pd_upper_series (double x, double y, int log_p)
333 {
334 double term = x / y;
335 double sum = term;
336
337 do {
338 y++;
339 term *= x / y;
340 sum += term;
341 } while (term > sum * DBL_EPSILON);
342
343 /* sum = \sum_{n=1}^ oo x^n / (y*(y+1)*...*(y+n-1))
344 * = \sum_{n=0}^ oo x^(n+1) / (y*(y+1)*...*(y+n))
345 * = x/y * (1 + \sum_{n=1}^oo x^n / ((y+1)*...*(y+n)))
346 * ~ x/y + o(x/y) {which happens when alph -> Inf}
347 */
348 return log_p ? log (sum) : sum;
349 }
350
351 /* Continued fraction for calculation of
352 * scaled upper-tail F_{gamma}
353 * ~= (y / d) * [1 + (1-y)/d + O( ((1-y)/d)^2 ) ]
354 */
355 static double
pd_lower_cf(double y,double d)356 pd_lower_cf (double y, double d)
357 {
358 double f= 0.0 /* -Wall */, of, f0;
359 double i, c2, c3, c4, a1, b1, a2, b2;
360
361 #define NEEDED_SCALE \
362 (b2 > scalefactor) { \
363 a1 /= scalefactor; \
364 b1 /= scalefactor; \
365 a2 /= scalefactor; \
366 b2 /= scalefactor; \
367 }
368
369 #define max_it 200000
370
371 #ifdef DEBUG_p
372 REprintf("pd_lower_cf(y=%.14g, d=%.14g)", y, d);
373 #endif
374 if (y == 0) return 0;
375
376 f0 = y/d;
377 /* Needed, e.g. for pgamma(10^c(100,295), shape= 1.1, log=TRUE): */
378 if(fabs(y - 1) < fabs(d) * DBL_EPSILON) { /* includes y < d = Inf */
379 #ifdef DEBUG_p
380 REprintf(" very small 'y' -> returning (y/d)\n");
381 #endif
382 return (f0);
383 }
384
385 if(f0 > 1.) f0 = 1.;
386 c2 = y;
387 c4 = d; /* original (y,d), *not* potentially scaled ones!*/
388
389 a1 = 0; b1 = 1;
390 a2 = y; b2 = d;
391
392 while NEEDED_SCALE
393
394 i = 0; of = -1.; /* far away */
395 while (i < max_it) {
396
397 i++; c2--; c3 = i * c2; c4 += 2;
398 /* c2 = y - i, c3 = i(y - i), c4 = d + 2i, for i odd */
399 a1 = c4 * a2 + c3 * a1;
400 b1 = c4 * b2 + c3 * b1;
401
402 i++; c2--; c3 = i * c2; c4 += 2;
403 /* c2 = y - i, c3 = i(y - i), c4 = d + 2i, for i even */
404 a2 = c4 * a1 + c3 * a2;
405 b2 = c4 * b1 + c3 * b2;
406
407 if NEEDED_SCALE
408
409 if (b2 != 0) {
410 f = a2 / b2;
411 /* convergence check: relative; "absolute" for very small f : */
412 if (fabs (f - of) <= DBL_EPSILON * fmax2(f0, fabs(f))) {
413 #ifdef DEBUG_p
414 REprintf(" %g iter.\n", i);
415 #endif
416 return f;
417 }
418 of = f;
419 }
420 }
421
422 MATHLIB_WARNING(" ** NON-convergence in pgamma()'s pd_lower_cf() f= %g.\n",
423 f);
424 return f;/* should not happen ... */
425 } /* pd_lower_cf() */
426 #undef NEEDED_SCALE
427
428
429 static double
pd_lower_series(double lambda,double y)430 pd_lower_series (double lambda, double y)
431 {
432 double term = 1, sum = 0;
433
434 #ifdef DEBUG_p
435 REprintf("pd_lower_series(lam=%.14g, y=%.14g) ...", lambda, y);
436 #endif
437 while (y >= 1 && term > sum * DBL_EPSILON) {
438 term *= y / lambda;
439 sum += term;
440 y--;
441 }
442 /* sum = \sum_{n=0}^ oo y*(y-1)*...*(y - n) / lambda^(n+1)
443 * = y/lambda * (1 + \sum_{n=1}^Inf (y-1)*...*(y-n) / lambda^n)
444 * ~ y/lambda + o(y/lambda)
445 */
446 #ifdef DEBUG_p
447 REprintf(" done: term=%g, sum=%g, y= %g\n", term, sum, y);
448 #endif
449
450 if (y != floor (y)) {
451 /*
452 * The series does not converge as the terms start getting
453 * bigger (besides flipping sign) for y < -lambda.
454 */
455 double f;
456 #ifdef DEBUG_p
457 REprintf(" y not int: add another term ");
458 #endif
459 /* FIXME: in quite few cases, adding term*f has no effect (f too small)
460 * and is unnecessary e.g. for pgamma(4e12, 121.1) */
461 f = pd_lower_cf (y, lambda + 1 - y);
462 #ifdef DEBUG_p
463 REprintf(" (= %.14g) * term = %.14g to sum %g\n", f, term * f, sum);
464 #endif
465 sum += term * f;
466 }
467
468 return sum;
469 } /* pd_lower_series() */
470
471 /*
472 * Compute the following ratio with higher accuracy that would be had
473 * from doing it directly.
474 *
475 * dnorm (x, 0, 1, FALSE)
476 * ----------------------------------
477 * pnorm (x, 0, 1, lower_tail, FALSE)
478 *
479 * Abramowitz & Stegun 26.2.12
480 */
481 static double
dpnorm(double x,int lower_tail,double lp)482 dpnorm (double x, int lower_tail, double lp)
483 {
484 /*
485 * So as not to repeat a pnorm call, we expect
486 *
487 * lp == pnorm (x, 0, 1, lower_tail, TRUE)
488 *
489 * but use it only in the non-critical case where either x is small
490 * or p==exp(lp) is close to 1.
491 */
492
493 if (x < 0) {
494 x = -x;
495 lower_tail = !lower_tail;
496 }
497
498 if (x > 10 && !lower_tail) {
499 double term = 1 / x;
500 double sum = term;
501 double x2 = x * x;
502 double i = 1;
503
504 do {
505 term *= -i / x2;
506 sum += term;
507 i += 2;
508 } while (fabs (term) > DBL_EPSILON * sum);
509
510 return 1 / sum;
511 } else {
512 double d = dnorm (x, 0., 1., FALSE);
513 return d / exp (lp);
514 }
515 }
516
517 /*
518 * Asymptotic expansion to calculate the probability that Poisson variate
519 * has value <= x.
520 * Various assertions about this are made (without proof) at
521 * http://members.aol.com/iandjmsmith/PoissonApprox.htm
522 */
523 static double
ppois_asymp(double x,double lambda,int lower_tail,int log_p)524 ppois_asymp (double x, double lambda, int lower_tail, int log_p)
525 {
526 static const double coefs_a[8] = {
527 -1e99, /* placeholder used for 1-indexing */
528 2/3.,
529 -4/135.,
530 8/2835.,
531 16/8505.,
532 -8992/12629925.,
533 -334144/492567075.,
534 698752/1477701225.
535 };
536
537 static const double coefs_b[8] = {
538 -1e99, /* placeholder */
539 1/12.,
540 1/288.,
541 -139/51840.,
542 -571/2488320.,
543 163879/209018880.,
544 5246819/75246796800.,
545 -534703531/902961561600.
546 };
547
548 double elfb, elfb_term;
549 double res12, res1_term, res1_ig, res2_term, res2_ig;
550 double dfm, pt_, s2pt, f, np;
551 int i;
552
553 dfm = lambda - x;
554 /* If lambda is large, the distribution is highly concentrated
555 about lambda. So representation error in x or lambda can lead
556 to arbitrarily large values of pt_ and hence divergence of the
557 coefficients of this approximation.
558 */
559 pt_ = - log1pmx (dfm / x);
560 s2pt = sqrt (2 * x * pt_);
561 if (dfm < 0) s2pt = -s2pt;
562
563 res12 = 0;
564 res1_ig = res1_term = sqrt (x);
565 res2_ig = res2_term = s2pt;
566 for (i = 1; i < 8; i++) {
567 res12 += res1_ig * coefs_a[i];
568 res12 += res2_ig * coefs_b[i];
569 res1_term *= pt_ / i ;
570 res2_term *= 2 * pt_ / (2 * i + 1);
571 res1_ig = res1_ig / x + res1_term;
572 res2_ig = res2_ig / x + res2_term;
573 }
574
575 elfb = x;
576 elfb_term = 1;
577 for (i = 1; i < 8; i++) {
578 elfb += elfb_term * coefs_b[i];
579 elfb_term /= x;
580 }
581 if (!lower_tail) elfb = -elfb;
582 #ifdef DEBUG_p
583 REprintf ("res12 = %.14g elfb=%.14g\n", elfb, res12);
584 #endif
585
586 f = res12 / elfb;
587
588 np = pnorm (s2pt, 0.0, 1.0, !lower_tail, log_p);
589
590 if (log_p) {
591 double n_d_over_p = dpnorm (s2pt, !lower_tail, np);
592 #ifdef DEBUG_p
593 REprintf ("pp*_asymp(): f=%.14g np=e^%.14g nd/np=%.14g f*nd/np=%.14g\n",
594 f, np, n_d_over_p, f * n_d_over_p);
595 #endif
596 return np + log1p (f * n_d_over_p);
597 } else {
598 double nd = dnorm (s2pt, 0., 1., log_p);
599
600 #ifdef DEBUG_p
601 REprintf ("pp*_asymp(): f=%.14g np=%.14g nd=%.14g f*nd=%.14g\n",
602 f, np, nd, f * nd);
603 #endif
604 return np + f * nd;
605 }
606 } /* ppois_asymp() */
607
608
pgamma_raw(double x,double alph,int lower_tail,int log_p)609 double pgamma_raw (double x, double alph, int lower_tail, int log_p)
610 {
611 /* Here, assume that (x,alph) are not NA & alph > 0 . */
612
613 double res;
614
615 #ifdef DEBUG_p
616 REprintf("pgamma_raw(x=%.14g, alph=%.14g, low=%d, log=%d)\n",
617 x, alph, lower_tail, log_p);
618 #endif
619 R_P_bounds_01(x, 0., ML_POSINF);
620
621 if (x < 1) {
622 res = pgamma_smallx (x, alph, lower_tail, log_p);
623 } else if (x <= alph - 1 && x < 0.8 * (alph + 50)) {
624 /* incl. large alph compared to x */
625 double sum = pd_upper_series (x, alph, log_p);/* = x/alph + o(x/alph) */
626 double d = dpois_wrap (alph, x, log_p);
627 #ifdef DEBUG_p
628 REprintf(" alph 'large': sum=pd_upper*()= %.12g, d=dpois_w(*)= %.12g\n",
629 sum, d);
630 #endif
631 if (!lower_tail)
632 res = log_p
633 ? R_Log1_Exp (d + sum)
634 : 1 - d * sum;
635 else
636 res = log_p ? sum + d : sum * d;
637 } else if (alph - 1 < x && alph < 0.8 * (x + 50)) {
638 /* incl. large x compared to alph */
639 double sum;
640 double d = dpois_wrap (alph, x, log_p);
641 #ifdef DEBUG_p
642 REprintf(" x 'large': d=dpois_w(*)= %.14g ", d);
643 #endif
644 if (alph < 1) {
645 if (x * DBL_EPSILON > 1 - alph)
646 sum = R_D__1;
647 else {
648 double f = pd_lower_cf (alph, x - (alph - 1)) * x / alph;
649 /* = [alph/(x - alph+1) + o(alph/(x-alph+1))] * x/alph = 1 + o(1) */
650 sum = log_p ? log (f) : f;
651 }
652 } else {
653 sum = pd_lower_series (x, alph - 1);/* = (alph-1)/x + o((alph-1)/x) */
654 sum = log_p ? log1p (sum) : 1 + sum;
655 }
656 #ifdef DEBUG_p
657 REprintf(", sum= %.14g\n", sum);
658 #endif
659 if (!lower_tail)
660 res = log_p ? sum + d : sum * d;
661 else
662 res = log_p
663 ? R_Log1_Exp (d + sum)
664 : 1 - d * sum;
665 } else { /* x >= 1 and x fairly near alph. */
666 #ifdef DEBUG_p
667 REprintf(" using ppois_asymp()\n");
668 #endif
669 res = ppois_asymp (alph - 1, x, !lower_tail, log_p);
670 }
671
672 /*
673 * We lose a fair amount of accuracy to underflow in the cases
674 * where the final result is very close to DBL_MIN. In those
675 * cases, simply redo via log space.
676 */
677 if (!log_p && res < DBL_MIN / DBL_EPSILON) {
678 /* with(.Machine, double.xmin / double.eps) #|-> 1.002084e-292 */
679 #ifdef DEBUG_p
680 REprintf(" very small res=%.14g; -> recompute via log\n", res);
681 #endif
682 return exp (pgamma_raw (x, alph, lower_tail, 1));
683 } else
684 return res;
685 }
686
687
pgamma(double x,double alph,double scale,int lower_tail,int log_p)688 double pgamma(double x, double alph, double scale, int lower_tail, int log_p)
689 {
690 #ifdef IEEE_754
691 if (ISNAN(x) || ISNAN(alph) || ISNAN(scale))
692 return x + alph + scale;
693 #endif
694 if(alph < 0. || scale <= 0.)
695 ML_ERR_return_NAN;
696 x /= scale;
697 #ifdef IEEE_754
698 if (ISNAN(x)) /* eg. original x = scale = +Inf */
699 return x;
700 #endif
701 if(alph == 0.) /* limit case; useful e.g. in pnchisq() */
702 return (x <= 0) ? R_DT_0: R_DT_1; /* <= assert pgamma(0,0) ==> 0 */
703 return pgamma_raw (x, alph, lower_tail, log_p);
704 }
705 /* From: terra@gnome.org (Morten Welinder)
706 * To: R-bugs@biostat.ku.dk
707 * Cc: maechler@stat.math.ethz.ch
708 * Subject: Re: [Rd] pgamma discontinuity (PR#7307)
709 * Date: Tue, 11 Jan 2005 13:57:26 -0500 (EST)
710
711 * this version of pgamma appears to be quite good and certainly a vast
712 * improvement over current R code. (I last looked at 2.0.1) Apart from
713 * type naming, this is what I have been using for Gnumeric 1.4.1.
714
715 * This could be included into R as-is, but you might want to benefit from
716 * making logcf, log1pmx, lgamma1p, and possibly logspace_add/logspace_sub
717 * available to other parts of R.
718
719 * MM: I've not (yet?) taken logcf(), but the other four
720 */
721
722