Description
LogBeta computes the log of the beta function using a log of a computation using the Gamma function when the arguments are small (a < 1). The code includes a comment that this is more accurate than using LogGamma:
// The original NSWC implementation was // LogGamma.value(a) + (LogGamma.value(b) - LogGamma.value(a + b)); // but the following command turned out to be more accurate. return Math.log(Gamma.value(a) * Gamma.value(b) / Gamma.value(a + b));
However Gamma(z) -> 1/z when z -> 0. Thus this can overflow, e.g. this test fails:
@Test void testTinyAB() { double a = 1.0e-155; double b = 1.5e-155; double beta = 2.5 / 1.5 * 1e155; Assertions.assertEquals(Math.log(beta), LogBeta.value(a, b)); }
This can be fixed by checking for overflow and reverting to the original NSWC implementation:
final double beta = Gamma.value(a) * Gamma.value(b) / Gamma.value(a + b); if (Double.isFinite(beta)) { return Math.log(beta); } return LogGamma.value(a) + (LogGamma.value(b) - LogGamma.value(a + b));