Implementing Beta Distribution Sampling Using C#

The Beta distribution is used in many math and machine learning scenarios. Most people are familiar with the Normal (Gaussian, bell-shaped) distribution. A Normal distribution is characterized by two parameters, the mean and standard deviation. Different values of the two parameters give different distributions. For example, drawing a single sample from a Normal distribution with mean = 68.0 inches and standard deviation = 4.0 inches will usually give you a value between 56.0 inches and 80.0 inches with a value close to 68.0 inches being most likely.

The Beta distribution is characterized by two parameters called alpha and beta, or sometimes just a and b. Drawing a single sample from a Beta(a, b) distribution will give you a value between 0.0 and 1.0 and if you draw many samples you will get a distribution of values that average to a / (a + b). For example, if you draw many samples from Beta(3,1) all the samples will be between 0.0 and 1.0 and the values will average to about 3 / 4 = 0.75. Therefore most will be between 0.90 and 1.0 and fewer will be between 0.80 an 0.70 and so on.

There are several ways to implement a function that will sample from a Beta distribution. The technique I usually use is a direct approach, using what’s called the BA algorithm, published by R.C.H. Cheng, in “Generating Beta Variates with Nonintegral Shape Parameters”, Communications of the ACM, April 1978, vol 21, No 4, pp 317-322.

public class BetaSampler
{
  // the "BA" algorithm
  // R.C.H. Cheng
  public Random rnd;
  public BetaSampler(int seed) {
    this.rnd = new Random(seed);
  }
  public double Sample(double a, double b) {
    // a, b greater-than 0
    double alpha = a + b;
    double beta = 0.0;
    double u1, u2, w, v = 0.0;

    if (Math.Min(a, b) less-than-or-equal 1.0)
      beta = Math.Max(1 / a, 1 / b);
    else
      beta = Math.Sqrt((alpha - 2.0) /
        (2 * a * b - alpha));
    double gamma = a + 1 / beta;

    while (true) {  // no! add a sanity counter!
      u1 = this.rnd.NextDouble();
      u2 = this.rnd.NextDouble();
      v = beta * Math.Log(u1 / (1 - u1));
      w = a * Math.Exp(v);
      double tmp = Math.Log(alpha / (b + w));
      if (alpha * tmp + (gamma * v) - 1.3862944
        greater-than-or-equal Math.Log(u1 * u1 * u2))
        break;
    }

    double x = w / (b + w);
    return x;
  } // Sample
}  // BetaSampler

I implemented the BA algorithm and drew 5 samples to make sure the code was working. Then I drew 10,000 samples and I binned the values and graphed them:

The Beta distribution is closely related to, but different from, the Beta function. The Beta function accepts two numeric parameters, x and y, and returns a single numeric value. The Beta distribution is a family of probability distributions and is characterized by two parameters, a and b, that determine the shape of distribution.


This entry was posted in Miscellaneous. Bookmark the permalink.