Class BoostBeta
- java.lang.Object
-
- org.apache.commons.numbers.gamma.BoostBeta
-
final class BoostBeta extends java.lang.Object
Implementation of the regularized beta functions and incomplete beta functions.This code has been adapted from the Boost
c++
implementation<boost/math/special_functions/beta.hpp>
. All work is copyright to the original authors and subject to the Boost Software License.- See Also:
- Boost C++ beta functions
-
-
Field Summary
Fields Modifier and Type Field Description private static double
EPSILON
Default epsilon value for relative error.private static double
HALF_PI
pi/2.private static int
LOG_MAX_VALUE
Approximate value for ln(Double.MAX_VALUE).private static int
LOG_MIN_VALUE
Approximate value for ln(Double.MIN_VALUE).private static int
MAX_FACTORIAL
The largest factorial that can be represented as a double.private static int
PN_SIZE
Size of the table of Pn's.private static double
TWO_POW_53
2^53.private static double
TWO_POW_M53
2^-53.
-
Constructor Summary
Constructors Modifier Constructor Description private
BoostBeta()
Private constructor.
-
Method Summary
All Methods Static Methods Concrete Methods Modifier and Type Method Description (package private) static double
beta(double p, double q)
Beta function.(package private) static double
beta(double a, double b, double x)
Full incomplete beta.(package private) static double
beta(double a, double b, double x, Policy policy)
Full incomplete beta.(package private) static double
betac(double a, double b, double x)
Complement of the incomplete beta.(package private) static double
betac(double a, double b, double x, Policy policy)
Complement of the incomplete beta.private static double
betaIncompleteImp(double a, double b, double x, Policy pol, boolean normalised, boolean inv)
Main incomplete beta entry point, handles all four incomplete betas.private static double
betaSmallBLargeASeries(double a, double b, double x, double y, double s0, double mult, Policy pol, boolean normalised)
Computes beta(a, b, x) using small b and large a.private static double
binomialCCdf(int n, int k, double x, double y)
Binomial complementary cdf.(package private) static double
binomialCoefficient(int n, int k)
Computes the binomial coefficient.private static double
binomialTerm(int n, int k, double x, double y)
Compute the binomial term.(package private) static double
ibeta(double a, double b, double x)
Regularised incomplete beta.(package private) static double
ibeta(double a, double b, double x, Policy policy)
Regularised incomplete beta.private static double
ibetaAStep(double a, double b, double x, double y, int k, boolean normalised)
Computes the difference between ibeta(a,b,x) and ibeta(a+k,b,x).(package private) static double
ibetac(double a, double b, double x)
Complement of the regularised incomplete beta.(package private) static double
ibetac(double a, double b, double x, Policy policy)
Complement of the regularised incomplete beta.(package private) static double
ibetaDerivative(double a, double b, double x)
Derivative of the regularised incomplete beta.(package private) static double
ibetaFraction(double a, double b, double x, double y, Policy pol, boolean normalised)
Evaluate the incomplete beta via the classic continued fraction representation.(package private) static double
ibetaFraction2(double a, double b, double x, double y, Policy pol, boolean normalised)
Evaluate the incomplete beta via a continued fraction representation.private static double
ibetaPowerTerms(double a, double b, double x, double y, boolean normalised)
Compute the leading power terms in the incomplete Beta.private static double
ibetaPowerTerms(double a, double b, double x, double y, boolean normalised, double prefix)
Compute the leading power terms in the incomplete Beta.private static double
ibetaSeries(double a, double b, double x, double s0, boolean normalised, Policy pol)
Series approximation to the incomplete beta.private static double
risingFactorialRatio(double a, double b, int k)
Rising factorial ratio.
-
-
-
Field Detail
-
EPSILON
private static final double EPSILON
Default epsilon value for relative error. This is equal to the Boost constantboost::math::EPSILON
.- See Also:
- Constant Field Values
-
LOG_MAX_VALUE
private static final int LOG_MAX_VALUE
Approximate value for ln(Double.MAX_VALUE). This is equal to the Boost constantboost::math::tools::log_max_value<double>()
. No termx
should be used inexp(x)
ifx > LOG_MAX_VALUE
to avoid overflow.- See Also:
- Constant Field Values
-
LOG_MIN_VALUE
private static final int LOG_MIN_VALUE
Approximate value for ln(Double.MIN_VALUE). This is equal to the Boost constantboost::math::tools::log_min_value<double>()
. No termx
should be used inexp(x)
ifx < LOG_MIN_VALUE
to avoid underflow to sub-normal or zero.- See Also:
- Constant Field Values
-
HALF_PI
private static final double HALF_PI
pi/2.- See Also:
- Constant Field Values
-
MAX_FACTORIAL
private static final int MAX_FACTORIAL
The largest factorial that can be represented as a double. This is equal to the Boost constantboost::math::max_factorial<double>::value
.- See Also:
- Constant Field Values
-
PN_SIZE
private static final int PN_SIZE
Size of the table of Pn's. Equal toPn_size<double>
suitable for 16-20 digit accuracy.- See Also:
- Constant Field Values
-
TWO_POW_53
private static final double TWO_POW_53
2^53. Used to scale sub-normal values.- See Also:
- Constant Field Values
-
TWO_POW_M53
private static final double TWO_POW_M53
2^-53. Used to rescale values.- See Also:
- Constant Field Values
-
-
Method Detail
-
beta
static double beta(double p, double q)
Beta function.\[ B(p, q) = \frac{\Gamma(p) \Gamma(q)}{\Gamma(p+q)} \]
- Parameters:
p
- Argument pq
- Argument q- Returns:
- beta value
-
ibetaDerivative
static double ibetaDerivative(double a, double b, double x)
Derivative of the regularised incomplete beta.\[ \frac{\delta}{\delta x} I_x(a, b) = \frac{(1-x)^{b-1} x^{a-1}}{\B(a, b)} \]
Adapted from
boost::math::ibeta_derivative
.- Parameters:
a
- Argument ab
- Argument bx
- Argument x- Returns:
- ibeta derivative
-
ibetaPowerTerms
private static double ibetaPowerTerms(double a, double b, double x, double y, boolean normalised)
Compute the leading power terms in the incomplete Beta.Utility function to call
ibetaPowerTerms(double, double, double, double, boolean, double)
using a multiplication prefix of1.0
.- Parameters:
a
- Argument ab
- Argument bx
- Argument xy
- Argument 1-xnormalised
- true to divide by beta(a, b)- Returns:
- incomplete beta power terms
-
ibetaPowerTerms
private static double ibetaPowerTerms(double a, double b, double x, double y, boolean normalised, double prefix)
Compute the leading power terms in the incomplete Beta.(x^a)(y^b)/Beta(a,b) when normalised, and (x^a)(y^b) otherwise.
Almost all of the error in the incomplete beta comes from this function: particularly when a and b are large. Computing large powers are *hard* though, and using logarithms just leads to horrendous cancellation errors.
- Parameters:
a
- Argument ab
- Argument bx
- Argument xy
- Argument 1-xnormalised
- true to divide by beta(a, b)prefix
- Prefix to multiply by the result- Returns:
- incomplete beta power terms
-
beta
static double beta(double a, double b, double x)
Full incomplete beta.\[ B_x(a,b) = \int_0^x t^{a-1}\,(1-t)^{b-1}\,dt \]
- Parameters:
a
- Argument ab
- Argument bx
- Argument x- Returns:
- lower beta value
-
beta
static double beta(double a, double b, double x, Policy policy)
Full incomplete beta.\[ B_x(a,b) = \int_0^x t^{a-1}\,(1-t)^{b-1}\,dt \]
- Parameters:
a
- Argument ab
- Argument bx
- Argument xpolicy
- Function evaluation policy- Returns:
- lower beta value
-
betac
static double betac(double a, double b, double x)
Complement of the incomplete beta.\[ betac(a, b, x) = beta(b, a, 1-x) = B_{1-x}(b,a) \]
- Parameters:
a
- Argument ab
- Argument bx
- Argument x- Returns:
- upper beta value
-
betac
static double betac(double a, double b, double x, Policy policy)
Complement of the incomplete beta.\[ betac(a, b, x) = beta(b, a, 1-x) = B_{1-x}(b,a) \]
- Parameters:
a
- Argument ab
- Argument bx
- Argument xpolicy
- Function evaluation policy- Returns:
- upper beta value
-
ibeta
static double ibeta(double a, double b, double x)
Regularised incomplete beta.\[ I_x(a,b) = \frac{1}{B(a, b)} \int_0^x t^{a-1}\,(1-t)^{b-1}\,dt \]
- Parameters:
a
- Argument ab
- Argument bx
- Argument x- Returns:
- p
-
ibeta
static double ibeta(double a, double b, double x, Policy policy)
Regularised incomplete beta.\[ I_x(a,b) = \frac{1}{B(a, b)} \int_0^x t^{a-1}\,(1-t)^{b-1}\,dt \]
- Parameters:
a
- Argument ab
- Argument bx
- Argument xpolicy
- Function evaluation policy- Returns:
- p
-
ibetac
static double ibetac(double a, double b, double x)
Complement of the regularised incomplete beta.\[ ibetac(a, b, x) = 1 - I_x(a,b) = I_{1-x}(b, a) \]
- Parameters:
a
- Argument ab
- Argument bx
- Argument x- Returns:
- q
-
ibetac
static double ibetac(double a, double b, double x, Policy policy)
Complement of the regularised incomplete beta.\[ ibetac(a, b, x) = 1 - I_x(a,b) = I_{1-x}(b, a) \]
- Parameters:
a
- Argument ab
- Argument bx
- Argument xpolicy
- Function evaluation policy- Returns:
- q
-
betaIncompleteImp
private static double betaIncompleteImp(double a, double b, double x, Policy pol, boolean normalised, boolean inv)
Main incomplete beta entry point, handles all four incomplete betas. Adapted fromboost::math::detail::ibeta_imp
.The Boost code has a pointer
p_derivative
that can be set to the value of the derivative. This is used for the inverse incomplete beta functions. It is not required for the forward evaluation functions.- Parameters:
a
- Argument ab
- Argument bx
- Argument xpol
- Function evaluation policynormalised
- true to compute the regularised valueinv
- true to compute the complement value- Returns:
- incomplete beta value
-
ibetaSeries
private static double ibetaSeries(double a, double b, double x, double s0, boolean normalised, Policy pol)
Series approximation to the incomplete beta.- Parameters:
a
- Argument ab
- Argument bx
- Argument xs0
- Initial sum for the seriesnormalised
- true to compute the regularised valuepol
- Function evaluation policy- Returns:
- incomplete beta series
-
risingFactorialRatio
private static double risingFactorialRatio(double a, double b, int k)
Rising factorial ratio. This function is only needed for the non-regular incomplete beta, it computes the delta in:beta(a,b,x) = prefix + delta * beta(a+k,b,x)
It is currently only called for small k.
- Parameters:
a
- Argument ab
- Argument bk
- Argument k- Returns:
- the rising factorial ratio
-
binomialCCdf
private static double binomialCCdf(int n, int k, double x, double y)
Binomial complementary cdf. For integer arguments we can relate the incomplete beta to the complement of the binomial distribution cdf and use this finite sum.- Parameters:
n
- Argument n (called with[2, 39] + k
)k
- Argument k (called withk in [1, Integer.MAX_VALUE - 102]
)x
- Argument xy
- Argument 1-x- Returns:
- Binomial complementary cdf
-
binomialTerm
private static double binomialTerm(int n, int k, double x, double y)
Compute the binomial term.x^k * (1-x)^(n-k) * binomial(n, k)
This is a helper function used to guard against extreme values generated in the term which can produce NaN from zero multiplied by infinity.
- Parameters:
n
- Argument nk
- Argument kx
- Argument xy
- Argument 1-x- Returns:
- Binomial term
-
binomialCoefficient
static double binomialCoefficient(int n, int k)
Computes the binomial coefficient.Adapted from
org.apache.commons.numbers.combinatorics.BinomialCoefficientDouble
.Note: This does not use
BinomialCoefficientDouble
to avoid a circular dependency as the combinatorics depends on the gamma package. No checks are made on the arguments.- Parameters:
n
- Size of the set (must be positive).k
- Size of the subsets to be counted (must be in [0, n]).- Returns:
n choose k
.
-
ibetaAStep
private static double ibetaAStep(double a, double b, double x, double y, int k, boolean normalised)
Computes the difference between ibeta(a,b,x) and ibeta(a+k,b,x).- Parameters:
a
- Argument ab
- Argument bx
- Argument xy
- Argument 1-xk
- Argument knormalised
- true to compute the regularised value- Returns:
- ibeta difference
-
betaSmallBLargeASeries
private static double betaSmallBLargeASeries(double a, double b, double x, double y, double s0, double mult, Policy pol, boolean normalised)
Computes beta(a, b, x) using small b and large a.- Parameters:
a
- Argument ab
- Argument bx
- Argument xy
- Argument 1-xs0
- Initial sum for the seriesmult
- Multiplication prefix factorpol
- Function evaluation policynormalised
- true to compute the regularised value- Returns:
- beta series
-
ibetaFraction2
static double ibetaFraction2(double a, double b, double x, double y, Policy pol, boolean normalised)
Evaluate the incomplete beta via a continued fraction representation.Note: This is not a generic continued fraction. The formula is from Didonato and Morris and is only used when a and b are above 1. See Incomplete Beta Function: Implementation.
- Parameters:
a
- Argument ab
- Argument bx
- Argument xy
- Argument 1-xpol
- Function evaluation policynormalised
- true to compute the regularised value- Returns:
- incomplete beta
-
ibetaFraction
static double ibetaFraction(double a, double b, double x, double y, Policy pol, boolean normalised)
Evaluate the incomplete beta via the classic continued fraction representation.Note: This is not part of the Boost C++ implementation. This is a generic continued fraction applicable to all arguments. It is used when alternative methods are unsuitable due to the presence of sub-normal terms and the result is expected to be sub-normal. In this case the original Boost implementation would return zero.
This continued fraction was the evaluation method used in Commons Numbers 1.0. This method has been updated to use the
ibetaPowerTerms
function to compute the power terms. Reversal of the arguments to call1 - ibetaFraction(b, a, 1 - x)
is not performed as the result is expected to be very small and this optimisation for accuracy is not required.- Parameters:
a
- Argument ab
- Argument bx
- Argument xy
- Argument 1-xpol
- Function evaluation policynormalised
- true to compute the regularised value- Returns:
- incomplete beta
-
-