|
Scythe-1.0.3
|
00001 /* 00002 * Scythe Statistical Library Copyright (C) 2000-2002 Andrew D. Martin 00003 * and Kevin M. Quinn; 2002-present Andrew D. Martin, Kevin M. Quinn, 00004 * and Daniel Pemstein. All Rights Reserved. 00005 * 00006 * This program is free software; you can redistribute it and/or 00007 * modify under the terms of the GNU General Public License as 00008 * published by Free Software Foundation; either version 2 of the 00009 * License, or (at your option) any later version. See the text files 00010 * COPYING and LICENSE, distributed with this source code, for further 00011 * information. 00012 * -------------------------------------------------------------------- 00013 * scythestat/distributions.h 00014 * 00015 */ 00016 00069 /* TODO: Figure out how to get doxygen to stop listing erroneous 00070 * variables at the end of the doc for this file. They stem from it 00071 * misreading the nested macro calls used to generate matrix procs. 00072 */ 00073 00074 /* TODO: Full R-style versions of these function including arbitrary 00075 * recycling of matrix arguments. This is going to have to wait for 00076 * variadic templates to be doable without a complete mess. There is 00077 * currently a variadic patch available for g++, perhaps we can add a 00078 * SCYTHE_VARIADIC flag and include these as option until they become 00079 * part of the standard in 2009. Something to come back to after 1.0. 00080 */ 00081 00082 #ifndef SCYTHE_DISTRIBUTIONS_H 00083 #define SCYTHE_DISTRIBUTIONS_H 00084 00085 #include <iostream> 00086 #include <cmath> 00087 #include <cfloat> 00088 #include <climits> 00089 #include <algorithm> 00090 #include <limits> 00091 00092 #ifdef HAVE_IEEEFP_H 00093 #include <ieeefp.h> 00094 #endif 00095 00096 #ifdef SCYTHE_COMPILE_DIRECT 00097 #include "matrix.h" 00098 #include "ide.h" 00099 #include "error.h" 00100 #else 00101 #include "scythestat/matrix.h" 00102 #include "scythestat/ide.h" 00103 #include "scythestat/error.h" 00104 #endif 00105 00106 /* Fill in some defs from R that aren't in math.h */ 00107 #ifndef M_PI 00108 #define M_PI 3.141592653589793238462643383280 00109 #endif 00110 #define M_LN_SQRT_2PI 0.918938533204672741780329736406 00111 #define M_LN_SQRT_PId2 0.225791352644727432363097614947 00112 #define M_1_SQRT_2PI 0.39894228040143267793994605993 00113 #define M_2PI 6.28318530717958647692528676655 00114 #define M_SQRT_32 5.656854249492380195206754896838 00115 00116 #ifndef HAVE_TRUNC 00117 00118 inline double trunc(double x) throw () 00119 { 00120 if (x >= 0) 00121 return std::floor(x); 00122 else 00123 return std::ceil(x); 00124 } 00126 #endif 00127 00128 /* Many random number generators, pdfs, cdfs, and functions (gamma, 00129 * etc) in this file are based on code from the R Project, version 00130 * 1.6.0-1.7.1. This code is available under the terms of the GNU 00131 * GPL. Original copyright: 00132 * 00133 * Copyright (C) 1998 Ross Ihaka 00134 * Copyright (C) 2000-2002 The R Development Core Team 00135 * Copyright (C) 2003 The R Foundation 00136 */ 00137 00138 namespace scythe { 00139 00142 /* Forward declarations */ 00143 double gammafn (double); 00144 double lngammafn (double); 00145 double lnbetafn (double, double); 00146 double pgamma (double, double, double); 00147 double dgamma(double, double, double); 00148 double pnorm (double, double, double); 00149 00152 /******************** 00153 * Helper Functions * 00154 ********************/ 00155 namespace { 00156 00157 /* Evaluate a Chebysheve series at a given point */ 00158 double 00159 chebyshev_eval (double x, const double *a, int n) 00160 { 00161 SCYTHE_CHECK_10(n < 1 || n > 1000, scythe_invalid_arg, 00162 "n not on [1, 1000]"); 00163 00164 SCYTHE_CHECK_10(x < -1.1 || x > 1.1, scythe_invalid_arg, 00165 "x not on [-1.1, 1.1]"); 00166 00167 double b0, b1, b2; 00168 b0 = b1 = b2 = 0; 00169 00170 double twox = x * 2; 00171 00172 for (int i = 1; i <= n; ++i) { 00173 b2 = b1; 00174 b1 = b0; 00175 b0 = twox * b1 - b2 + a[n - i]; 00176 } 00177 00178 return (b0 - b2) * 0.5; 00179 } 00180 00181 /* Computes the log gamma correction factor for x >= 10 */ 00182 double 00183 lngammacor(double x) 00184 { 00185 const double algmcs[15] = { 00186 +.1666389480451863247205729650822e+0, 00187 -.1384948176067563840732986059135e-4, 00188 +.9810825646924729426157171547487e-8, 00189 -.1809129475572494194263306266719e-10, 00190 +.6221098041892605227126015543416e-13, 00191 }; 00192 00193 SCYTHE_CHECK_10(x < 10, scythe_invalid_arg, 00194 "This function requires x >= 10"); 00195 SCYTHE_CHECK_10(x >= 3.745194030963158e306, scythe_range_error, 00196 "Underflow"); 00197 00198 if (x < 94906265.62425156) { 00199 double tmp = 10 / x; 00200 return chebyshev_eval(tmp * tmp * 2 - 1, algmcs, 5) / x; 00201 } 00202 00203 return 1 / (x * 12); 00204 } 00205 00206 /* Evaluates the "deviance part" */ 00207 double 00208 bd0(double x, double np) 00209 { 00210 00211 if(std::fabs(x - np) < 0.1 * (x + np)) { 00212 double v = (x - np) / (x + np); 00213 double s = (x - np) * v; 00214 double ej = 2 * x * v; 00215 v = v * v; 00216 for (int j = 1; ; j++) { 00217 ej *= v; 00218 double s1 = s + ej / ((j << 1) + 1); 00219 if (s1 == s) 00220 return s1; 00221 s = s1; 00222 } 00223 } 00224 00225 return x * std::log(x / np) + np - x; 00226 } 00227 00228 /* Computes the log of the error term in Stirling's formula */ 00229 double 00230 stirlerr(double n) 00231 { 00232 #define S0 0.083333333333333333333 /* 1/12 */ 00233 #define S1 0.00277777777777777777778 /* 1/360 */ 00234 #define S2 0.00079365079365079365079365 /* 1/1260 */ 00235 #define S3 0.000595238095238095238095238 /* 1/1680 */ 00236 #define S4 0.0008417508417508417508417508/* 1/1188 */ 00237 00238 /* error for 0, 0.5, 1.0, 1.5, ..., 14.5, 15.0 */ 00239 const double sferr_halves[31] = { 00240 0.0, /* n=0 - wrong, place holder only */ 00241 0.1534264097200273452913848, /* 0.5 */ 00242 0.0810614667953272582196702, /* 1.0 */ 00243 0.0548141210519176538961390, /* 1.5 */ 00244 0.0413406959554092940938221, /* 2.0 */ 00245 0.03316287351993628748511048, /* 2.5 */ 00246 0.02767792568499833914878929, /* 3.0 */ 00247 0.02374616365629749597132920, /* 3.5 */ 00248 0.02079067210376509311152277, /* 4.0 */ 00249 0.01848845053267318523077934, /* 4.5 */ 00250 0.01664469118982119216319487, /* 5.0 */ 00251 0.01513497322191737887351255, /* 5.5 */ 00252 0.01387612882307074799874573, /* 6.0 */ 00253 0.01281046524292022692424986, /* 6.5 */ 00254 0.01189670994589177009505572, /* 7.0 */ 00255 0.01110455975820691732662991, /* 7.5 */ 00256 0.010411265261972096497478567, /* 8.0 */ 00257 0.009799416126158803298389475, /* 8.5 */ 00258 0.009255462182712732917728637, /* 9.0 */ 00259 0.008768700134139385462952823, /* 9.5 */ 00260 0.008330563433362871256469318, /* 10.0 */ 00261 0.007934114564314020547248100, /* 10.5 */ 00262 0.007573675487951840794972024, /* 11.0 */ 00263 0.007244554301320383179543912, /* 11.5 */ 00264 0.006942840107209529865664152, /* 12.0 */ 00265 0.006665247032707682442354394, /* 12.5 */ 00266 0.006408994188004207068439631, /* 13.0 */ 00267 0.006171712263039457647532867, /* 13.5 */ 00268 0.005951370112758847735624416, /* 14.0 */ 00269 0.005746216513010115682023589, /* 14.5 */ 00270 0.005554733551962801371038690 /* 15.0 */ 00271 }; 00272 double nn; 00273 00274 if (n <= 15.0) { 00275 nn = n + n; 00276 if (nn == (int)nn) 00277 return(sferr_halves[(int)nn]); 00278 return (scythe::lngammafn(n + 1.) - (n + 0.5) * std::log(n) + n - 00279 std::log(std::sqrt(2 * M_PI))); 00280 } 00281 00282 nn = n*n; 00283 if (n > 500) 00284 return((S0 - S1 / nn) / n); 00285 if (n > 80) 00286 return((S0 - (S1 - S2 / nn) / nn) / n); 00287 if (n > 35) 00288 return((S0 - (S1 - (S2 - S3 / nn) / nn) / nn) / n); 00289 /* 15 < n <= 35 : */ 00290 return((S0 - (S1 - (S2 - (S3 - S4 / nn) / nn) / nn) / nn) / n); 00291 #undef S1 00292 #undef S2 00293 #undef S3 00294 #undef S4 00295 } 00296 00297 00298 /* Helper for dpois and dgamma */ 00299 double 00300 dpois_raw (double x, double lambda) 00301 { 00302 if (lambda == 0) 00303 return ( (x == 0) ? 1.0 : 0.0); 00304 00305 if (x == 0) 00306 return std::exp(-lambda); 00307 00308 if (x < 0) 00309 return 0.0; 00310 00311 return std::exp(-stirlerr(x) - bd0(x, lambda)) 00312 / std::sqrt(2 * M_PI * x); 00313 } 00314 00315 00316 /* helper for pbeta */ 00317 double 00318 pbeta_raw(double x, double pin, double qin) 00319 { 00320 double ans, c, finsum, p, ps, p1, q, term, xb, xi, y; 00321 int n, i, ib, swap_tail; 00322 00323 const double eps = .5 * DBL_EPSILON; 00324 const double sml = DBL_MIN; 00325 const double lneps = std::log(eps); 00326 const double lnsml = std::log(eps); 00327 00328 if (pin / (pin + qin) < x) { 00329 swap_tail = 1; 00330 y = 1 - x; 00331 p = qin; 00332 q = pin; 00333 } else { 00334 swap_tail=0; 00335 y = x; 00336 p = pin; 00337 q = qin; 00338 } 00339 00340 if ((p + q) * y / (p + 1) < eps) { 00341 ans = 0; 00342 xb = p * std::log(std::max(y,sml)) - std::log(p) - lnbetafn(p,q); 00343 if (xb > lnsml && y != 0) 00344 ans = std::exp(xb); 00345 if (swap_tail) 00346 ans = 1-ans; 00347 } else { 00348 ps = q - std::floor(q); 00349 if (ps == 0) 00350 ps = 1; 00351 xb = p * std::log(y) - lnbetafn(ps, p) - std::log(p); 00352 ans = 0; 00353 if (xb >= lnsml) { 00354 ans = std::exp(xb); 00355 term = ans * p; 00356 if (ps != 1) { 00357 n = (int)std::max(lneps/std::log(y), 4.0); 00358 for(i = 1; i <= n; i++){ 00359 xi = i; 00360 term *= (xi-ps)*y/xi; 00361 ans += term/(p+xi); 00362 } 00363 } 00364 } 00365 if (q > 1) { 00366 xb = p * std::log(y) + q * std::log(1 - y) 00367 - lnbetafn(p, q) - std::log(q); 00368 ib = (int) std::max(xb / lnsml, 0.0); 00369 term = std::exp(xb - ib * lnsml); 00370 c = 1 / (1 - y); 00371 p1 = q * c / (p + q - 1); 00372 00373 finsum = 0; 00374 n = (int) q; 00375 if(q == n) 00376 n--; 00377 for (i = 1; i <= n; i++) { 00378 if(p1 <= 1 && term / eps <= finsum) 00379 break; 00380 xi = i; 00381 term = (q -xi + 1) * c * term / (p + q - xi); 00382 if (term > 1) { 00383 ib--; 00384 term *= sml; 00385 } 00386 if (ib == 0) 00387 finsum += term; 00388 } 00389 ans += finsum; 00390 } 00391 00392 if(swap_tail) 00393 ans = 1-ans; 00394 ans = std::max(std::min(ans,1.),0.); 00395 } 00396 return ans; 00397 } 00398 00399 /* Helper for dbinom */ 00400 double 00401 dbinom_raw (double x, double n, double p, double q) 00402 { 00403 double f, lc; 00404 00405 if (p == 0) 00406 return((x == 0) ? 1.0 : 0.0); 00407 if (q == 0) 00408 return((x == n) ? 1.0 : 0.0); 00409 00410 if (x == 0) { 00411 if(n == 0) 00412 return 1.0; 00413 00414 lc = (p < 0.1) ? -bd0(n, n * q) - n * p : n * std::log(q); 00415 return(std::exp(lc)); 00416 } 00417 if (x == n) { 00418 lc = (q < 0.1) ? -bd0(n,n * p) - n * q : n * std::log(p); 00419 return(std::exp(lc)); 00420 } 00421 00422 if (x < 0 || x > n) 00423 return 0.0; 00424 00425 lc = stirlerr(n) - stirlerr(x) - stirlerr(n-x) - bd0(x,n*p) - 00426 bd0(n - x, n * q); 00427 00428 f = (M_2PI * x * (n-x)) / n; 00429 00430 return (std::exp(lc) / std::sqrt(f)); 00431 } 00432 00433 /* The normal probability density function implementation. */ 00434 00435 #define SIXTEN 16 00436 #define do_del(X) \ 00437 xsq = trunc(X * SIXTEN) / SIXTEN; \ 00438 del = (X - xsq) * (X + xsq); \ 00439 if(log_p) { \ 00440 *cum = (-xsq * xsq * 0.5) + (-del * 0.5) + std::log(temp); \ 00441 if((lower && x > 0.) || (upper && x <= 0.)) \ 00442 *ccum = ::log1p(-std::exp(-xsq * xsq * 0.5) * \ 00443 std::exp(-del * 0.5) * temp); \ 00444 } \ 00445 else { \ 00446 *cum = std::exp(-xsq * xsq * 0.5) * std::exp(-del * 0.5) * temp; \ 00447 *ccum = 1.0 - *cum; \ 00448 } 00449 00450 #define swap_tail \ 00451 if (x > 0.) {/* swap ccum <--> cum */ \ 00452 temp = *cum; if(lower) *cum = *ccum; *ccum = temp; \ 00453 } 00454 00455 void 00456 pnorm_both(double x, double *cum, double *ccum, int i_tail, 00457 bool log_p) 00458 { 00459 const double a[5] = { 00460 2.2352520354606839287, 00461 161.02823106855587881, 00462 1067.6894854603709582, 00463 18154.981253343561249, 00464 0.065682337918207449113 00465 }; 00466 const double b[4] = { 00467 47.20258190468824187, 00468 976.09855173777669322, 00469 10260.932208618978205, 00470 45507.789335026729956 00471 }; 00472 const double c[9] = { 00473 0.39894151208813466764, 00474 8.8831497943883759412, 00475 93.506656132177855979, 00476 597.27027639480026226, 00477 2494.5375852903726711, 00478 6848.1904505362823326, 00479 11602.651437647350124, 00480 9842.7148383839780218, 00481 1.0765576773720192317e-8 00482 }; 00483 const double d[8] = { 00484 22.266688044328115691, 00485 235.38790178262499861, 00486 1519.377599407554805, 00487 6485.558298266760755, 00488 18615.571640885098091, 00489 34900.952721145977266, 00490 38912.003286093271411, 00491 19685.429676859990727 00492 }; 00493 const double p[6] = { 00494 0.21589853405795699, 00495 0.1274011611602473639, 00496 0.022235277870649807, 00497 0.001421619193227893466, 00498 2.9112874951168792e-5, 00499 0.02307344176494017303 00500 }; 00501 const double q[5] = { 00502 1.28426009614491121, 00503 0.468238212480865118, 00504 0.0659881378689285515, 00505 0.00378239633202758244, 00506 7.29751555083966205e-5 00507 }; 00508 00509 double xden, xnum, temp, del, eps, xsq, y; 00510 int i, lower, upper; 00511 00512 /* Consider changing these : */ 00513 eps = DBL_EPSILON * 0.5; 00514 00515 /* i_tail in {0,1,2} =^= {lower, upper, both} */ 00516 lower = i_tail != 1; 00517 upper = i_tail != 0; 00518 00519 y = std::fabs(x); 00520 if (y <= 0.67448975) { 00521 /* qnorm(3/4) = .6744.... -- earlier had 0.66291 */ 00522 if (y > eps) { 00523 xsq = x * x; 00524 xnum = a[4] * xsq; 00525 xden = xsq; 00526 for (i = 0; i < 3; ++i) { 00527 xnum = (xnum + a[i]) * xsq; 00528 xden = (xden + b[i]) * xsq; 00529 } 00530 } else xnum = xden = 0.0; 00531 00532 temp = x * (xnum + a[3]) / (xden + b[3]); 00533 if(lower) *cum = 0.5 + temp; 00534 if(upper) *ccum = 0.5 - temp; 00535 if(log_p) { 00536 if(lower) *cum = std::log(*cum); 00537 if(upper) *ccum = std::log(*ccum); 00538 } 00539 } else if (y <= M_SQRT_32) { 00540 /* Evaluate pnorm for 0.674.. = qnorm(3/4) < |x| <= sqrt(32) 00541 * ~= 5.657 */ 00542 00543 xnum = c[8] * y; 00544 xden = y; 00545 for (i = 0; i < 7; ++i) { 00546 xnum = (xnum + c[i]) * y; 00547 xden = (xden + d[i]) * y; 00548 } 00549 temp = (xnum + c[7]) / (xden + d[7]); 00550 do_del(y); 00551 swap_tail; 00552 } else if (log_p 00553 || (lower && -37.5193 < x && x < 8.2924) 00554 || (upper && -8.2929 < x && x < 37.5193) 00555 ) { 00556 /* Evaluate pnorm for x in (-37.5, -5.657) union (5.657, 37.5) */ 00557 xsq = 1.0 / (x * x); 00558 xnum = p[5] * xsq; 00559 xden = xsq; 00560 for (i = 0; i < 4; ++i) { 00561 xnum = (xnum + p[i]) * xsq; 00562 xden = (xden + q[i]) * xsq; 00563 } 00564 temp = xsq * (xnum + p[4]) / (xden + q[4]); 00565 temp = (M_1_SQRT_2PI - temp) / y; 00566 do_del(x); 00567 swap_tail; 00568 } else { 00569 if (x > 0) { 00570 *cum = 1.; 00571 *ccum = 0.; 00572 } else { 00573 *cum = 0.; 00574 *ccum = 1.; 00575 } 00576 SCYTHE_THROW_10(scythe_convergence_error, "Did not converge"); 00577 } 00578 00579 return; 00580 } 00581 #undef SIXTEN 00582 #undef do_del 00583 #undef swap_tail 00584 00585 /* The standard normal distribution function */ 00586 double 00587 pnorm1 (double x, bool lower_tail, bool log_p) 00588 { 00589 SCYTHE_CHECK_10(! finite(x), scythe_invalid_arg, 00590 "Quantile x is inifinte (+/-Inf) or NaN"); 00591 00592 double p, cp; 00593 pnorm_both(x, &p, &cp, (lower_tail ? 0 : 1), log_p); 00594 00595 return (lower_tail ? p : cp); 00596 } 00597 } // anonymous namespace 00598 00599 /************* 00600 * Functions * 00601 *************/ 00602 00603 /* The gamma function */ 00604 00619 inline double 00620 gammafn (double x) 00621 { 00622 const double gamcs[22] = { 00623 +.8571195590989331421920062399942e-2, 00624 +.4415381324841006757191315771652e-2, 00625 +.5685043681599363378632664588789e-1, 00626 -.4219835396418560501012500186624e-2, 00627 +.1326808181212460220584006796352e-2, 00628 -.1893024529798880432523947023886e-3, 00629 +.3606925327441245256578082217225e-4, 00630 -.6056761904460864218485548290365e-5, 00631 +.1055829546302283344731823509093e-5, 00632 -.1811967365542384048291855891166e-6, 00633 +.3117724964715322277790254593169e-7, 00634 -.5354219639019687140874081024347e-8, 00635 +.9193275519859588946887786825940e-9, 00636 -.1577941280288339761767423273953e-9, 00637 +.2707980622934954543266540433089e-10, 00638 -.4646818653825730144081661058933e-11, 00639 +.7973350192007419656460767175359e-12, 00640 -.1368078209830916025799499172309e-12, 00641 +.2347319486563800657233471771688e-13, 00642 -.4027432614949066932766570534699e-14, 00643 +.6910051747372100912138336975257e-15, 00644 -.1185584500221992907052387126192e-15, 00645 }; 00646 00647 00648 double y = std::fabs(x); 00649 00650 if (y <= 10) { 00651 00652 /* Compute gamma(x) for -10 <= x <= 10 00653 * Reduce the interval and find gamma(1 + y) for 0 <= y < 1 00654 * first of all. */ 00655 00656 int n = (int) x; 00657 if (x < 0) 00658 --n; 00659 00660 y = x - n;/* n = floor(x) ==> y in [ 0, 1 ) */ 00661 --n; 00662 double value = chebyshev_eval(y * 2 - 1, gamcs, 22) + .9375; 00663 00664 if (n == 0) 00665 return value;/* x = 1.dddd = 1+y */ 00666 00667 if (n < 0) { 00668 /* compute gamma(x) for -10 <= x < 1 */ 00669 00670 /* If the argument is exactly zero or a negative integer */ 00671 /* then return NaN. */ 00672 SCYTHE_CHECK_10(x == 0 || (x < 0 && x == n + 2), 00673 scythe_range_error, "x is 0 or a negative integer"); 00674 00675 /* The answer is less than half precision */ 00676 /* because x too near a negative integer. */ 00677 SCYTHE_CHECK_10(x < -0.5 && 00678 std::fabs(x - (int)(x - 0.5) / x) < 67108864.0, 00679 scythe_precision_error, 00680 "Answer < 1/2 precision because x is too near" << 00681 " a negative integer"); 00682 00683 /* The argument is so close to 0 that the result 00684 * * would overflow. */ 00685 SCYTHE_CHECK_10(y < 2.2474362225598545e-308, scythe_range_error, 00686 "x too close to 0"); 00687 00688 n = -n; 00689 00690 for (int i = 0; i < n; i++) 00691 value /= (x + i); 00692 00693 return value; 00694 } else { 00695 /* gamma(x) for 2 <= x <= 10 */ 00696 00697 for (int i = 1; i <= n; i++) { 00698 value *= (y + i); 00699 } 00700 return value; 00701 } 00702 } else { 00703 /* gamma(x) for y = |x| > 10. */ 00704 00705 /* Overflow */ 00706 SCYTHE_CHECK_10(x > 171.61447887182298, 00707 scythe_range_error,"Overflow"); 00708 00709 /* Underflow */ 00710 SCYTHE_CHECK_10(x < -170.5674972726612, 00711 scythe_range_error, "Underflow"); 00712 00713 double value = std::exp((y - 0.5) * std::log(y) - y 00714 + M_LN_SQRT_2PI + lngammacor(y)); 00715 00716 if (x > 0) 00717 return value; 00718 00719 SCYTHE_CHECK_10(std::fabs((x - (int)(x - 0.5))/x) < 67108864.0, 00720 scythe_precision_error, 00721 "Answer < 1/2 precision because x is " << 00722 "too near a negative integer"); 00723 00724 double sinpiy = std::sin(M_PI * y); 00725 00726 /* Negative integer arg - overflow */ 00727 SCYTHE_CHECK_10(sinpiy == 0, scythe_range_error, "Overflow"); 00728 00729 return -M_PI / (y * sinpiy * value); 00730 } 00731 } 00732 00733 /* The natural log of the absolute value of the gamma function */ 00750 inline double 00751 lngammafn(double x) 00752 { 00753 SCYTHE_CHECK_10(x <= 0 && x == (int) x, scythe_range_error, 00754 "x is 0 or a negative integer"); 00755 00756 double y = std::fabs(x); 00757 00758 if (y <= 10) 00759 return std::log(std::fabs(gammafn(x))); 00760 00761 SCYTHE_CHECK_10(y > 2.5327372760800758e+305, scythe_range_error, 00762 "Overflow"); 00763 00764 if (x > 0) /* i.e. y = x > 10 */ 00765 return M_LN_SQRT_2PI + (x - 0.5) * std::log(x) - x 00766 + lngammacor(x); 00767 00768 /* else: x < -10; y = -x */ 00769 double sinpiy = std::fabs(std::sin(M_PI * y)); 00770 00771 if (sinpiy == 0) /* Negative integer argument */ 00772 throw scythe_exception("UNEXPECTED ERROR", 00773 __FILE__, __func__, __LINE__, 00774 "ERROR: Should never happen!"); 00775 00776 double ans = M_LN_SQRT_PId2 + (x - 0.5) * std::log(y) - x - std::log(sinpiy) 00777 - lngammacor(y); 00778 00779 SCYTHE_CHECK_10(std::fabs((x - (int)(x - 0.5)) * ans / x) 00780 < 1.490116119384765696e-8, scythe_precision_error, 00781 "Answer < 1/2 precision because x is " 00782 << "too near a negative integer"); 00783 00784 return ans; 00785 } 00786 00787 /* The beta function */ 00804 inline double 00805 betafn(double a, double b) 00806 { 00807 SCYTHE_CHECK_10(a <= 0 || b <= 0, scythe_invalid_arg, "a or b < 0"); 00808 00809 if (a + b < 171.61447887182298) /* ~= 171.61 for IEEE */ 00810 return gammafn(a) * gammafn(b) / gammafn(a+b); 00811 00812 double val = lnbetafn(a, b); 00813 SCYTHE_CHECK_10(val < -708.39641853226412, scythe_range_error, 00814 "Underflow"); 00815 00816 return std::exp(val); 00817 } 00818 00819 /* The natural log of the beta function */ 00837 inline double 00838 lnbetafn (double a, double b) 00839 { 00840 double p, q; 00841 00842 p = q = a; 00843 if(b < p) p = b;/* := min(a,b) */ 00844 if(b > q) q = b;/* := max(a,b) */ 00845 00846 SCYTHE_CHECK_10(p <= 0 || q <= 0,scythe_invalid_arg, "a or b <= 0"); 00847 00848 if (p >= 10) { 00849 /* p and q are big. */ 00850 double corr = lngammacor(p) + lngammacor(q) - lngammacor(p + q); 00851 return std::log(q) * -0.5 + M_LN_SQRT_2PI + corr 00852 + (p - 0.5) * std::log(p / (p + q)) + q * std::log(1 + (-p / (p + q))); 00853 } else if (q >= 10) { 00854 /* p is small, but q is big. */ 00855 double corr = lngammacor(q) - lngammacor(p + q); 00856 return lngammafn(p) + corr + p - p * std::log(p + q) 00857 + (q - 0.5) * std::log(1 + (-p / (p + q))); 00858 } 00859 00860 /* p and q are small: p <= q > 10. */ 00861 return std::log(gammafn(p) * (gammafn(q) / gammafn(p + q))); 00862 } 00863 00864 /* Compute the factorial of a non-negative integer */ 00874 inline int 00875 factorial (unsigned int n) 00876 { 00877 if (n == 0) 00878 return 1; 00879 00880 return n * factorial(n - 1); 00881 } 00882 00883 /* Compute the natural log of the factorial of a non-negative 00884 * integer 00885 */ 00895 inline double 00896 lnfactorial (unsigned int n) 00897 { 00898 double x = n+1; 00899 double cof[6] = { 00900 76.18009172947146, -86.50532032941677, 00901 24.01409824083091, -1.231739572450155, 00902 0.1208650973866179e-2, -0.5395239384953e-5 00903 }; 00904 double y = x; 00905 double tmp = x + 5.5 - (x + 0.5) * std::log(x + 5.5); 00906 double ser = 1.000000000190015; 00907 for (int j = 0; j <= 5; j++) { 00908 ser += (cof[j] / ++y); 00909 } 00910 return(std::log(2.5066282746310005 * ser / x) - tmp); 00911 } 00912 00913 /********************************* 00914 * Fully Specified Distributions * 00915 *********************************/ 00916 00917 /* These macros provide a nice shorthand for the matrix versions of 00918 * the pdf and cdf functions. 00919 */ 00920 00921 #define SCYTHE_ARGSET(...) __VA_ARGS__ 00922 00923 #define SCYTHE_DISTFUN_MATRIX(NAME, XTYPE, ARGNAMES, ...) \ 00924 template <matrix_order RO, matrix_style RS, \ 00925 matrix_order PO, matrix_style PS> \ 00926 Matrix<double, RO, RS> \ 00927 NAME (const Matrix<XTYPE, PO, PS>& X, __VA_ARGS__) \ 00928 { \ 00929 Matrix<double, RO, Concrete> ret(X.rows(), X.cols(), false); \ 00930 const_matrix_forward_iterator<XTYPE,RO,PO,PS> xit; \ 00931 const_matrix_forward_iterator<XTYPE,RO,PO,PS> xlast \ 00932 = X.template end_f<RO>(); \ 00933 typename Matrix<double,RO,Concrete>::forward_iterator rit \ 00934 = ret.begin_f(); \ 00935 for (xit = X.template begin_f<RO>(); xit != xlast; ++xit) { \ 00936 *rit = NAME (*xit, ARGNAMES); \ 00937 ++rit; \ 00938 } \ 00939 SCYTHE_VIEW_RETURN(double, RO, RS, ret) \ 00940 } \ 00941 \ 00942 template <matrix_order O, matrix_style S> \ 00943 Matrix<double, O, Concrete> \ 00944 NAME (const Matrix<XTYPE, O, S>& X, __VA_ARGS__) \ 00945 { \ 00946 return NAME <O, Concrete, O, S> (X, ARGNAMES); \ 00947 } 00948 00949 /**** The Beta Distribution ****/ 00950 00951 /* CDFs */ 00952 00981 inline double 00982 pbeta(double x, double a, double b) 00983 { 00984 SCYTHE_CHECK_10(a <= 0 || b <= 0,scythe_invalid_arg, "a or b <= 0"); 00985 00986 if (x <= 0) 00987 return 0.; 00988 if (x >= 1) 00989 return 1.; 00990 00991 return pbeta_raw(x,a,b); 00992 } 00993 00994 SCYTHE_DISTFUN_MATRIX(pbeta, double, SCYTHE_ARGSET(a, b), double a, double b) 00995 00996 /* PDFs */ 01025 inline double 01026 dbeta(double x, double a, double b) 01027 { 01028 SCYTHE_CHECK_10((x < 0.0) || (x > 1.0), scythe_invalid_arg, 01029 "x not in [0,1]"); 01030 SCYTHE_CHECK_10(a < 0.0, scythe_invalid_arg, "a < 0"); 01031 SCYTHE_CHECK_10(b < 0.0, scythe_invalid_arg, "b < 0"); 01032 01033 return (std::pow(x, (a-1.0)) * std::pow((1.0-x), (b-1.0)) ) 01034 / betafn(a,b); 01035 } 01036 01037 SCYTHE_DISTFUN_MATRIX(dbeta, double, SCYTHE_ARGSET(a, b), double a, double b) 01038 01039 /* Returns the natural log of the ordinate of the Beta density 01040 * evaluated at x with Shape1 a, and Shape2 b 01041 */ 01042 01043 01070 inline double 01071 lndbeta1(double x, double a, double b) 01072 { 01073 SCYTHE_CHECK_10((x < 0.0) || (x > 1.0), scythe_invalid_arg, 01074 "x not in [0,1]"); 01075 SCYTHE_CHECK_10(a < 0.0, scythe_invalid_arg, "a < 0"); 01076 SCYTHE_CHECK_10(b < 0.0, scythe_invalid_arg, "b < 0"); 01077 01078 return (a-1.0) * std::log(x) + (b-1) * std::log(1.0-x) 01079 - lnbetafn(a,b); 01080 } 01081 01082 SCYTHE_DISTFUN_MATRIX(lndbeta1, double, SCYTHE_ARGSET(a, b), double a, double b) 01083 01084 01085 /**** The Binomial Distribution ****/ 01086 01087 /* CDFs */ 01088 01089 01115 inline double 01116 pbinom(double x, unsigned int n, double p) 01117 { 01118 01119 SCYTHE_CHECK_10(p < 0 || p > 1, scythe_invalid_arg, "p not in [0,1]"); 01120 double X = std::floor(x); 01121 01122 if (X < 0.0) 01123 return 0; 01124 01125 if (n <= X) 01126 return 1; 01127 01128 return pbeta(1 - p, n - X, X + 1); 01129 } 01130 01131 SCYTHE_DISTFUN_MATRIX(pbinom, double, SCYTHE_ARGSET(n, p), unsigned int n, double p) 01132 01133 /* PDFs */ 01160 inline double 01161 dbinom(double x, unsigned int n, double p) 01162 { 01163 SCYTHE_CHECK_10(p < 0 || p > 1, scythe_invalid_arg, "p not in [0, 1]"); 01164 double X = std::floor(x + 0.5); 01165 return dbinom_raw(X, n, p, 1 - p); 01166 } 01167 01168 SCYTHE_DISTFUN_MATRIX(dbinom, double, SCYTHE_ARGSET(n, p), unsigned int n, double p) 01169 01170 /**** The Chi Squared Distribution ****/ 01171 01172 /* CDFs */ 01200 inline double 01201 pchisq(double x, double df) 01202 { 01203 return pgamma(x, df/2.0, 2.0); 01204 } 01205 01206 SCYTHE_DISTFUN_MATRIX(pchisq, double, df, double df) 01207 01208 /* PDFs */ 01236 inline double 01237 dchisq(double x, double df) 01238 { 01239 return dgamma(x, df / 2.0, 2.0); 01240 } 01241 01242 SCYTHE_DISTFUN_MATRIX(dchisq, double, df, double df) 01243 01244 /**** The Exponential Distribution ****/ 01245 01246 /* CDFs */ 01270 inline double 01271 pexp(double x, double scale) 01272 { 01273 SCYTHE_CHECK_10(scale <= 0, scythe_invalid_arg, "scale <= 0"); 01274 01275 if (x <= 0) 01276 return 0; 01277 01278 return (1 - std::exp(-x*scale)); 01279 } 01280 01281 SCYTHE_DISTFUN_MATRIX(pexp, double, scale, double scale) 01282 01283 /* PDFs */ 01307 inline double 01308 dexp(double x, double scale) 01309 { 01310 SCYTHE_CHECK_10(scale <= 0, scythe_invalid_arg, "scale <= 0"); 01311 01312 if (x < 0) 01313 return 0; 01314 01315 return std::exp(-x * scale) * scale; 01316 } 01317 01318 SCYTHE_DISTFUN_MATRIX(dexp, double, scale, double scale) 01319 01320 /**** The f Distribution ****/ 01321 01322 /* CDFs */ 01353 inline double 01354 pf(double x, double df1, double df2) 01355 { 01356 SCYTHE_CHECK_10(df1 <= 0 || df2 <= 0, scythe_invalid_arg, 01357 "df1 or df2 <= 0"); 01358 01359 if (x <= 0) 01360 return 0; 01361 01362 if (df2 > 4e5) 01363 return pchisq(x*df1,df1); 01364 if (df1 > 4e5) 01365 return 1-pchisq(df2/x,df2); 01366 01367 return (1-pbeta(df2 / (df2 + df1 * x), df2 / 2.0, df1 / 2.0)); 01368 } 01369 01370 SCYTHE_DISTFUN_MATRIX(pf, double, SCYTHE_ARGSET(df1, df2), double df1, double df2) 01371 01372 /* PDFs */ 01373 01374 01403 inline double 01404 df(double x, double df1, double df2) 01405 { 01406 double dens; 01407 01408 SCYTHE_CHECK_10(df1 <= 0 || df2 <= 0, scythe_invalid_arg, 01409 "df1 or df2 <= 0"); 01410 01411 if (x <= 0) 01412 return 0; 01413 01414 double f = 1 / (df2 + x * df1); 01415 double q = df2 * f; 01416 double p = x * df1 * f; 01417 01418 if (df1 >= 2) { 01419 f = df1 * q / 2; 01420 dens = dbinom_raw((df1 - 2) / 2,(df1 + df2 - 2) / 2, p, q); 01421 } else { 01422 f = (df1 * df1 * q) /(2 * p * (df1 + df2)); 01423 dens = dbinom_raw(df1 / 2,(df1 + df2)/ 2, p, q); 01424 } 01425 01426 return f*dens; 01427 } 01428 01429 SCYTHE_DISTFUN_MATRIX(df, double, SCYTHE_ARGSET(df1, df2), double df1, double df2) 01430 01431 /**** The Gamma Distribution ****/ 01432 01433 /* CDFs */ 01463 inline double 01464 pgamma (double x, double shape, double scale) 01465 { 01466 const double xbig = 1.0e+8, xlarge = 1.0e+37, 01467 alphlimit = 1000.;/* normal approx. for shape > alphlimit */ 01468 01469 int lower_tail = 1; 01470 01471 double pn1, pn2, pn3, pn4, pn5, pn6, arg, a, b, c, an, osum, sum; 01472 long n; 01473 int pearson; 01474 01475 /* check that we have valid values for x and shape */ 01476 01477 SCYTHE_CHECK_10(shape <= 0. || scale <= 0., scythe_invalid_arg, 01478 "shape or scale <= 0"); 01479 01480 x /= scale; 01481 01482 if (x <= 0.) 01483 return 0.0; 01484 01485 /* use a normal approximation if shape > alphlimit */ 01486 01487 if (shape > alphlimit) { 01488 pn1 = std::sqrt(shape) * 3. * (std::pow(x/shape, 1./3.) + 1. 01489 / (9. * shape) - 1.); 01490 return pnorm(pn1, 0., 1.); 01491 } 01492 01493 /* if x is extremely large __compared to shape__ then return 1 */ 01494 01495 if (x > xbig * shape) 01496 return 1.0; 01497 01498 if (x <= 1. || x < shape) { 01499 pearson = 1;/* use pearson's series expansion. */ 01500 arg = shape * std::log(x) - x - lngammafn(shape + 1.); 01501 c = 1.; 01502 sum = 1.; 01503 a = shape; 01504 do { 01505 a += 1.; 01506 c *= x / a; 01507 sum += c; 01508 } while (c > DBL_EPSILON); 01509 arg += std::log(sum); 01510 } 01511 else { /* x >= max( 1, shape) */ 01512 pearson = 0;/* use a continued fraction expansion */ 01513 arg = shape * std::log(x) - x - lngammafn(shape); 01514 a = 1. - shape; 01515 b = a + x + 1.; 01516 pn1 = 1.; 01517 pn2 = x; 01518 pn3 = x + 1.; 01519 pn4 = x * b; 01520 sum = pn3 / pn4; 01521 for (n = 1; ; n++) { 01522 a += 1.;/* = n+1 -shape */ 01523 b += 2.;/* = 2(n+1)-shape+x */ 01524 an = a * n; 01525 pn5 = b * pn3 - an * pn1; 01526 pn6 = b * pn4 - an * pn2; 01527 if (std::fabs(pn6) > 0.) { 01528 osum = sum; 01529 sum = pn5 / pn6; 01530 if (std::fabs(osum - sum) <= DBL_EPSILON * std::min(1., sum)) 01531 break; 01532 } 01533 pn1 = pn3; 01534 pn2 = pn4; 01535 pn3 = pn5; 01536 pn4 = pn6; 01537 if (std::fabs(pn5) >= xlarge) { 01538 /* re-scale terms in continued fraction if they are large */ 01539 pn1 /= xlarge; 01540 pn2 /= xlarge; 01541 pn3 /= xlarge; 01542 pn4 /= xlarge; 01543 } 01544 } 01545 arg += std::log(sum); 01546 } 01547 01548 lower_tail = (lower_tail == pearson); 01549 01550 sum = std::exp(arg); 01551 01552 return (lower_tail) ? sum : 1 - sum; 01553 } 01554 01555 SCYTHE_DISTFUN_MATRIX(pgamma, double, SCYTHE_ARGSET(shape, scale), double shape, double scale) 01556 01557 /* PDFs */ 01587 inline double 01588 dgamma(double x, double shape, double scale) 01589 { 01590 SCYTHE_CHECK_10(shape <= 0 || scale <= 0,scythe_invalid_arg, 01591 "shape or scale <= 0"); 01592 01593 if (x < 0) 01594 return 0.0; 01595 01596 if (x == 0) { 01597 SCYTHE_CHECK_10(shape < 1,scythe_invalid_arg, 01598 "x == 0 and shape < 1"); 01599 01600 if (shape > 1) 01601 return 0.0; 01602 01603 return 1 / scale; 01604 } 01605 01606 if (shape < 1) { 01607 double pr = dpois_raw(shape, x/scale); 01608 return pr * shape / x; 01609 } 01610 01611 /* else shape >= 1 */ 01612 double pr = dpois_raw(shape - 1, x / scale); 01613 return pr / scale; 01614 } 01615 01616 SCYTHE_DISTFUN_MATRIX(dgamma, double, SCYTHE_ARGSET(shape, scale), double shape, double scale) 01617 01618 /**** The Logistic Distribution ****/ 01619 01620 /* CDFs */ 01645 inline double 01646 plogis (double x, double location, double scale) 01647 { 01648 SCYTHE_CHECK_10(scale <= 0.0, scythe_invalid_arg, "scale <= 0"); 01649 01650 double X = (x-location) / scale; 01651 01652 X = std::exp(-X); 01653 01654 return 1 / (1+X); 01655 } 01656 01657 SCYTHE_DISTFUN_MATRIX(plogis, double, SCYTHE_ARGSET(location, scale), double location, double scale) 01658 01659 /* PDFs */ 01684 inline double 01685 dlogis(double x, double location, double scale) 01686 { 01687 SCYTHE_CHECK_10(scale <= 0.0, scythe_invalid_arg, "scale <= 0"); 01688 01689 double X = (x - location) / scale; 01690 double e = std::exp(-X); 01691 double f = 1.0 + e; 01692 01693 return e / (scale * f * f); 01694 } 01695 01696 SCYTHE_DISTFUN_MATRIX(dlogis, double, SCYTHE_ARGSET(location, scale), double location, double scale) 01697 01698 /**** The Log Normal Distribution ****/ 01699 01700 /* CDFs */ 01727 inline double 01728 plnorm (double x, double logmean, double logsd) 01729 { 01730 SCYTHE_CHECK_10(logsd <= 0, scythe_invalid_arg, "logsd <= 0"); 01731 01732 if (x > 0) 01733 return pnorm(std::log(x), logmean, logsd); 01734 01735 return 0; 01736 } 01737 01738 SCYTHE_DISTFUN_MATRIX(plnorm, double, SCYTHE_ARGSET(logmean, logsd), double logmean, double logsd) 01739 01740 /* PDFs */ 01766 inline double 01767 dlnorm(double x, double logmean, double logsd) 01768 { 01769 SCYTHE_CHECK_10(logsd <= 0, scythe_invalid_arg, "logsd <= 0"); 01770 01771 if (x == 0) 01772 return 0; 01773 01774 double y = (std::log(x) - logmean) / logsd; 01775 01776 return (1 / (std::sqrt(2 * M_PI))) * std::exp(-0.5 * y * y) / (x * logsd); 01777 } 01778 01779 SCYTHE_DISTFUN_MATRIX(dlnorm, double, SCYTHE_ARGSET(logmean, logsd), double logmean, double logsd) 01780 01781 /**** The Negative Binomial Distribution ****/ 01782 01783 /* CDFs */ 01812 inline double 01813 pnbinom(unsigned int x, double n, double p) 01814 { 01815 SCYTHE_CHECK_10(n == 0 || p <= 0 || p >= 1, scythe_invalid_arg, 01816 "n == 0 or p not in (0,1)"); 01817 01818 return pbeta(p, n, x + 1); 01819 } 01820 01821 SCYTHE_DISTFUN_MATRIX(pnbinom, unsigned int, SCYTHE_ARGSET(n, p), double n, double p) 01822 01823 /* PDFs */ 01852 inline double 01853 dnbinom(unsigned int x, double n, double p) 01854 { 01855 SCYTHE_CHECK_10(n == 0 || p <= 0 || p >= 1, scythe_invalid_arg, 01856 "n == 0 or p not in (0,1)"); 01857 01858 double prob = dbinom_raw(n, x + n, p, 1 - p); 01859 double P = (double) n / (n + x); 01860 01861 return P * prob; 01862 } 01863 01864 SCYTHE_DISTFUN_MATRIX(dnbinom, unsigned int, SCYTHE_ARGSET(n, p), double n, double p) 01865 01866 /**** The Normal Distribution ****/ 01867 01868 /* CDFs */ 01894 inline double 01895 pnorm (double x, double mean, double sd) 01896 01897 { 01898 SCYTHE_CHECK_10(sd <= 0, scythe_invalid_arg, 01899 "negative standard deviation"); 01900 01901 return pnorm1((x - mean) / sd, true, false); 01902 } 01903 01904 SCYTHE_DISTFUN_MATRIX(pnorm, double, SCYTHE_ARGSET(mean, sd), double mean, double sd) 01905 01906 01907 /* PDFs */ 01932 inline double 01933 dnorm(double x, double mean, double sd) 01934 { 01935 SCYTHE_CHECK_10(sd <= 0, scythe_invalid_arg, 01936 "negative standard deviation"); 01937 01938 double X = (x - mean) / sd; 01939 01940 return (M_1_SQRT_2PI * std::exp(-0.5 * X * X) / sd); 01941 } 01942 01943 SCYTHE_DISTFUN_MATRIX(dnorm, double, SCYTHE_ARGSET(mean, sd), double mean, double sd) 01944 01945 01946 /* Return the natural log of the normal PDF */ 01972 inline double 01973 lndnorm (double x, double mean, double sd) 01974 { 01975 SCYTHE_CHECK_10(sd <= 0, scythe_invalid_arg, 01976 "negative standard deviation"); 01977 01978 double X = (x - mean) / sd; 01979 01980 return -(M_LN_SQRT_2PI + 0.5 * X * X + std::log(sd)); 01981 } 01982 01983 SCYTHE_DISTFUN_MATRIX(lndnorm, double, SCYTHE_ARGSET(mean, sd), double mean, double sd) 01984 01985 /* Quantile functions */ 02008 inline double 02009 qnorm1 (double in_p) 02010 { 02011 double p0 = -0.322232431088; 02012 double q0 = 0.0993484626060; 02013 double p1 = -1.0; 02014 double q1 = 0.588581570495; 02015 double p2 = -0.342242088547; 02016 double q2 = 0.531103462366; 02017 double p3 = -0.0204231210245; 02018 double q3 = 0.103537752850; 02019 double p4 = -0.453642210148e-4; 02020 double q4 = 0.38560700634e-2; 02021 double xp = 0.0; 02022 double p = in_p; 02023 02024 if (p > 0.5) 02025 p = 1 - p; 02026 02027 SCYTHE_CHECK_10(p < 10e-20, scythe_range_error, 02028 "p outside accuracy limit"); 02029 02030 if (p == 0.5) 02031 return xp; 02032 02033 double y = std::sqrt (std::log (1.0 / std::pow (p, 2))); 02034 xp = y + ((((y * p4 + p3) * y + p2) * y + p1) * y + p0) / 02035 ((((y * q4 + q3) * y + q2) * y + q1) * y + q0); 02036 02037 if (in_p < 0.5) 02038 xp = -1 * xp; 02039 02040 return xp; 02041 } 02042 02043 SCYTHE_DISTFUN_MATRIX(qnorm1, double, in_p, double in_p) 02044 02045 /**** The Poisson Distribution ****/ 02046 02047 /* CDFs */ 02074 inline double 02075 ppois(unsigned int x, double lambda) 02076 { 02077 SCYTHE_CHECK_10(lambda<=0.0, scythe_invalid_arg, "lambda <= 0"); 02078 02079 if (lambda == 1) 02080 return 1; 02081 02082 return 1 - pgamma(lambda, x + 1, 1.0); 02083 } 02084 02085 SCYTHE_DISTFUN_MATRIX(ppois, unsigned int, lambda, double lambda) 02086 02087 /* PDFs */ 02111 inline double 02112 dpois(unsigned int x, double lambda) 02113 { 02114 SCYTHE_CHECK_10(lambda<=0.0, scythe_invalid_arg, "lambda <= 0"); 02115 02116 // compute log(x!) 02117 double xx = x+1; 02118 double cof[6] = { 02119 76.18009172947146, -86.50532032941677, 02120 24.01409824083091, -1.231739572450155, 02121 0.1208650973866179e-2, -0.5395239384953e-5 02122 }; 02123 double y = xx; 02124 double tmp = xx + 5.5 - (xx + 0.5) * std::log(xx + 5.5); 02125 double ser = 1.000000000190015; 02126 for (int j = 0; j <= 5; j++) { 02127 ser += (cof[j] / ++y); 02128 } 02129 double lnfactx = std::log(2.5066282746310005 * ser / xx) - tmp; 02130 02131 return (std::exp( -1*lnfactx + x * std::log(lambda) - lambda)); 02132 } 02133 02134 SCYTHE_DISTFUN_MATRIX(dpois, unsigned int, lambda, double lambda) 02135 02136 /**** The t Distribution ****/ 02137 02138 /* CDFs */ 02165 inline double 02166 pt(double x, double n) 02167 { 02168 double val; 02169 02170 SCYTHE_CHECK_10(n <= 0, scythe_invalid_arg, "n <= 0"); 02171 02172 if (n > 4e5) { 02173 val = 1/(4*n); 02174 return pnorm1(x * (1 - val) / ::sqrt(1 + x * x * 2. * val), 02175 true, false); 02176 } 02177 02178 val = pbeta(n / (n + x * x), n / 2.0, 0.5); 02179 02180 val /= 2; 02181 02182 if (x <= 0) 02183 return val; 02184 else 02185 return 1 - val; 02186 } 02187 02188 SCYTHE_DISTFUN_MATRIX(pt, double, n, double n) 02189 02190 /* PDFs */ 02216 inline double 02217 dt(double x, double n) 02218 { 02219 double u; 02220 02221 SCYTHE_CHECK_10(n <= 0, scythe_invalid_arg, "n <= 0"); 02222 02223 double t = -bd0(n/2., (n + 1) / 2.) 02224 + stirlerr((n + 1) / 2.) 02225 - stirlerr(n / 2.); 02226 if(x*x > 0.2*n) 02227 u = std::log(1+x*x/n)*n/2; 02228 else 02229 u = -bd0(n/2., (n+x*x)/2.) + x*x/2; 02230 02231 return std::exp(t-u)/std::sqrt(2*M_PI*(1+x*x/n)); 02232 } 02233 02234 SCYTHE_DISTFUN_MATRIX(dt, double, n, double n) 02235 02236 /* Returns the univariate Student-t density evaluated at x 02237 * with mean mu, scale sigma^2, and nu degrees of freedom. 02238 * 02239 * TODO: Do we want a pt for this distribution? 02240 */ 02241 02242 02270 inline double 02271 dt1(double x, double mu, double sigma2, double nu) 02272 { 02273 double logdens = lngammafn((nu + 1.0) /2.0) 02274 - std::log(std::sqrt(nu * M_PI)) 02275 - lngammafn(nu / 2.0) - std::log(std::sqrt(sigma2)) 02276 - (nu + 1.0) / 2.0 * std::log(1.0 + (std::pow((x - mu), 2.0)) 02277 / (nu * sigma2)); 02278 02279 return(std::exp(logdens)); 02280 } 02281 02282 SCYTHE_DISTFUN_MATRIX(dt1, double, SCYTHE_ARGSET(mu, sigma2, nu), double mu, double sigma2, double nu) 02283 02284 /* Returns the natural log of the univariate Student-t density 02285 * evaluated at x with mean mu, scale sigma^2, and nu 02286 * degrees of freedom 02287 */ 02288 02289 02318 inline double 02319 lndt1(double x, double mu, double sigma2, double nu) 02320 { 02321 double logdens = lngammafn((nu+1.0)/2.0) 02322 - std::log(std::sqrt(nu*M_PI)) 02323 - lngammafn(nu/2.0) - std::log(std::sqrt(sigma2)) 02324 - (nu+1.0)/2.0 * std::log(1.0 + (std::pow((x-mu),2.0)) 02325 /(nu * sigma2)); 02326 02327 return(logdens); 02328 } 02329 02330 SCYTHE_DISTFUN_MATRIX(lndt1, double, SCYTHE_ARGSET(mu, sigma2, nu), double mu, double sigma2, double nu) 02331 02332 /**** The Uniform Distribution ****/ 02333 02334 /* CDFs */ 02359 inline double 02360 punif(double x, double a, double b) 02361 { 02362 SCYTHE_CHECK_10(b <= a, scythe_invalid_arg, "b <= a"); 02363 02364 if (x <= a) 02365 return 0.0; 02366 02367 if (x >= b) 02368 return 1.0; 02369 02370 return (x - a) / (b - a); 02371 } 02372 02373 SCYTHE_DISTFUN_MATRIX(punif, double, SCYTHE_ARGSET(a, b), double a, double b) 02374 02375 /* PDFs */ 02400 inline double 02401 dunif(double x, double a, double b) 02402 { 02403 SCYTHE_CHECK_10(b <= a, scythe_invalid_arg, "b <= a"); 02404 02405 if (a <= x && x <= b) 02406 return 1.0 / (b - a); 02407 02408 return 0.0; 02409 } 02410 02411 SCYTHE_DISTFUN_MATRIX(dunif, double, SCYTHE_ARGSET(a, b), double a, double b) 02412 02413 /**** The Weibull Distribution ****/ 02414 02415 /* CDFs */ 02440 inline double 02441 pweibull(double x, double shape, double scale) 02442 { 02443 SCYTHE_CHECK_10(shape <= 0 || scale <= 0, scythe_invalid_arg, 02444 "shape or scale <= 0"); 02445 02446 if (x <= 0) 02447 return 0.0; 02448 02449 return 1 - std::exp(-std::pow(x / scale, shape)); 02450 } 02451 02452 SCYTHE_DISTFUN_MATRIX(pweibull, double, SCYTHE_ARGSET(shape, scale), double shape, double scale) 02453 02454 /* PDFs */ 02479 inline double 02480 dweibull(double x, double shape, double scale) 02481 { 02482 SCYTHE_CHECK_10(shape <= 0 || scale <= 0, scythe_invalid_arg, 02483 "shape or scale <= 0"); 02484 02485 if (x < 0) 02486 return 0.; 02487 02488 double tmp1 = std::pow(x / scale, shape - 1); 02489 double tmp2 = tmp1*(x / scale); 02490 02491 return shape * tmp1 * std::exp(-tmp2) / scale; 02492 } 02493 02494 SCYTHE_DISTFUN_MATRIX(dweibull, double, SCYTHE_ARGSET(shape, scale), double shape, double scale) 02495 02496 /* Multivariate Normal */ 02497 02498 // TODO: distribution function. Plain old (non-logged) dmvnorm. 02499 02500 02518 template <matrix_order O1, matrix_style S1, 02519 matrix_order O2, matrix_style S2, 02520 matrix_order O3, matrix_style S3> 02521 double lndmvn (const Matrix<double, O1, S1>& x, 02522 const Matrix<double, O2, S2>& mu, 02523 const Matrix<double, O3, S3>& Sigma) 02524 { 02525 SCYTHE_CHECK_10(! x.isColVector(), scythe_dimension_error, 02526 "x is not a column vector"); 02527 SCYTHE_CHECK_10(! mu.isColVector(), scythe_dimension_error, 02528 "mu is not a column vector"); 02529 SCYTHE_CHECK_10(! Sigma.isSquare(), scythe_dimension_error, 02530 "Sigma is not square"); 02531 SCYTHE_CHECK_10(mu.rows()!=Sigma.rows() || x.rows()!=Sigma.rows(), 02532 scythe_conformation_error, 02533 "mu, x and Sigma have mismatched row lengths") 02534 int k = (int) mu.rows(); 02535 return ( (-k/2.0)*std::log(2*M_PI) -0.5 * std::log(det(Sigma)) 02536 -0.5 * (t(x - mu)) * invpd(Sigma) * (x-mu) )[0]; 02537 } 02538 02539 } // end namespace scythe 02540 02541 02542 #endif /* SCYTHE_DISTRIBUTIONS_H */
1.7.6.1