Class BoostBeta

    • 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.
      • Methods inherited from class java.lang.Object

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

      • EPSILON

        private static final double EPSILON
        Default epsilon value for relative error. This is equal to the Boost constant boost::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 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:
        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 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:
        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 constant boost::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 to Pn_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
    • Constructor Detail

      • BoostBeta

        private BoostBeta()
        Private constructor.
    • 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 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