stat.js

// https://github.com/numpy/numpy/blob/master/numpy/random/mtrand/distributions.c
// https://github.com/numpy/numpy/blob/master/numpy/random/mtrand/randomkit.c

function DistributionSampling(seed) {
  // xorshift: https://sbfl.net/blog/2017/06/01/javascript-reproducible-random/
  var state = {
      x: 123456789,
      y: 362436069,
      z: 521288629,
      w: seed || 88675123,
      hasGauss: false,
      gauss: 0
    };

  var self = this;

  function setSeed(s) {
    state.w = s;
  }

  function next() {
    var t = state.x ^ (state.x << 11);
    state.x = state.y; state.y = state.z; state.z = state.w;
    state.w = (state.w ^ (state.w >>> 19)) ^ (t ^ (t >>> 8));
    return state.w;
  }

  function random() {
    var y = next();
    return (y % 0xFFFFFFF) / 0xFFFFFFF;
  }

  function gauss() {
    if (state.hasGauss) {
      var tmp = state.gauss;
      state.gauss = 0;
      state.hasGauss = false;
      return tmp;
    } else {
      var f, x1, x2, r2;
      do {
        x1 = 2.0 * random() - 1;
        x2 = 2.0 * random() - 1;
        r2 = x1*x1 + x2*x2;
      } while (r2 >= 1.0 || r2 == 0.0);

      f = Math.sqrt(-2.0*Math.log(r2)/r2);
      state.gauss = f * x1;
      state.hasGauss = true;
      return f * x2;
    }
  }

  function standardExponential() {
    return -Math.log(1.0 - random());
  }

  function standardGamma(shape) {
    var b, c;
    var U, V, X, Y;
    if (shape < 1.0) {  // not implemented
      throw "can't calculate shape <= 1.0";
    }
    if (shape == 1.0) {
      return standardExponential();
    }

    b = shape - 1.0/3.0;
    c = 1.0/Math.sqrt(9*b);
    while(true) {
      do {
        X = gauss();
        V = 1.0 + c*X;
      } while (V <= 0.0);

      V = V*V*V;
      U = random();
      if (U < 1.0 - 0.0331*(X*X)*(X*X)) {
        return b*V;
      }
      if (Math.log(U) < 0.5*X*X + b*(1.0 - V+Math.log(V))) {
        return b*V;
      }
    }
  }

  function gamma(shape, scale) {
    return scale * standardGamma(shape);
  }

  function beta(a, b) {
    if (a < 1.0 && b < 1.0) {  // not implemented
      throw "can't calc beta sampling when a <= 1.0 && b <= 1.0";
    }

    var ga, gb;
    ga = standardGamma(a);
    gb = standardGamma(b);
    return ga / (ga + gb);
  }

  function betaSampling(a, b, n) {
    n = n || 1;
    var ret = [];
    for (var i=0; i<n; i++) {
      ret.push(beta(a, b));
    }
    return ret;
  }

  self.beta = betaSampling;
  return self;
}