Class BoostBeta

java.lang.Object
org.apache.commons.numbers.gamma.BoostBeta

final class BoostBeta extends 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:
  • Field Summary

    Fields
    Modifier and Type
    Field
    Description
    private static final double
    Default epsilon value for relative error.
    private static final double
    pi/2.
    private static final int
    Approximate value for ln(Double.MAX_VALUE).
    private static final int
    Approximate value for ln(Double.MIN_VALUE).
    private static final int
    The largest factorial that can be represented as a double.
    private static final int
    Size of the table of Pn's.
    private static final double
    2^53.
    private static final double
    2^-53.
  • Constructor Summary

    Constructors
    Modifier
    Constructor
    Description
    private
    Private constructor.
  • Method Summary

    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.

    Methods inherited from class java.lang.Object

    clone, equals, finalize, getClass, hashCode, notify, notifyAll, toString, wait, wait, wait
  • Field Details

    • EPSILON

      private static final double EPSILON
      Default epsilon value for relative error. This is equal to the Boost constant boost::math::EPSILON.
      See Also:
    • LOG_MAX_VALUE

      private static final int LOG_MAX_VALUE
      Approximate value for ln(Double.MAX_VALUE). This is equal to the Boost constant boost::math::tools::log_max_value<double>(). No term x should be used in exp(x) if x > LOG_MAX_VALUE to avoid overflow.
      See Also:
    • LOG_MIN_VALUE

      private static final int LOG_MIN_VALUE
      Approximate value for ln(Double.MIN_VALUE). This is equal to the Boost constant boost::math::tools::log_min_value<double>(). No term x should be used in exp(x) if x < LOG_MIN_VALUE to avoid underflow to sub-normal or zero.
      See Also:
    • HALF_PI

      private static final double HALF_PI
      pi/2.
      See Also:
    • MAX_FACTORIAL

      private static final int MAX_FACTORIAL
      The largest factorial that can be represented as a double. This is equal to the Boost constant boost::math::max_factorial<double>::value.
      See Also:
    • PN_SIZE

      private static final int PN_SIZE
      Size of the table of Pn's. Equal to Pn_size<double> suitable for 16-20 digit accuracy.
      See Also:
    • TWO_POW_53

      private static final double TWO_POW_53
      2^53. Used to scale sub-normal values.
      See Also:
    • TWO_POW_M53

      private static final double TWO_POW_M53
      2^-53. Used to rescale values.
      See Also:
  • Constructor Details

    • BoostBeta

      private BoostBeta()
      Private constructor.
  • Method Details

    • beta

      static double beta(double p, double q)
      Beta function.

      \[ B(p, q) = \frac{\Gamma(p) \Gamma(q)}{\Gamma(p+q)} \]

      Parameters:
      p - Argument p
      q - 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 a
      b - Argument b
      x - 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 of 1.0.

      Parameters:
      a - Argument a
      b - Argument b
      x - Argument x
      y - Argument 1-x
      normalised - 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 a
      b - Argument b
      x - Argument x
      y - Argument 1-x
      normalised - 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 a
      b - Argument b
      x - 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 a
      b - Argument b
      x - Argument x
      policy - 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 a
      b - Argument b
      x - 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 a
      b - Argument b
      x - Argument x
      policy - 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 a
      b - Argument b
      x - 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 a
      b - Argument b
      x - Argument x
      policy - 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 a
      b - Argument b
      x - 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 a
      b - Argument b
      x - Argument x
      policy - 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 from boost::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 a
      b - Argument b
      x - Argument x
      pol - Function evaluation policy
      normalised - true to compute the regularised value
      inv - 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 a
      b - Argument b
      x - Argument x
      s0 - Initial sum for the series
      normalised - true to compute the regularised value
      pol - 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 a
      b - Argument b
      k - 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 with k in [1, Integer.MAX_VALUE - 102])
      x - Argument x
      y - 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 n
      k - Argument k
      x - Argument x
      y - 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 a
      b - Argument b
      x - Argument x
      y - Argument 1-x
      k - Argument k
      normalised - 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 a
      b - Argument b
      x - Argument x
      y - Argument 1-x
      s0 - Initial sum for the series
      mult - Multiplication prefix factor
      pol - Function evaluation policy
      normalised - 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 a
      b - Argument b
      x - Argument x
      y - Argument 1-x
      pol - Function evaluation policy
      normalised - 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 call 1 - 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 a
      b - Argument b
      x - Argument x
      y - Argument 1-x
      pol - Function evaluation policy
      normalised - true to compute the regularised value
      Returns:
      incomplete beta