13

I need to generate random numbers from JavaScript within the Beta probability distribution. I've Googled but can't find any libraries that appear to support this.

Can anyone suggest where I can find a library or code snippet that will do this?

3
  • The Wikipedia article explains how to do just that. But you'll have to generate Gamma-distributed RVs to do so
    – Blender
    Mar 6, 2012 at 19:13
  • I was hoping for a library, there are certainly examples in other languages
    – sanity
    Mar 6, 2012 at 19:16
  • If you can't come up with a solution, at least try converting the other examples into JS. This isn't that complicated.
    – Blender
    Mar 6, 2012 at 19:17

5 Answers 5

12

The jStat library has functions to sample from a beta distribution, as well as many other distributions.

var random_num = jStat.beta.sample( alpha, beta );
11

My translation. It's pretty much word for word, so it's probably not the most idiomatic javascript.

// javascript shim for Python's built-in 'sum'
function sum(nums) {
  var accumulator = 0;
  for (var i = 0, l = nums.length; i < l; i++)
    accumulator += nums[i];
  return accumulator;
}

// In case you were wondering, the nice functional version is slower.
// function sum_slow(nums) {
//   return nums.reduce(function(a, b) { return a + b; }, 0);
// }
// var tenmil = _.range(1e7); sum(tenmil); sum_slow(tenmil);

// like betavariate, but more like R's name
function rbeta(alpha, beta) {
  var alpha_gamma = rgamma(alpha, 1);
  return alpha_gamma / (alpha_gamma + rgamma(beta, 1));
}

// From Python source, so I guess it's PSF Licensed
var SG_MAGICCONST = 1 + Math.log(4.5);
var LOG4 = Math.log(4.0);

function rgamma(alpha, beta) {
  // does not check that alpha > 0 && beta > 0
  if (alpha > 1) {
    // Uses R.C.H. Cheng, "The generation of Gamma variables with non-integral
    // shape parameters", Applied Statistics, (1977), 26, No. 1, p71-74
    var ainv = Math.sqrt(2.0 * alpha - 1.0);
    var bbb = alpha - LOG4;
    var ccc = alpha + ainv;

    while (true) {
      var u1 = Math.random();
      if (!((1e-7 < u1) && (u1 < 0.9999999))) {
        continue;
      }
      var u2 = 1.0 - Math.random();
      v = Math.log(u1/(1.0-u1))/ainv;
      x = alpha*Math.exp(v);
      var z = u1*u1*u2;
      var r = bbb+ccc*v-x;
      if (r + SG_MAGICCONST - 4.5*z >= 0.0 || r >= Math.log(z)) {
        return x * beta;
      }
    }
  }
  else if (alpha == 1.0) {
    var u = Math.random();
    while (u <= 1e-7) {
      u = Math.random();
    }
    return -Math.log(u) * beta;
  }
  else { // 0 < alpha < 1
    // Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle
    while (true) {
      var u3 = Math.random();
      var b = (Math.E + alpha)/Math.E;
      var p = b*u3;
      if (p <= 1.0) {
        x = Math.pow(p, (1.0/alpha));
      }
      else {
        x = -Math.log((b-p)/alpha);
      }
      var u4 = Math.random();
      if (p > 1.0) {
        if (u4 <= Math.pow(x, (alpha - 1.0))) {
          break;
        }
      }
      else if (u4 <= Math.exp(-x)) {
        break;
      }
    }
    return x * beta;
  }
}

Partially testable with means, which are easily calculated:

function testbeta(a, b, N) {
  var sample_mean = sum(_.range(N).map(function() { return rbeta(a, b); })) / N;
  var analytic_mean = a / (a + b);
  console.log(sample_mean, "~", analytic_mean);
}
testbeta(5, 1, 100000);
1
  • P.S. Thanks to @Blender below for digging up the Python source code.
    – chbrown
    Nov 26, 2012 at 18:24
3

