Beta Distribution in Javascript

An implementation of the beta distribution probability density function in Javascript.
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:

$$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$.