gamma.js 1.1 KB

12345678910111213141516171819202122232425262728293031323334
  1. import defaultSource from "./defaultSource.js";
  2. import normal from "./normal.js";
  3. export default (function sourceRandomGamma(source) {
  4. var randomNormal = normal.source(source)();
  5. function randomGamma(k, theta) {
  6. if ((k = +k) < 0) throw new RangeError("invalid k");
  7. // degenerate distribution if k === 0
  8. if (k === 0) return () => 0;
  9. theta = theta == null ? 1 : +theta;
  10. // exponential distribution if k === 1
  11. if (k === 1) return () => -Math.log1p(-source()) * theta;
  12. var d = (k < 1 ? k + 1 : k) - 1 / 3,
  13. c = 1 / (3 * Math.sqrt(d)),
  14. multiplier = k < 1 ? () => Math.pow(source(), 1 / k) : () => 1;
  15. return function() {
  16. do {
  17. do {
  18. var x = randomNormal(),
  19. v = 1 + c * x;
  20. } while (v <= 0);
  21. v *= v * v;
  22. var u = 1 - source();
  23. } while (u >= 1 - 0.0331 * x * x * x * x && Math.log(u) >= 0.5 * x * x + d * (1 - v + Math.log(v)));
  24. return d * v * multiplier() * theta;
  25. };
  26. }
  27. randomGamma.source = sourceRandomGamma;
  28. return randomGamma;
  29. })(defaultSource);