in commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/BoostBeta.java [523:822]
private static double betaIncompleteImp(double a, double b, double x,
Policy pol, boolean normalised, boolean inv) {
//
// The incomplete beta function implementation:
// This is just a big bunch of spaghetti code to divide up the
// input range and select the right implementation method for
// each domain:
//
if (!(x >= 0 && x <= 1)) {
// Domain error
return Double.NaN;
}
if (normalised) {
if (!(a >= 0 && b >= 0)) {
// Domain error
return Double.NaN;
}
// extend to a few very special cases:
if (a == 0) {
if (b == 0) {
// a and b cannot both be zero
return Double.NaN;
}
// Assume b > 0
return inv ? 0 : 1;
} else if (b == 0) {
// assume a > 0
return inv ? 1 : 0;
}
} else {
if (!(a > 0 && b > 0)) {
// Domain error
return Double.NaN;
}
}
if (x == 0) {
if (inv) {
return normalised ? 1 : beta(a, b);
}
return 0;
}
if (x == 1) {
if (!inv) {
return normalised ? 1 : beta(a, b);
}
return 0;
}
if (a == 0.5f && b == 0.5f) {
// We have an arcsine distribution:
final double z = inv ? 1 - x : x;
final double asin = Math.asin(Math.sqrt(z));
return normalised ? asin / HALF_PI : 2 * asin;
}
boolean invert = inv;
double y = 1 - x;
if (a == 1) {
// swap(a, b)
double tmp = a;
a = b;
b = tmp;
// swap(x, y)
tmp = x;
x = y;
y = tmp;
invert = !invert;
}
if (b == 1) {
//
// Special case see:
// http://functions.wolfram.com/GammaBetaErf/BetaRegularized/03/01/01/
//
if (a == 1) {
return invert ? y : x;
}
double p;
if (y < 0.5) {
p = invert ? -Math.expm1(a * Math.log1p(-y)) : Math.exp(a * Math.log1p(-y));
} else {
p = invert ? -BoostMath.powm1(x, a) : Math.pow(x, a);
}
if (!normalised) {
p /= a;
}
return p;
}
double fract;
if (Math.min(a, b) <= 1) {
if (x > 0.5) {
// swap(a, b)
double tmp = a;
a = b;
b = tmp;
// swap(x, y)
tmp = x;
x = y;
y = tmp;
invert = !invert;
}
if (Math.max(a, b) <= 1) {
// Both a,b < 1:
if (a >= Math.min(0.2, b) || Math.pow(x, a) <= 0.9) {
if (invert) {
fract = -(normalised ? 1 : beta(a, b));
invert = false;
fract = -ibetaSeries(a, b, x, fract, normalised, pol);
} else {
fract = ibetaSeries(a, b, x, 0, normalised, pol);
}
} else {
// swap(a, b)
double tmp = a;
a = b;
b = tmp;
// swap(x, y)
tmp = x;
x = y;
y = tmp;
invert = !invert;
if (y >= 0.3) {
if (invert) {
fract = -(normalised ? 1 : beta(a, b));
invert = false;
fract = -ibetaSeries(a, b, x, fract, normalised, pol);
} else {
fract = ibetaSeries(a, b, x, 0, normalised, pol);
}
} else {
// Sidestep on a, and then use the series representation:
final double prefix;
if (normalised) {
prefix = 1;
} else {
prefix = risingFactorialRatio(a + b, a, 20);
}
fract = ibetaAStep(a, b, x, y, 20, normalised);
if (invert) {
fract -= normalised ? 1 : beta(a, b);
invert = false;
fract = -betaSmallBLargeASeries(a + 20, b, x, y, fract, prefix, pol, normalised);
} else {
fract = betaSmallBLargeASeries(a + 20, b, x, y, fract, prefix, pol, normalised);
}
}
}
} else {
// One of a, b < 1 only:
if (b <= 1 || (x < 0.1 && Math.pow(b * x, a) <= 0.7)) {
if (invert) {
fract = -(normalised ? 1 : beta(a, b));
invert = false;
fract = -ibetaSeries(a, b, x, fract, normalised, pol);
} else {
fract = ibetaSeries(a, b, x, 0, normalised, pol);
}
} else {
// swap(a, b)
double tmp = a;
a = b;
b = tmp;
// swap(x, y)
tmp = x;
x = y;
y = tmp;
invert = !invert;
if (y >= 0.3) {
if (invert) {
fract = -(normalised ? 1 : beta(a, b));
invert = false;
fract = -ibetaSeries(a, b, x, fract, normalised, pol);
} else {
fract = ibetaSeries(a, b, x, 0, normalised, pol);
}
} else if (a >= 15) {
if (invert) {
fract = -(normalised ? 1 : beta(a, b));
invert = false;
fract = -betaSmallBLargeASeries(a, b, x, y, fract, 1, pol, normalised);
} else {
fract = betaSmallBLargeASeries(a, b, x, y, 0, 1, pol, normalised);
}
} else {
// Sidestep to improve errors:
final double prefix;
if (normalised) {
prefix = 1;
} else {
prefix = risingFactorialRatio(a + b, a, 20);
}
fract = ibetaAStep(a, b, x, y, 20, normalised);
if (invert) {
fract -= normalised ? 1 : beta(a, b);
invert = false;
fract = -betaSmallBLargeASeries(a + 20, b, x, y, fract, prefix, pol, normalised);
} else {
fract = betaSmallBLargeASeries(a + 20, b, x, y, fract, prefix, pol, normalised);
}
}
}
}
} else {
// Both a,b >= 1:
// Note:
// median ~ (a - 1/3) / (a + b - 2/3) ~ a / (a + b)
// if x > a / (a + b) => a - (a + b) * x < 0
final double lambda;
if (a < b) {
lambda = a - (a + b) * x;
} else {
lambda = (a + b) * y - b;
}
if (lambda < 0) {
// swap(a, b)
double tmp = a;
a = b;
b = tmp;
// swap(x, y)
tmp = x;
x = y;
y = tmp;
invert = !invert;
}
if (b < 40) {
// Note: y != 1 check is required for non-zero x < epsilon
if (Math.rint(a) == a && Math.rint(b) == b && a < (Integer.MAX_VALUE - 100) && y != 1) {
// Here: a in [2, 2^31 - 102] && b in [2, 39]
// relate to the binomial distribution and use a finite sum:
final int k = (int) (a - 1);
final int n = (int) (b + k);
fract = binomialCCdf(n, k, x, y);
if (!normalised) {
fract *= beta(a, b);
}
} else if (b * x <= 0.7) {
if (invert) {
fract = -(normalised ? 1 : beta(a, b));
invert = false;
fract = -ibetaSeries(a, b, x, fract, normalised, pol);
} else {
fract = ibetaSeries(a, b, x, 0, normalised, pol);
}
} else if (a > 15) {
// sidestep so we can use the series representation:
int n = (int) b;
if (n == b) {
--n;
}
final double bbar = b - n;
final double prefix;
if (normalised) {
prefix = 1;
} else {
prefix = risingFactorialRatio(a + bbar, bbar, n);
}
fract = ibetaAStep(bbar, a, y, x, n, normalised);
fract = betaSmallBLargeASeries(a, bbar, x, y, fract, 1, pol, normalised);
fract /= prefix;
} else if (normalised) {
// The formula here for the non-normalised case is tricky to figure
// out (for me!!), and requires two pochhammer calculations rather
// than one, so leave it for now and only use this in the normalized case....
int n = (int) Math.floor(b);
double bbar = b - n;
if (bbar <= 0) {
--n;
bbar += 1;
}
fract = ibetaAStep(bbar, a, y, x, n, normalised);
fract += ibetaAStep(a, bbar, x, y, 20, normalised);
if (invert) {
// Note this line would need changing if we ever enable this branch in
// non-normalized case
fract -= 1;
}
fract = betaSmallBLargeASeries(a + 20, bbar, x, y, fract, 1, pol, normalised);
if (invert) {
fract = -fract;
invert = false;
}
} else {
fract = ibetaFraction2(a, b, x, y, pol, normalised);
}
} else {
fract = ibetaFraction2(a, b, x, y, pol, normalised);
}
}
if (invert) {
return (normalised ? 1 : beta(a, b)) - fract;
}
return fract;
}