Gamma { public static double beta(double alpha, double beta) { double y; if (alpha < 40 && beta < 40) { y = gamma(alpha + beta); if (y == 0.0) { return 1.0; } if (alpha > beta) { y = gamma(alpha) / y; y *= gamma(beta); } else { y = gamma(beta) / y; y *= gamma(alpha); } } else { y = Math.exp(logGamma(alpha) + logGamma(beta) - logGamma(alpha + beta)); } return y; } private  Gamma(); static double beta(double alpha, double beta); static double gamma(double x); static double incompleteBeta(double alpha, double beta, double xx); static double incompleteGamma(double alpha, double x); static double incompleteGammaComplement(double alpha, double x); static double logGamma(double x);  }
@Test public void beta() { Random r = RandomUtils.getRandom(); for (int i = 0; i < 200; i++) { double alpha = -50 * Math.log1p(-r.nextDouble()); double beta = -50 * Math.log1p(-r.nextDouble()); double ref = Math.exp(Gamma.logGamma(alpha) + Gamma.logGamma(beta) - Gamma.logGamma(alpha + beta)); double actual = Gamma.beta(alpha, beta); double err = (ref - actual) / ref; assertEquals("beta at (" + alpha + ", " + beta + ") relative error = " + err, 0, err, 1.0e-10); } }