You can convert this Python code to JS:

SG_MAGICCONST = 1.0 + _log(4.5)
LOG4 = log(4.0)

def gamma(z, sqrt2pi=(2.0*pi)**0.5):
  # Reflection to right half of complex plane
  if z < 0.5:
      return pi / sin(pi*z) / gamma(1.0-z)
  # Lanczos approximation with g=7
  az = z + (7.0 - 0.5)
  return az ** (z-0.5) / exp(az) * sqrt2pi * fsum([
    0.9999999999995183,
    676.5203681218835 / z,
    -1259.139216722289 / (z+1.0),
    771.3234287757674 / (z+2.0),
    -176.6150291498386 / (z+3.0),
    12.50734324009056 / (z+4.0),
    -0.1385710331296526 / (z+5.0),
    0.9934937113930748e-05 / (z+6.0),
    0.1659470187408462e-06 / (z+7.0),
  ])



def gammavariate(self, alpha, beta):
  """Gamma distribution.  Not the gamma function!

  Conditions on the parameters are alpha > 0 and beta > 0.

  The probability distribution function is:

        x ** (alpha - 1) * math.exp(-x / beta)
    pdf(x) =  --------------------------------------
          math.gamma(alpha) * beta ** alpha

  """

  # alpha > 0, beta > 0, mean is alpha*beta, variance is alpha*beta**2

  # Warning: a few older sources define the gamma distribution in terms
  # of alpha > -1.0
  if alpha <= 0.0 or beta <= 0.0:
    raise ValueError, 'gammavariate: alpha and beta must be > 0.0'

  random = self.random
  if alpha > 1.0:

    # Uses R.C.H. Cheng, "The generation of Gamma
    # variables with non-integral shape parameters",
    # Applied Statistics, (1977), 26, No. 1, p71-74

    ainv = _sqrt(2.0 * alpha - 1.0)
    bbb = alpha - LOG4
    ccc = alpha + ainv

    while 1:
      u1 = random()
      if not 1e-7 < u1 < .9999999:
        continue
      u2 = 1.0 - random()
      v = _log(u1/(1.0-u1))/ainv
      x = alpha*_exp(v)
      z = u1*u1*u2
      r = bbb+ccc*v-x
      if r + SG_MAGICCONST - 4.5*z >= 0.0 or r >= _log(z):
        return x * beta

  elif alpha == 1.0:
    # expovariate(1)
    u = random()
    while u <= 1e-7:
      u = random()
    return -_log(u) * beta

  else:   # alpha is between 0 and 1 (exclusive)

    # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle

    while 1:
      u = random()
      b = (_e + alpha)/_e
      p = b*u
      if p <= 1.0:
        x = p ** (1.0/alpha)
      else:
        x = -_log((b-p)/alpha)
      u1 = random()
      if p > 1.0:
        if u1 <= x ** (alpha - 1.0):
          break
      elif u1 <= _exp(-x):
        break
    return x * beta



def betavariate(alpha, beta):
  if y == 0:
    return 0.0
  else:
    return y / (y + gammavariate(beta, 1.0))

It's directly from the Python source code (with slight modifications), but it should be easy to convert.

2

See stdlib, which includes seedable PRNGs for many distributions, including the beta distribution. For example, within the stdlib development environment,

var beta = require( '@stdlib/random/base/beta' );

var r = beta( 2.0, 5.0 );
// returns <number>

Otherwise, see the source code which is released under an Apache license.

2
  • Nice work, this will come in handy for much more than beta distributions. Aug 26, 2020 at 21:45
  • 1
    @PhilipKirkbride Thanks for updating! I also updated the link to point to the most recent URL on GitHub.
    – kgryte
    Aug 27, 2020 at 22:15
-2

http://mbostock.github.com/d3/ex/stack.html

You have an example of generating stacked bar chart layers using gamma distributions. Quiet close to what you are looking for I guess.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.