probability.js 3.2 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192
  1. "use strict";
  2. Object.defineProperty(exports, "__esModule", {
  3. value: true
  4. });
  5. exports.gammaG = void 0;
  6. exports.gammaNumber = gammaNumber;
  7. exports.lgammaN = exports.lgammaG = exports.gammaP = void 0;
  8. exports.lgammaNumber = lgammaNumber;
  9. exports.lnSqrt2PI = exports.lgammaSeries = void 0;
  10. var _number = require("../../utils/number.js");
  11. var _product = require("../../utils/product.js");
  12. /* eslint-disable no-loss-of-precision */
  13. function gammaNumber(n) {
  14. var x;
  15. if ((0, _number.isInteger)(n)) {
  16. if (n <= 0) {
  17. return isFinite(n) ? Infinity : NaN;
  18. }
  19. if (n > 171) {
  20. return Infinity; // Will overflow
  21. }
  22. return (0, _product.product)(1, n - 1);
  23. }
  24. if (n < 0.5) {
  25. return Math.PI / (Math.sin(Math.PI * n) * gammaNumber(1 - n));
  26. }
  27. if (n >= 171.35) {
  28. return Infinity; // will overflow
  29. }
  30. if (n > 85.0) {
  31. // Extended Stirling Approx
  32. var twoN = n * n;
  33. var threeN = twoN * n;
  34. var fourN = threeN * n;
  35. var fiveN = fourN * n;
  36. return Math.sqrt(2 * Math.PI / n) * Math.pow(n / Math.E, n) * (1 + 1 / (12 * n) + 1 / (288 * twoN) - 139 / (51840 * threeN) - 571 / (2488320 * fourN) + 163879 / (209018880 * fiveN) + 5246819 / (75246796800 * fiveN * n));
  37. }
  38. --n;
  39. x = gammaP[0];
  40. for (var i = 1; i < gammaP.length; ++i) {
  41. x += gammaP[i] / (n + i);
  42. }
  43. var t = n + gammaG + 0.5;
  44. return Math.sqrt(2 * Math.PI) * Math.pow(t, n + 0.5) * Math.exp(-t) * x;
  45. }
  46. gammaNumber.signature = 'number';
  47. // TODO: comment on the variables g and p
  48. var gammaG = 4.7421875;
  49. exports.gammaG = gammaG;
  50. var gammaP = [0.99999999999999709182, 57.156235665862923517, -59.597960355475491248, 14.136097974741747174, -0.49191381609762019978, 0.33994649984811888699e-4, 0.46523628927048575665e-4, -0.98374475304879564677e-4, 0.15808870322491248884e-3, -0.21026444172410488319e-3, 0.21743961811521264320e-3, -0.16431810653676389022e-3, 0.84418223983852743293e-4, -0.26190838401581408670e-4, 0.36899182659531622704e-5];
  51. // lgamma implementation ref: https://mrob.com/pub/ries/lanczos-gamma.html#code
  52. // log(2 * pi) / 2
  53. exports.gammaP = gammaP;
  54. var lnSqrt2PI = 0.91893853320467274178;
  55. exports.lnSqrt2PI = lnSqrt2PI;
  56. var lgammaG = 5; // Lanczos parameter "g"
  57. exports.lgammaG = lgammaG;
  58. var lgammaN = 7; // Range of coefficients "n"
  59. exports.lgammaN = lgammaN;
  60. var lgammaSeries = [1.000000000190015, 76.18009172947146, -86.50532032941677, 24.01409824083091, -1.231739572450155, 0.1208650973866179e-2, -0.5395239384953e-5];
  61. exports.lgammaSeries = lgammaSeries;
  62. function lgammaNumber(n) {
  63. if (n < 0) return NaN;
  64. if (n === 0) return Infinity;
  65. if (!isFinite(n)) return n;
  66. if (n < 0.5) {
  67. // Use Euler's reflection formula:
  68. // gamma(z) = PI / (sin(PI * z) * gamma(1 - z))
  69. return Math.log(Math.PI / Math.sin(Math.PI * n)) - lgammaNumber(1 - n);
  70. }
  71. // Compute the logarithm of the Gamma function using the Lanczos method
  72. n = n - 1;
  73. var base = n + lgammaG + 0.5; // Base of the Lanczos exponential
  74. var sum = lgammaSeries[0];
  75. // We start with the terms that have the smallest coefficients and largest denominator
  76. for (var i = lgammaN - 1; i >= 1; i--) {
  77. sum += lgammaSeries[i] / (n + i);
  78. }
  79. return lnSqrt2PI + (n + 0.5) * Math.log(base) - base + Math.log(sum);
  80. }
  81. lgammaNumber.signature = 'number';