1 module rmathd.gamma; 2 3 public import rmathd.common; 4 public import rmathd.normal; 5 //import rmathd.uniform; 6 7 /* 8 ** dmd gamma.d uniform.d normal.d common.d && ./gamma 9 */ 10 11 T dgamma(T: double)(T x, T shape, T scale, int give_log) 12 { 13 T pr; 14 if (isNaN(x) || isNaN(shape) || isNaN(scale)) 15 return x + shape + scale; 16 17 mixin R_D__0!give_log; 18 19 if (shape < 0 || scale <= 0) 20 return T.nan; 21 if (x < 0) 22 return R_D__0; 23 if (shape == 0) /* point mass at 0 */ 24 return (x == 0)? T.infinity : R_D__0; 25 if (x == 0) { 26 if (shape < 1) 27 return T.infinity; 28 if (shape > 1) 29 return R_D__0; 30 return give_log ? -log(scale) : 1 / scale; 31 } 32 33 if (shape < 1) { 34 pr = dpois_raw!T(shape, x/scale, give_log); 35 return give_log ? pr + log(shape/x) : pr*shape/x; 36 } 37 /* else shape >= 1 */ 38 pr = dpois_raw!T(shape - 1, x/scale, give_log); 39 return give_log ? pr - log(scale) : pr/scale; 40 } 41 42 /* Continued fraction for calculation of 43 * 1/i + x/(i+d) + x^2/(i+2*d) + x^3/(i+3*d) + ... = sum_{k=0}^Inf x^k/(i+k*d) 44 * 45 * auxilary in log1pmx() and lgamma1p() 46 */ 47 static T logcf(T)(T x, T i, T d, T eps /* ~ relative tolerance */) 48 { 49 T c1 = 2 * d; 50 T c2 = i + d; 51 T c4 = c2 + d; 52 T a1 = c2; 53 T b1 = i * (c2 - i * x); 54 T b2 = d * d * x; 55 T a2 = c4 * c2 - b2; 56 57 //#if 0 58 // assert(i > 0); 59 // assert(d >= 0); 60 //#endif 61 62 b2 = c4 * b1 - i * b2; 63 64 while (fabs(a2 * b1 - a1 * b2) > fabs(eps * b1 * b2)) { 65 T c3 = c2*c2*x; 66 c2 += d; 67 c4 += d; 68 a1 = c4 * a2 - c3 * a1; 69 b1 = c4 * b2 - c3 * b1; 70 71 c3 = c1 * c1 * x; 72 c1 += d; 73 c4 += d; 74 a2 = c4 * a1 - c3 * a2; 75 b2 = c4 * b1 - c3 * b2; 76 77 if (fabs (b2) > scalefactor) { 78 a1 /= scalefactor; 79 b1 /= scalefactor; 80 a2 /= scalefactor; 81 b2 /= scalefactor; 82 } else if (fabs (b2) < 1 / scalefactor) { 83 a1 *= scalefactor; 84 b1 *= scalefactor; 85 a2 *= scalefactor; 86 b2 *= scalefactor; 87 } 88 } 89 90 return a2 / b2; 91 } 92 93 /* Accurate calculation of log(1+x)-x, particularly for small x. */ 94 T log1pmx(T)(T x) 95 { 96 static const(T) minLog1Value = -0.79149064; 97 98 if (x > 1 || x < minLog1Value) 99 return log1p(x) - x; 100 else { /* -.791 <= x <= 1 -- expand in [x/(2+x)]^2 =: y : 101 * log(1+x) - x = x/(2+x) * [ 2 * y * S(y) - x], with 102 * --------------------------------------------- 103 * S(y) = 1/3 + y/5 + y^2/7 + ... = \sum_{k=0}^\infty y^k / (2k + 3) 104 */ 105 T r = x / (2 + x), y = r * r; 106 if (fabs(x) < 1e-2) { 107 static const(T) two = 2; 108 return r * ((((two / 9 * y + two / 7) * y + two / 5) * y + 109 two / 3) * y - x); 110 } else { 111 static const(T) tol_logcf = 1e-14; 112 return r * (2 * y * logcf!T(y, 3, 2, tol_logcf) - x); 113 } 114 } 115 } 116 117 /* Compute log(gamma(a+1)) accurately also for small a (0 < a < 0.5). */ 118 T lgamma1p(T)(T a) 119 { 120 const(T) eulers_const = 0.5772156649015328606065120900824024; 121 122 /* coeffs[i] holds (zeta(i+2)-1)/(i+2) , i = 0:(N-1), N = 40 : */ 123 const(int) N = 40; 124 static const(T)[40] coeffs = [ 125 0.3224670334241132182362075833230126e-0,/* = (zeta(2)-1)/2 */ 126 0.6735230105319809513324605383715000e-1,/* = (zeta(3)-1)/3 */ 127 0.2058080842778454787900092413529198e-1, 128 0.7385551028673985266273097291406834e-2, 129 0.2890510330741523285752988298486755e-2, 130 0.1192753911703260977113935692828109e-2, 131 0.5096695247430424223356548135815582e-3, 132 0.2231547584535793797614188036013401e-3, 133 0.9945751278180853371459589003190170e-4, 134 0.4492623673813314170020750240635786e-4, 135 0.2050721277567069155316650397830591e-4, 136 0.9439488275268395903987425104415055e-5, 137 0.4374866789907487804181793223952411e-5, 138 0.2039215753801366236781900709670839e-5, 139 0.9551412130407419832857179772951265e-6, 140 0.4492469198764566043294290331193655e-6, 141 0.2120718480555466586923135901077628e-6, 142 0.1004322482396809960872083050053344e-6, 143 0.4769810169363980565760193417246730e-7, 144 0.2271109460894316491031998116062124e-7, 145 0.1083865921489695409107491757968159e-7, 146 0.5183475041970046655121248647057669e-8, 147 0.2483674543802478317185008663991718e-8, 148 0.1192140140586091207442548202774640e-8, 149 0.5731367241678862013330194857961011e-9, 150 0.2759522885124233145178149692816341e-9, 151 0.1330476437424448948149715720858008e-9, 152 0.6422964563838100022082448087644648e-10, 153 0.3104424774732227276239215783404066e-10, 154 0.1502138408075414217093301048780668e-10, 155 0.7275974480239079662504549924814047e-11, 156 0.3527742476575915083615072228655483e-11, 157 0.1711991790559617908601084114443031e-11, 158 0.8315385841420284819798357793954418e-12, 159 0.4042200525289440065536008957032895e-12, 160 0.1966475631096616490411045679010286e-12, 161 0.9573630387838555763782200936508615e-13, 162 0.4664076026428374224576492565974577e-13, 163 0.2273736960065972320633279596737272e-13, 164 0.1109139947083452201658320007192334e-13/* = (zeta(40+1)-1)/(40+1) */ 165 ]; 166 167 const(T) c = 0.2273736845824652515226821577978691e-12;/* zeta(N+2)-1 */ 168 const(T) tol_logcf = 1e-14; 169 T lgam; 170 int i; 171 172 if (fabs (a) >= 0.5) 173 return lgammafn!T(a + 1); 174 175 /* Abramowitz & Stegun 6.1.33 : for |x| < 2, 176 * <==> log(gamma(1+x)) = -(log(1+x) - x) - gamma*x + x^2 * \sum_{n=0}^\infty c_n (-x)^n 177 * where c_n := (Zeta(n+2) - 1)/(n+2) = coeffs[n] 178 * 179 * Here, another convergence acceleration trick is used to compute 180 * lgam(x) := sum_{n=0..Inf} c_n (-x)^n 181 */ 182 lgam = c * logcf!T(-a / 2, N + 2, 1, tol_logcf); 183 for (i = N - 1; i >= 0; i--) 184 lgam = coeffs[i] - a * lgam; 185 186 return (a * lgam - eulers_const) * a - log1pmx!T(a); 187 } /* lgamma1p */ 188 189 190 /* 191 * Compute the log of a sum from logs of terms, i.e., 192 * 193 * log (exp (logx) + exp (logy)) 194 * 195 * without causing overflows and without throwing away large handfuls 196 * of accuracy. 197 */ 198 T logspace_add(T)(T logx, T logy) 199 { 200 return fmax2!T(logx, logy) + log1p (exp (-fabs (logx - logy))); 201 } 202 203 /* 204 * Compute the log of a difference from logs of terms, i.e., 205 * 206 * log (exp (logx) - exp (logy)) 207 * 208 * without causing overflows and without throwing away large handfuls 209 * of accuracy. 210 */ 211 T logspace_sub(T)(T logx, T logy) 212 { 213 return logx + R_Log1_Exp!T(logy - logx); 214 } 215 216 /* 217 * Compute the log of a sum from logs of terms, i.e., 218 * 219 * log (sum_i exp (logx[i]) ) = 220 * log (e^M * sum_i e^(logx[i] - M) ) = 221 * M + log( sum_i e^(logx[i] - M) 222 * 223 * without causing overflows or throwing much accuracy. 224 */ 225 T logspace_sum(T)(const(T)* logx, int n) 226 { 227 if(n == 0) 228 return -T.infinity; // = log( sum(<empty>) ) 229 if(n == 1) 230 return logx[0]; 231 if(n == 2) 232 return logspace_add(logx[0], logx[1]); 233 // else (n >= 3) : 234 int i; 235 // Mx := max_i log(x_i) 236 T Mx = logx[0]; 237 for(i = 1; i < n; i++) 238 if(Mx < logx[i]) Mx = logx[i]; 239 real s = 0.; 240 for(i = 0; i < n; i++) 241 s += exp(logx[i] - Mx); 242 return Mx + cast(T)log(s); 243 } 244 245 /* dpois_wrap (x__1, lambda) := dpois(x__1 - 1, lambda); where 246 * dpois(k, L) := exp(-L) L^k / gamma(k+1) {the usual Poisson probabilities} 247 * 248 * and dpois*(.., give_log = TRUE) := log( dpois*(..) ) 249 */ 250 static T dpois_wrap(T)(T x_plus_1, T lambda, int give_log) 251 { 252 mixin R_D__0!(give_log); 253 mixin M_cutoff!T; 254 if (!isFinite(lambda)) 255 return R_D__0; 256 if (x_plus_1 > 1) 257 return dpois_raw!T(x_plus_1 - 1, lambda, give_log); 258 if (lambda > fabs(x_plus_1 - 1) * M_cutoff) 259 return R_D_exp!T(-lambda - lgammafn!T(x_plus_1), give_log); 260 else { 261 T d = dpois_raw!T(x_plus_1, lambda, give_log); 262 return give_log 263 ? d + log(x_plus_1 / lambda) 264 : d * (x_plus_1 / lambda); 265 } 266 } 267 268 /* 269 * Abramowitz and Stegun 6.5.29 [right] 270 */ 271 static T pgamma_smallx(T)(T x, T alph, int lower_tail, int log_p) 272 { 273 T sum = 0, c = alph, n = 0, term; 274 275 /* 276 * Relative to 6.5.29 all terms have been multiplied by alph 277 * and the first, thus being 1, is omitted. 278 */ 279 280 do { 281 n++; 282 c *= -x / n; 283 term = c / (alph + n); 284 sum += term; 285 } while (fabs(term) > EPSILON!T * fabs(sum)); 286 287 if (lower_tail) { 288 T f1 = log_p ? log1p(sum) : 1 + sum; 289 T f2; 290 if (alph > 1) { 291 f2 = dpois_raw!T(alph, x, log_p); 292 f2 = log_p ? f2 + x : f2 * exp (x); 293 } else if (log_p) 294 f2 = alph * log (x) - lgamma1p!T(alph); 295 else 296 f2 = pow(x, alph) / exp (lgamma1p!T(alph)); 297 298 return log_p ? f1 + f2 : f1 * f2; 299 } else { 300 T lf2 = alph * log (x) - lgamma1p!T(alph); 301 302 if (log_p) 303 return R_Log1_Exp!T(log1p(sum) + lf2); 304 else { 305 T f1m1 = sum; 306 T f2m1 = expm1 (lf2); 307 return -(f1m1 + f2m1 + f1m1 * f2m1); 308 } 309 } 310 } /* pgamma_smallx() */ 311 312 313 static T pd_upper_series(T)(T x, T y, int log_p) 314 { 315 T term = x / y; 316 T sum = term; 317 318 do { 319 y++; 320 term *= x / y; 321 sum += term; 322 } while (term > sum * EPSILON!T); 323 324 /* sum = \sum_{n=1}^ oo x^n / (y*(y+1)*...*(y+n-1)) 325 * = \sum_{n=0}^ oo x^(n+1) / (y*(y+1)*...*(y+n)) 326 * = x/y * (1 + \sum_{n=1}^oo x^n / ((y+1)*...*(y+n))) 327 * ~ x/y + o(x/y) {which happens when alph -> Inf} 328 */ 329 return log_p ? log(sum) : sum; 330 } 331 332 333 immutable(string) NEEDED_SCALE(string ctrl)() 334 { 335 immutable(string) output = ctrl ~ `(b2 > scalefactor) { 336 a1 /= scalefactor; 337 b1 /= scalefactor; 338 a2 /= scalefactor; 339 b2 /= scalefactor; 340 }`; 341 return output; 342 } 343 344 345 static T pd_lower_cf (T)(T y, T d) 346 { 347 T f = 0.0 /* -Wall */, of, f0; 348 T i, c2, c3, c4, a1, b1, a2, b2; 349 350 enum max_it = 200000; 351 352 if (y == 0) 353 return 0; 354 355 f0 = y/d; 356 /* Needed, e.g. for pgamma(10^c(100,295), shape= 1.1, log=TRUE): */ 357 if(fabs(y - 1) < fabs(d) * EPSILON!T) { /* includes y < d = Inf */ 358 return (f0); 359 } 360 361 if(f0 > 1.) f0 = 1.; 362 c2 = y; 363 c4 = d; /* original (y,d), *not* potentially scaled ones!*/ 364 365 a1 = 0; b1 = 1; 366 a2 = y; b2 = d; 367 368 mixin (NEEDED_SCALE!("while")()); 369 370 i = 0; of = -1.; /* far away */ 371 while (i < max_it) { 372 i++; c2--; c3 = i * c2; c4 += 2; 373 /* c2 = y - i, c3 = i(y - i), c4 = d + 2i, for i odd */ 374 a1 = c4 * a2 + c3 * a1; 375 b1 = c4 * b2 + c3 * b1; 376 377 i++; c2--; c3 = i * c2; c4 += 2; 378 /* c2 = y - i, c3 = i(y - i), c4 = d + 2i, for i even */ 379 a2 = c4 * a1 + c3 * a2; 380 b2 = c4 * b1 + c3 * b2; 381 382 mixin (NEEDED_SCALE!("if")()); 383 384 if (b2 != 0) { 385 f = a2 / b2; 386 /* convergence check: relative; "absolute" for very small f : */ 387 if (fabs (f - of) <= EPSILON!T * fmax2(f0, fabs(f))) { 388 return f; 389 } 390 of = f; 391 } 392 } 393 394 return f;/* should not happen ... */ 395 } 396 397 398 static T pd_lower_series(T)(T lambda, T y) 399 { 400 T term = 1, sum = 0; 401 402 while (y >= 1 && term > sum * EPSILON!T) { 403 term *= y / lambda; 404 sum += term; 405 y--; 406 } 407 /* sum = \sum_{n=0}^ oo y*(y-1)*...*(y - n) / lambda^(n+1) 408 * = y/lambda * (1 + \sum_{n=1}^Inf (y-1)*...*(y-n) / lambda^n) 409 * ~ y/lambda + o(y/lambda) 410 */ 411 412 if (y != floor(y)) { 413 /* 414 * The series does not converge as the terms start getting 415 * bigger (besides flipping sign) for y < -lambda. 416 */ 417 T f; 418 /* FIXME: in quite few cases, adding term*f has no effect (f too small) 419 * and is unnecessary e.g. for pgamma(4e12, 121.1) */ 420 f = pd_lower_cf!T(y, lambda + 1 - y); 421 sum += term * f; 422 } 423 424 return sum; 425 } /* pd_lower_series() */ 426 427 428 /* 429 * Compute the following ratio with higher accuracy that would be had 430 * from doing it directly. 431 * 432 * dnorm (x, 0, 1, FALSE) 433 * ---------------------------------- 434 * pnorm (x, 0, 1, lower_tail, FALSE) 435 * 436 * Abramowitz & Stegun 26.2.12 437 */ 438 static T dpnorm(T)(T x, int lower_tail, T lp) 439 { 440 /* 441 * So as not to repeat a pnorm call, we expect 442 * 443 * lp == pnorm (x, 0, 1, lower_tail, TRUE) 444 * 445 * but use it only in the non-critical case where either x is small 446 * or p==exp(lp) is close to 1. 447 */ 448 449 if (x < 0) { 450 x = -x; 451 lower_tail = !lower_tail; 452 } 453 454 if (x > 10 && !lower_tail) { 455 T term = 1 / x; 456 T sum_ = term; 457 T x2 = x * x; 458 T i = 1; 459 460 do { 461 term *= -i / x2; 462 sum_ += term; 463 i += 2; 464 } while(fabs(term) > EPSILON!T*sum_); 465 466 return 1/sum_; 467 } else { 468 T d = dnorm(x, 0., 1., 0); 469 return d/exp(lp); 470 } 471 } 472 473 /* 474 * Asymptotic expansion to calculate the probability that Poisson variate 475 * has value <= x. 476 * Various assertions about this are made (without proof) at 477 * http://members.aol.com/iandjmsmith/PoissonApprox.htm 478 */ 479 static T ppois_asymp(T)(T x, T lambda, int lower_tail, int log_p) 480 { 481 static const(T)[8] coefs_a = [ 482 -1e99, /* placeholder used for 1-indexing */ 483 2/3., 484 -4/135., 485 8/2835., 486 16/8505., 487 -8992/12629925., 488 -334144/492567075., 489 698752/1477701225. 490 ]; 491 492 static const(T)[8] coefs_b = [ 493 -1e99, /* placeholder */ 494 1/12., 495 1/288., 496 -139/51840., 497 -571/2488320., 498 163879/209018880., 499 5246819/75246796800., 500 -534703531/902961561600. 501 ]; 502 503 T elfb, elfb_term; 504 T res12, res1_term, res1_ig, res2_term, res2_ig; 505 T dfm, pt_, s2pt, f, np; 506 int i; 507 508 dfm = lambda - x; 509 /* If lambda is large, the distribution is highly concentrated 510 about lambda. So representation error in x or lambda can lead 511 to arbitrarily large values of pt_ and hence divergence of the 512 coefficients of this approximation. 513 */ 514 pt_ = - log1pmx (dfm / x); 515 s2pt = sqrt (2 * x * pt_); 516 if (dfm < 0) 517 s2pt = -s2pt; 518 519 res12 = 0; 520 res1_ig = res1_term = sqrt (x); 521 res2_ig = res2_term = s2pt; 522 for (i = 1; i < 8; i++) { 523 res12 += res1_ig * coefs_a[i]; 524 res12 += res2_ig * coefs_b[i]; 525 res1_term *= pt_ / i ; 526 res2_term *= 2 * pt_ / (2 * i + 1); 527 res1_ig = res1_ig / x + res1_term; 528 res2_ig = res2_ig / x + res2_term; 529 } 530 531 elfb = x; 532 elfb_term = 1; 533 for (i = 1; i < 8; i++) { 534 elfb += elfb_term * coefs_b[i]; 535 elfb_term /= x; 536 } 537 if (!lower_tail) 538 elfb = -elfb; 539 540 f = res12 / elfb; 541 542 np = pnorm(s2pt, 0.0, 1.0, !lower_tail, log_p); 543 544 if (log_p) { 545 T n_d_over_p = dpnorm (s2pt, !lower_tail, np); 546 return np + log1p (f * n_d_over_p); 547 } else { 548 T nd = dnorm (s2pt, 0., 1., log_p); 549 return np + f * nd; 550 } 551 } /* ppois_asymp() */ 552 553 T pgamma_raw(T)(T x, T alph, int lower_tail, int log_p) 554 { 555 /* Here, assume that (x,alph) are not NA & alph > 0 . */ 556 557 T res; 558 mixin R_D__1!log_p; 559 560 mixin R_P_bounds_01!(x, 0., T.infinity); 561 562 if (x < 1) { 563 res = pgamma_smallx!T(x, alph, lower_tail, log_p); 564 } else if (x <= alph - 1 && x < 0.8 * (alph + 50)) { 565 /* incl. large alph compared to x */ 566 T sum = pd_upper_series!T(x, alph, log_p);/* = x/alph + o(x/alph) */ 567 T d = dpois_wrap!T(alph, x, log_p); 568 if (!lower_tail) 569 res = log_p ? R_Log1_Exp!T(d + sum) : 1 - d * sum; 570 else 571 res = log_p ? sum + d : sum * d; 572 } else if (alph - 1 < x && alph < 0.8 * (x + 50)) { 573 /* incl. large x compared to alph */ 574 T sum; 575 T d = dpois_wrap!T(alph, x, log_p); 576 if (alph < 1) { 577 if (x * EPSILON!T > 1 - alph) 578 sum = R_D__1; 579 else { 580 T f = pd_lower_cf!T(alph, x - (alph - 1)) * x / alph; 581 /* = [alph/(x - alph+1) + o(alph/(x-alph+1))] * x/alph = 1 + o(1) */ 582 sum = log_p ? log (f) : f; 583 } 584 } else { 585 sum = pd_lower_series!T(x, alph - 1);/* = (alph-1)/x + o((alph-1)/x) */ 586 sum = log_p ? log1p (sum) : 1 + sum; 587 } 588 if (!lower_tail) 589 res = log_p ? sum + d : sum * d; 590 else 591 res = log_p ? R_Log1_Exp!T(d + sum) : 1 - d * sum; 592 } else { /* x >= 1 and x fairly near alph. */ 593 res = ppois_asymp!T(alph - 1, x, !lower_tail, log_p); 594 } 595 596 /* 597 * We lose a fair amount of accuracy to underflow in the cases 598 * where the final result is very close to DBL_MIN. In those 599 * cases, simply redo via log space. 600 */ 601 if (!log_p && res < MIN!T / EPSILON!T) { 602 /* with(.Machine, double.xmin / double.eps) #|-> 1.002084e-292 */ 603 return exp(pgamma_raw!T(x, alph, lower_tail, 1)); 604 } else 605 return res; 606 } 607 608 609 T pgamma(T: double)(T x, T alph, T scale, int lower_tail, int log_p) 610 { 611 if (isNaN(x) || isNaN(alph) || isNaN(scale)) 612 return x + alph + scale; 613 if(alph < 0. || scale <= 0.) 614 return T.nan; 615 x /= scale; 616 617 if (isNaN(x)) /* eg. original x = scale = +Inf */ 618 return x; 619 620 if(alph == 0.) /* limit case; useful e.g. in pnchisq() */ 621 return (x <= 0) ? R_DT_0!T(lower_tail, log_p): R_DT_1!T(lower_tail, log_p); /* <= assert pgamma(0,0) ==> 0 */ 622 return pgamma_raw!T(x, alph, lower_tail, log_p); 623 } 624 625 626 T qchisq_appr(T)(T p, T nu, T g /* = log Gamma(nu/2) */, 627 int lower_tail, int log_p, T tol /* EPS1 */) 628 { 629 enum C7 = 4.67; 630 enum C8 = 6.66; 631 enum C9 = 6.73; 632 enum C10 = 13.32; 633 634 T alpha, a, c, ch, p1; 635 T p2, q, t, x; 636 637 /* test arguments and initialise */ 638 639 if (isNaN(p) || isNaN(nu)) 640 return p + nu; 641 mixin (R_Q_P01_check!(p)()); 642 if (nu <= 0) 643 return T.nan; 644 645 alpha = 0.5 * nu;/* = [pq]gamma() shape */ 646 c = alpha - 1; 647 648 if(nu < (-1.24)*(p1 = R_DT_log!T(p, lower_tail))) { /* for small chi-squared */ 649 /* log(alpha) + g = log(alpha) + log(gamma(alpha)) = 650 * = log(alpha*gamma(alpha)) = lgamma(alpha+1) suffers from 651 * catastrophic cancellation when alpha << 1 652 */ 653 T lgam1pa = (alpha < 0.5) ? lgamma1p!T(alpha) : (log(alpha) + g); 654 ch = exp((lgam1pa + p1)/alpha + M_LN2); 655 } else if(nu > 0.32) { /* using Wilson and Hilferty estimate */ 656 657 x = qnorm!T(p, 0, 1, lower_tail, log_p); 658 p1 = 2./(9*nu); 659 ch = nu*pow(x*sqrt(p1) + 1-p1, 3); 660 661 /* approximation for p tending to 1: */ 662 if( ch > 2.2*nu + 6 ) 663 ch = -2*(R_DT_Clog!T(p, lower_tail, log_p) - c*log(0.5*ch) + g); 664 665 } else { /* "small nu" : 1.24*(-log(p)) <= nu <= 0.32 */ 666 667 ch = 0.4; 668 a = R_DT_Clog!T(p, lower_tail, log_p) + g + c*M_LN2; 669 670 do { 671 q = ch; 672 p1 = 1. / (1 + ch*(C7 + ch)); 673 p2 = ch*(C9 + ch*(C8 + ch)); 674 t = -0.5 +(C7 + 2*ch)*p1 - (C9 + ch*(C10 + 3*ch))/p2; 675 ch -= (1 - exp(a + 0.5*ch)*p2*p1)/t; 676 } while(fabs(q - ch) > tol * fabs(ch)); 677 } 678 679 return ch; 680 } 681 682 T qgamma(T: double)(T p, T alpha, T scale, int lower_tail, int log_p) 683 /* shape = alpha */ 684 { 685 enum EPS1 = 1e-2; 686 enum EPS2 = 5e-7; /* final precision of AS 91 */ 687 enum EPS_N = 1e-15; /* precision of Newton step / iterations */ 688 enum LN_EPS = -36.043653389117156; /* = log(.Machine$double.eps) iff IEEE_754 */ 689 690 enum MAXIT = 1000; /* was 20 */ 691 692 enum pMIN = 1e-100; /* was 0.000002 = 2e-6 */ 693 enum pMAX = (1-1e-14); /* was (1-1e-12) and 0.999998 = 1 - 2e-6 */ 694 695 mixin R_D__0!log_p; 696 697 const static T 698 i420 = 1./ 420., 699 i2520 = 1./ 2520., 700 i5040 = 1./ 5040; 701 702 T p_, a, b, c, g, ch, ch0, p1; 703 T p2, q, s1, s2, s3, s4, s5, s6, t, x; 704 int i, max_it_Newton = 1; 705 706 /* test arguments and initialise */ 707 708 if (isNaN(p) || isNaN(alpha) || isNaN(scale)) 709 return p + alpha + scale; 710 711 mixin R_Q_P01_boundaries!(p, 0., T.infinity); 712 713 if (alpha < 0 || scale <= 0) 714 return T.nan; 715 716 if (alpha == 0) /* all mass at 0 : */ 717 return 0.; 718 719 if (alpha < 1e-10) { 720 /* Warning seems unnecessary now: */ 721 max_it_Newton = 7;/* may still be increased below */ 722 } 723 724 mixin R_DT_qIv!(p); 725 p_ = R_DT_qIv;/* lower_tail prob (in any case) */ 726 727 g = lgammafn!T(alpha);/* log Gamma(v/2) */ 728 729 /*----- Phase I : Starting Approximation */ 730 ch = qchisq_appr!T(p, /* nu= 'df' = */ 2*alpha, /* lgamma(nu/2)= */ g, 731 lower_tail, log_p, /* tol= */ EPS1); 732 if(!isFinite(ch)) { 733 /* forget about all iterations! */ 734 max_it_Newton = 0; goto END; 735 } 736 if(ch < EPS2) {/* Corrected according to AS 91; MM, May 25, 1999 */ 737 max_it_Newton = 20; 738 goto END;/* and do Newton steps */ 739 } 740 741 /* FIXME: This (cutoff to {0, +Inf}) is far from optimal 742 * ----- when log_p or !lower_tail, but NOT doing it can be even worse */ 743 if(p_ > pMAX || p_ < pMIN) { 744 /* did return ML_POSINF or 0.; much better: */ 745 max_it_Newton = 20; 746 goto END;/* and do Newton steps */ 747 } 748 749 /*----- Phase II: Iteration 750 * Call pgamma() [AS 239] and calculate seven term taylor series 751 */ 752 c = alpha - 1; 753 s6 = (120 + c*(346 + 127*c)) * i5040; /* used below, is "const" */ 754 755 ch0 = ch;/* save initial approx. */ 756 for(i=1; i <= MAXIT; i++ ) { 757 q = ch; 758 p1 = 0.5*ch; 759 p2 = p_ - pgamma_raw!T(p1, alpha, /*lower_tail*/1, /*log_p*/0); 760 761 if(!isFinite(p2) || ch <= 0) 762 { 763 ch = ch0; max_it_Newton = 27; goto END; 764 }/*was return ML_NAN;*/ 765 766 t = p2*exp(alpha*M_LN2 + g + p1 - c*log(ch)); 767 b = t/ch; 768 a = 0.5*t - b*c; 769 s1 = (210 + a*(140 + a*(105 + a*(84 + a*(70 + 60*a))))) * i420; 770 s2 = (420 + a*(735 + a*(966 + a*(1141 + 1278*a)))) * i2520; 771 s3 = (210 + a*(462 + a*(707 + 932*a))) * i2520; 772 s4 = (252 + a*(672 + 1182*a) + c*(294 + a*(889 + 1740*a))) * i5040; 773 s5 = (84 + 2264*a + c*(1175 + 606*a)) * i2520; 774 775 ch += t*(1 + 0.5*t*s1 - b*c*(s1 - b*(s2 - b*(s3 - b*(s4 - b*(s5 - b*s6)))))); 776 if(fabs(q - ch) < EPS2*ch) 777 goto END; 778 if(fabs(q - ch) > 0.1*ch) {/* diverging? -- also forces ch > 0 */ 779 if(ch < q) ch = 0.9 * q; else ch = 1.1 * q; 780 } 781 } 782 /* no convergence in MAXIT iterations -- but we add Newton now... */ 783 /* was 784 * ML_ERROR(ME_PRECISION, "qgamma"); 785 * does nothing in R !*/ 786 787 END: 788 /* PR# 2214 : From: Morten Welinder <terra@diku.dk>, Fri, 25 Oct 2002 16:50 789 -------- To: R-bugs@biostat.ku.dk Subject: qgamma precision 790 791 * With a final Newton step, double accuracy, e.g. for (p= 7e-4; nu= 0.9) 792 * 793 * Improved (MM): - only if rel.Err > EPS_N (= 1e-15); 794 * - also for lower_tail = FALSE or log_p = TRUE 795 * - optionally *iterate* Newton 796 */ 797 x = 0.5*scale*ch; 798 if(max_it_Newton) { 799 /* always use log scale */ 800 if (!log_p) { 801 p = log(p); 802 log_p = 1; 803 } 804 if(x == 0) { 805 const T _1_p = 1. + 1e-7; 806 const T _1_m = 1. - 1e-7; 807 x = MIN!T; 808 p_ = pgamma!T(x, alpha, scale, lower_tail, log_p); 809 if(( lower_tail && p_ > p * _1_p) || 810 (!lower_tail && p_ < p * _1_m)) 811 return(0.); 812 /* else: continue, using x = DBL_MIN instead of 0 */ 813 } 814 else 815 p_ = pgamma!T(x, alpha, scale, lower_tail, log_p); 816 if(p_ == -T.infinity) 817 return 0; /* PR#14710 */ 818 for(i = 1; i <= max_it_Newton; i++) { 819 p1 = p_ - p; 820 if(fabs(p1) < fabs(EPS_N * p)) 821 break; 822 /* else */ 823 if((g = dgamma!T(x, alpha, scale, log_p)) == R_D__0) { 824 break; 825 } 826 /* else : 827 * delta x = f(x)/f'(x); 828 * if(log_p) f(x) := log P(x) - p; f'(x) = d/dx log P(x) = P' / P 829 * ==> f(x)/f'(x) = f*P / P' = f*exp(p_) / P' (since p_ = log P(x)) 830 */ 831 t = log_p ? p1*exp(p_ - g) : p1/g ;/* = "delta x" */ 832 t = lower_tail ? x - t : x + t; 833 p_ = pgamma!T(t, alpha, scale, lower_tail, log_p); 834 if (fabs(p_ - p) > fabs(p1) || 835 (i > 1 && fabs(p_ - p) == fabs(p1)) /* <- against flip-flop */) { 836 /* no improvement */ 837 break; 838 } /* else : */ 839 //#ifdef Harmful_notably_if_max_it_Newton_is_1 840 /* control step length: this could have started at 841 the initial approximation */ 842 //if(t > 1.1*x) t = 1.1*x; 843 //else if(t < 0.9*x) t = 0.9*x; 844 //#endif 845 x = t; 846 } 847 } 848 849 return x; 850 } 851 852 853 T rgamma(T: double)(T a, T scale) 854 { 855 /* Constants : */ 856 const static T sqrt32 = 5.656854; 857 const static T exp_m1 = 0.36787944117144232159;/* exp(-1) = 1/e */ 858 859 /* Coefficients q[k] - for q0 = sum(q[k]*a^(-k)) 860 * Coefficients a[k] - for q = q0+(t*t/2)*sum(a[k]*v^k) 861 * Coefficients e[k] - for exp(q)-1 = sum(e[k]*q^k) 862 */ 863 const static T q1 = 0.04166669; 864 const static T q2 = 0.02083148; 865 const static T q3 = 0.00801191; 866 const static T q4 = 0.00144121; 867 const static T q5 = -7.388e-5; 868 const static T q6 = 2.4511e-4; 869 const static T q7 = 2.424e-4; 870 871 const static T a1 = 0.3333333; 872 const static T a2 = -0.250003; 873 const static T a3 = 0.2000062; 874 const static T a4 = -0.1662921; 875 const static T a5 = 0.1423657; 876 const static T a6 = -0.1367177; 877 const static T a7 = 0.1233795; 878 879 /* State variables [FIXME for threading!] :*/ 880 static T aa = 0.; 881 static T aaa = 0.; 882 static T s, s2, d; /* no. 1 (step 1) */ 883 static T q0, b, si, c;/* no. 2 (step 4) */ 884 885 T e, p, q, r, t, u, v, w, x, ret_val; 886 887 if (isNaN(a) || isNaN(scale)) 888 return T.nan; 889 if (a <= 0.0 || scale <= 0.0) { 890 if(scale == 0. || a == 0.) return 0.; 891 return T.nan; 892 } 893 if(!isFinite(a) || !isFinite(scale)) 894 return T.infinity; 895 896 if (a < 1.) { /* GS algorithm for parameters a < 1 */ 897 e = 1.0 + exp_m1 * a; 898 for(;;) { 899 p = e * unif_rand!T(); 900 if (p >= 1.0) { 901 x = -log((e - p) / a); 902 if (exp_rand!T() >= (1.0 - a) * log(x)) 903 break; 904 } else { 905 x = exp(log(p) / a); 906 if (exp_rand!T() >= x) 907 break; 908 } 909 } 910 return scale * x; 911 } 912 913 /* --- a >= 1 : GD algorithm --- */ 914 915 /* Step 1: Recalculations of s2, s, d if a has changed */ 916 if (a != aa) { 917 aa = a; 918 s2 = a - 0.5; 919 s = sqrt(s2); 920 d = sqrt32 - s * 12.0; 921 } 922 /* Step 2: t = standard normal deviate, 923 x = (s,1/2) -normal deviate. */ 924 925 /* immediate acceptance (i) */ 926 t = norm_rand!T(); 927 x = s + 0.5 * t; 928 ret_val = x * x; 929 if (t >= 0.0) 930 return scale * ret_val; 931 932 /* Step 3: u = 0,1 - uniform sample. squeeze acceptance (s) */ 933 u = unif_rand!T(); 934 if (d * u <= t * t * t) 935 return scale * ret_val; 936 937 /* Step 4: recalculations of q0, b, si, c if necessary */ 938 939 if (a != aaa) { 940 aaa = a; 941 r = 1.0 / a; 942 q0 = ((((((q7 * r + q6) * r + q5) * r + q4) * r + q3) * r 943 + q2) * r + q1) * r; 944 945 /* Approximation depending on size of parameter a */ 946 /* The constants in the expressions for b, si and c */ 947 /* were established by numerical experiments */ 948 949 if (a <= 3.686) { 950 b = 0.463 + s + 0.178 * s2; 951 si = 1.235; 952 c = 0.195 / s - 0.079 + 0.16 * s; 953 } else if (a <= 13.022) { 954 b = 1.654 + 0.0076 * s2; 955 si = 1.68 / s + 0.275; 956 c = 0.062 / s + 0.024; 957 } else { 958 b = 1.77; 959 si = 0.75; 960 c = 0.1515 / s; 961 } 962 } 963 /* Step 5: no quotient test if x not positive */ 964 965 if (x > 0.0) { 966 /* Step 6: calculation of v and quotient q */ 967 v = t / (s + s); 968 if (fabs(v) <= 0.25) 969 q = q0 + 0.5 * t * t * ((((((a7 * v + a6) * v + a5) * v + a4) * v 970 + a3) * v + a2) * v + a1) * v; 971 else 972 q = q0 - s * t + 0.25 * t * t + (s2 + s2) * log(1.0 + v); 973 974 975 /* Step 7: quotient acceptance (q) */ 976 if (log(1.0 - u) <= q) 977 return scale * ret_val; 978 } 979 980 for(;;) { 981 /* Step 8: e = standard exponential deviate 982 * u = 0,1 -uniform deviate 983 * t = (b,si)-double exponential (laplace) sample */ 984 e = exp_rand!T(); 985 u = unif_rand!T(); 986 u = u + u - 1.0; 987 if (u < 0.0) 988 t = b - si * e; 989 else 990 t = b + si * e; 991 /* Step 9: rejection if t < tau(1) = -0.71874483771719 */ 992 if (t >= -0.71874483771719) { 993 /* Step 10: calculation of v and quotient q */ 994 v = t / (s + s); 995 if (fabs(v) <= 0.25) 996 q = q0 + 0.5 * t * t * 997 ((((((a7 * v + a6) * v + a5) * v + a4) * v + a3) * v 998 + a2) * v + a1) * v; 999 else 1000 q = q0 - s * t + 0.25 * t * t + (s2 + s2) * log(1.0 + v); 1001 /* Step 11: hat acceptance (h) */ 1002 /* (if q not positive go to step 8) */ 1003 if (q > 0.0) { 1004 w = expm1(q); 1005 /* ^^^^^ original code had approximation with rel.err < 2e-7 */ 1006 /* if t is rejected sample again at step 8 */ 1007 if (c * fabs(u) <= w * exp(e - 0.5 * t * t)) 1008 break; 1009 } 1010 } 1011 } /* repeat .. until `t' is accepted */ 1012 x = s + 0.5 * t; 1013 return scale * x * x; 1014 } 1015 1016 1017 void test_gamma() 1018 { 1019 import std.stdio: write, writeln; 1020 writeln("dgamma: ", dgamma(1., 1., 1., 0)); 1021 writeln("pgamma: ", pgamma(1., 1., 1., 1, 0)); 1022 writeln("qgamma: ", qgamma(0.5, 1., 1., 1, 0)); 1023 writeln("rgamma: ", rgamma(1., 1.), ", rgamma: ", rgamma(1., 1.) , ", rgamma: ", rgamma(1., 1.)); 1024 }