Beta Distribution in Javascript

An implementation of the beta distribution probability density function in Javascript. <br> This implementation overcomes the problem of large numbers being generated by the Beta function which can cause JS to return inf values.

Beta Distribution

X \sim Beta(\alpha, \beta)

Probability Density Function

\frac{x^{\alpha-1}(1-x)^{\beta-1}}{B(\alpha,\beta)},  x \in [0,1]

where $B(\alpha,\beta)$ is the Beta function: <br>

B(\alpha,\beta) = \frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha+\beta)}

B(\alpha,\beta) =\frac{(x-1)!(y-1)!}{(x+y-1)!}

A naive implementation of the beta distribution using the equations above will work as intended for smaller values of $\alpha$ and $\beta$ but will fail for larger values due to its use of factorials:

function naiveBetaPDF(x, a, b) {
    // Naive implementation of the beta pdf function
    // Using factorials
    return Math.pow(x, a-1)*Math.pow(1-x, b-1)/naiveBetaFunc(a,b)
}
function naiveBetaFunc(a ,b) {
    // Naive implementation of the beta function
    // using factorials
    return factorial(a-1)*factorial(b-1)/factorial(a+b-1)
}
function factorial(x) { 
    if (x === 0) {
        return 1;
    }
    return x * factorial(x-1);     
}

console.log(naiveBetaPDF(x=0.5, a=10, b=10))     // 3.5239410400390625
console.log(naiveBetaPDF(x=0.5, a=100, b=100))   // NaN
console.log(naiveBetaPDF(x=0.5, a=1000, b=1000)) // NaN

The function fails at $\alpha=100$ and $\beta=100$ because the naive beta function returns an Infinity type value on its numerator. This happens because in Javascript, the largest possible numeric value is 1.79E+308 or $2^{308}$. Values larger than that are represented as Infinity. Therefore, $99!99!$ cannot be computed. This is quite limiting for the function and another approach is required.

To overcome this problem, we can use the log beta function instead to calculate $B(\alpha, \beta)$

log(Beta(\alpha, \beta)) = log((\alpha-1)!) + ln((\beta-1)!) - ln((\alpha+\beta-1)!)

log(Beta(\alpha, \beta)) = \sum_{i=0}^{\alpha-2} log(\alpha-1-i) + \sum_{i=0}^{\beta-2} log(\beta-1-i) + \sum_{i=0}^{\alpha+\beta-2} log(\alpha+\beta-1-i)

Therefore, the log beta pdf function becomes:

(\alpha-1)log(x) + (\beta-1)log(1-x) - \sum_{i=0}^{\alpha-2} log(\alpha-1-i) - \sum_{i=0}^{\beta-2} log(\beta-1-i) - \sum_{i=0}^{\alpha+\beta-2} log(\alpha+\beta-1-i)

In code:

function betaPDF(x, a, b) {
    // Beta probability density function impementation
    // using logarithms, no factorials involved.
    // Overcomes the problem with large integers
    return Math.exp(lnBetaPDF(x, a, b))
}
function lnBetaPDF(x, a, b) {
        // Log of the Beta Probability Density Function
    return ((a-1)*Math.log(x) + (b-1)*Math.log(1-x)) - lnBetaFunc(a,b)
}
function lnBetaFunc(a, b) {
		// Log Beta Function
	  // ln(Beta(x,y))
    foo = 0.0;

    for (i=0; i<a-2; i++) {
        foo += Math.log(a-1-i);
    }
    for (i=0; i<b-2; i++) {
        foo += Math.log(b-1-i);
    }
    for (i=0; i<a+b-2; i++) {
        foo -= Math.log(a+b-1-i);
    }
    return foo
}

console.log(betaPDF(x=0.5, a=10, b=10))       // 3.5239410400390625
console.log(betaPDF(x=0.5, a=100, b=100))     // 11.269695801846181
console.log(betaPDF(x=0.5, a=1000, b=1000))   // 35.67802229201808
console.log(betaPDF(x=0.5, a=10000, b=10000)) // 112.83650628497722

Now the function is able to handle larger values of $\alpha$ and $\beta$.