Beta.java

  1. /*
  2.  * Licensed to the Apache Software Foundation (ASF) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * The ASF licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *      https://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */

  18. /*
  19.  * This is not the original file distributed by the Apache Software Foundation
  20.  * It has been modified by the Hipparchus project
  21.  */
  22. package org.hipparchus.special;

  24. import org.hipparchus.exception.LocalizedCoreFormats;
  25. import org.hipparchus.exception.MathIllegalArgumentException;
  26. import org.hipparchus.util.ContinuedFraction;
  27. import org.hipparchus.util.FastMath;
  28. import org.hipparchus.util.MathUtils;

  30. /**
  31.  * <p>
  32.  * This is a utility class that provides computation methods related to the
  33.  * Beta family of functions.
  34.  * </p>
  35.  * <p>
  36.  * Implementation of {@link #logBeta(double, double)} is based on the
  37.  * algorithms described in
  38.  * </p>
  39.  * <ul>
  40.  * <li><a href="http://dx.doi.org/10.1145/22721.23109">Didonato and Morris
  41.  *     (1986)</a>, <em>Computation of the Incomplete Gamma Function Ratios
  42.  *     and their Inverse</em>, TOMS 12(4), 377-393,</li>
  43.  * <li><a href="http://dx.doi.org/10.1145/131766.131776">Didonato and Morris
  44.  *     (1992)</a>, <em>Algorithm 708: Significant Digit Computation of the
  45.  *     Incomplete Beta Function Ratios</em>, TOMS 18(3), 360-373,</li>
  46.  * </ul>
  47.  * <p>
  48.  * and implemented in the
  49.  * <a href="http://www.dtic.mil/docs/citations/ADA476840">NSWC Library of Mathematical Functions</a>,
  50.  * available
  51.  * <a href="http://www.ualberta.ca/CNS/RESEARCH/Software/NumericalNSWC/site.html">here</a>.
  52.  * This library is "approved for public release", and the
  53.  * <a href="http://www.dtic.mil/dtic/pdf/announcements/CopyrightGuidance.pdf">Copyright guidance</a>
  54.  * indicates that unless otherwise stated in the code, all FORTRAN functions in
  55.  * this library are license free. Since no such notice appears in the code these
  56.  * functions can safely be ported to Hipparchus.
  57.  * </p>
  58.  */
  59. public class Beta {
  60.     /** Maximum allowed numerical error. */
  61.     private static final double DEFAULT_EPSILON = 1E-14;

  63.     /** The constant value of ½log 2π. */
  64.     private static final double HALF_LOG_TWO_PI = .9189385332046727;

  66.     /**
  67.      * <p>
  68.      * The coefficients of the series expansion of the Δ function. This function
  69.      * is defined as follows
  70.      * </p>
  71.      * <center>Δ(x) = log Γ(x) - (x - 0.5) log a + a - 0.5 log 2π,</center>
  72.      * <p>
  73.      * see equation (23) in Didonato and Morris (1992). The series expansion,
  74.      * which applies for x ≥ 10, reads
  75.      * </p>
  76.      * <pre>
  77.      *                 14
  78.      *                ====
  79.      *             1  \                2 n
  80.      *     Δ(x) = ---  >    d  (10 / x)
  81.      *             x  /      n
  82.      *                ====
  83.      *                n = 0
  84.      * <pre>
  85.      */
  86.     private static final double[] DELTA = {
  87.         .833333333333333333333333333333E-01,
  88.         -.277777777777777777777777752282E-04,
  89.         .793650793650793650791732130419E-07,
  90.         -.595238095238095232389839236182E-09,
  91.         .841750841750832853294451671990E-11,
  92.         -.191752691751854612334149171243E-12,
  93.         .641025640510325475730918472625E-14,
  94.         -.295506514125338232839867823991E-15,
  95.         .179643716359402238723287696452E-16,
  96.         -.139228964661627791231203060395E-17,
  97.         .133802855014020915603275339093E-18,
  98.         -.154246009867966094273710216533E-19,
  99.         .197701992980957427278370133333E-20,
  100.         -.234065664793997056856992426667E-21,
  101.         .171348014966398575409015466667E-22
  102.     };

  104.     /**
  105.      * Default constructor.  Prohibit instantiation.
  106.      */
  107.     private Beta() {}

  109.     /**
  110.      * Returns the
  111.      * <a href="http://mathworld.wolfram.com/RegularizedBetaFunction.html">
  112.      * regularized beta function</a> I(x, a, b).
  113.      *
  114.      * @param x Value.
  115.      * @param a Parameter {@code a}.
  116.      * @param b Parameter {@code b}.
  117.      * @return the regularized beta function I(x, a, b).
  118.      * @throws org.hipparchus.exception.MathIllegalStateException
  119.      * if the algorithm fails to converge.
  120.      */
  121.     public static double regularizedBeta(double x, double a, double b) {
  122.         return regularizedBeta(x, a, b, DEFAULT_EPSILON, Integer.MAX_VALUE);
  123.     }

  125.     /**
  126.      * Returns the
  127.      * <a href="http://mathworld.wolfram.com/RegularizedBetaFunction.html">
  128.      * regularized beta function</a> I(x, a, b).
  129.      *
  130.      * @param x Value.
  131.      * @param a Parameter {@code a}.
  132.      * @param b Parameter {@code b}.
  133.      * @param epsilon When the absolute value of the nth item in the
  134.      * series is less than epsilon the approximation ceases to calculate
  135.      * further elements in the series.
  136.      * @return the regularized beta function I(x, a, b)
  137.      * @throws org.hipparchus.exception.MathIllegalStateException
  138.      * if the algorithm fails to converge.
  139.      */
  140.     public static double regularizedBeta(double x,
  141.                                          double a, double b,
  142.                                          double epsilon) {
  143.         return regularizedBeta(x, a, b, epsilon, Integer.MAX_VALUE);
  144.     }

  146.     /**
  147.      * Returns the regularized beta function I(x, a, b).
  148.      *
  149.      * @param x the value.
  150.      * @param a Parameter {@code a}.
  151.      * @param b Parameter {@code b}.
  152.      * @param maxIterations Maximum number of "iterations" to complete.
  153.      * @return the regularized beta function I(x, a, b)
  154.      * @throws org.hipparchus.exception.MathIllegalStateException
  155.      * if the algorithm fails to converge.
  156.      */
  157.     public static double regularizedBeta(double x,
  158.                                          double a, double b,
  159.                                          int maxIterations) {
  160.         return regularizedBeta(x, a, b, DEFAULT_EPSILON, maxIterations);
  161.     }

  163.     /**
  164.      * Returns the regularized beta function I(x, a, b).
  165.      *
  166.      * The implementation of this method is based on:
  167.      * <ul>
  168.      * <li>
  169.      * <a href="http://mathworld.wolfram.com/RegularizedBetaFunction.html">
  170.      * Regularized Beta Function</a>.</li>
  171.      * <li>
  172.      * <a href="http://functions.wolfram.com/06.21.10.0001.01">
  173.      * Regularized Beta Function</a>.</li>
  174.      * </ul>
  175.      *
  176.      * @param x the value.
  177.      * @param a Parameter {@code a}.
  178.      * @param b Parameter {@code b}.
  179.      * @param epsilon When the absolute value of the nth item in the
  180.      * series is less than epsilon the approximation ceases to calculate
  181.      * further elements in the series.
  182.      * @param maxIterations Maximum number of "iterations" to complete.
  183.      * @return the regularized beta function I(x, a, b)
  184.      * @throws org.hipparchus.exception.MathIllegalStateException
  185.      * if the algorithm fails to converge.
  186.      */
  187.     public static double regularizedBeta(double x,
  188.                                          final double a, final double b,
  189.                                          double epsilon, int maxIterations) {
  190.         double ret;

  192.         if (Double.isNaN(x) ||
  193.             Double.isNaN(a) ||
  194.             Double.isNaN(b) ||
  195.             x < 0 ||
  196.             x > 1 ||
  197.             a <= 0 ||
  198.             b <= 0) {
  199.             ret = Double.NaN;
  200.         } else if (x > (a + 1) / (2 + b + a) &&
  201.                    1 - x <= (b + 1) / (2 + b + a)) {
  202.             ret = 1 - regularizedBeta(1 - x, b, a, epsilon, maxIterations);
  203.         } else {
  204.             ContinuedFraction fraction = new ContinuedFraction() {

  206.                 /** {@inheritDoc} */
  207.                 @Override
  208.                 protected double getB(int n, double x) {
  209.                     double ret;
  210.                     double m;
  211.                     if (n % 2 == 0) { // even
  212.                         m = n / 2.0;
  213.                         ret = (m * (b - m) * x) /
  214.                             ((a + (2 * m) - 1) * (a + (2 * m)));
  215.                     } else {
  216.                         m = (n - 1.0) / 2.0;
  217.                         ret = -((a + m) * (a + b + m) * x) /
  218.                                 ((a + (2 * m)) * (a + (2 * m) + 1.0));
  219.                     }
  220.                     return ret;
  221.                 }

  223.                 /** {@inheritDoc} */
  224.                 @Override
  225.                 protected double getA(int n, double x) {
  226.                     return 1.0;
  227.                 }
  228.             };
  229.             ret = FastMath.exp((a * FastMath.log(x)) + (b * FastMath.log1p(-x)) -
  230.                 FastMath.log(a) - logBeta(a, b)) *
  231.                 1.0 / fraction.evaluate(x, epsilon, maxIterations);
  232.         }

  234.         return ret;
  235.     }

  237.     /**
  238.      * Returns the value of log Γ(a + b) for 1 ≤ a, b ≤ 2. Based on the
  239.      * <em>NSWC Library of Mathematics Subroutines</em> double precision
  240.      * implementation, {@code DGSMLN}. In {@code BetaTest.testLogGammaSum()},
  241.      * this private method is accessed through reflection.
  242.      *
  243.      * @param a First argument.
  244.      * @param b Second argument.
  245.      * @return the value of {@code log(Gamma(a + b))}.
  246.      * @throws MathIllegalArgumentException if {@code a} or {@code b} is lower than
  247.      * {@code 1.0} or greater than {@code 2.0}.
  248.      */
  249.     private static double logGammaSum(final double a, final double b)
  250.         throws MathIllegalArgumentException {

  252.         MathUtils.checkRangeInclusive(a, 1, 2);
  253.         MathUtils.checkRangeInclusive(b, 1, 2);

  255.         final double x = (a - 1.0) + (b - 1.0);
  256.         if (x <= 0.5) {
  257.             return Gamma.logGamma1p(1.0 + x);
  258.         } else if (x <= 1.5) {
  259.             return Gamma.logGamma1p(x) + FastMath.log1p(x);
  260.         } else {
  261.             return Gamma.logGamma1p(x - 1.0) + FastMath.log(x * (1.0 + x));
  262.         }
  263.     }

  265.     /**
  266.      * Returns the value of log[Γ(b) / Γ(a + b)] for a ≥ 0 and b ≥ 10. Based on
  267.      * the <em>NSWC Library of Mathematics Subroutines</em> double precision
  268.      * implementation, {@code DLGDIV}. In
  269.      * {@code BetaTest.testLogGammaMinusLogGammaSum()}, this private method is
  270.      * accessed through reflection.
  271.      *
  272.      * @param a First argument.
  273.      * @param b Second argument.
  274.      * @return the value of {@code log(Gamma(b) / Gamma(a + b))}.
  275.      * @throws MathIllegalArgumentException if {@code a < 0.0} or {@code b < 10.0}.
  276.      */
  277.     private static double logGammaMinusLogGammaSum(final double a,
  278.                                                    final double b)
  279.         throws MathIllegalArgumentException {

  281.         if (a < 0.0) {
  282.             throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_SMALL,
  283.                                                    a, 0.0);
  284.         }
  285.         if (b < 10.0) {
  286.             throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_SMALL,
  287.                                                    b, 10.0);
  288.         }

  290.         /*
  291.          * d = a + b - 0.5
  292.          */
  293.         final double d;
  294.         final double w;
  295.         if (a <= b) {
  296.             d = b + (a - 0.5);
  297.             w = deltaMinusDeltaSum(a, b);
  298.         } else {
  299.             d = a + (b - 0.5);
  300.             w = deltaMinusDeltaSum(b, a);
  301.         }

  303.         final double u = d * FastMath.log1p(a / b);
  304.         final double v = a * (FastMath.log(b) - 1.0);

  306.         return u <= v ? (w - u) - v : (w - v) - u;
  307.     }

  309.     /**
  310.      * Returns the value of Δ(b) - Δ(a + b), with 0 ≤ a ≤ b and b ≥ 10. Based
  311.      * on equations (26), (27) and (28) in Didonato and Morris (1992).
  312.      *
  313.      * @param a First argument.
  314.      * @param b Second argument.
  315.      * @return the value of {@code Delta(b) - Delta(a + b)}
  316.      * @throws MathIllegalArgumentException if {@code a < 0} or {@code a > b}
  317.      * @throws MathIllegalArgumentException if {@code b < 10}
  318.      */
  319.     private static double deltaMinusDeltaSum(final double a,
  320.                                              final double b)
  321.         throws MathIllegalArgumentException {

  323.         MathUtils.checkRangeInclusive(a, 0, b);
  324.         if (b < 10) {
  325.             throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_SMALL,
  326.                                                    b, 10);
  327.         }

  329.         final double h = a / b;
  330.         final double p = h / (1.0 + h);
  331.         final double q = 1.0 / (1.0 + h);
  332.         final double q2 = q * q;
  333.         /*
  334.          * s[i] = 1 + q + ... - q**(2 * i)
  335.          */
  336.         final double[] s = new double[DELTA.length];
  337.         s[0] = 1.0;
  338.         for (int i = 1; i < s.length; i++) {
  339.             s[i] = 1.0 + (q + q2 * s[i - 1]);
  340.         }
  341.         /*
  342.          * w = Delta(b) - Delta(a + b)
  343.          */
  344.         final double sqrtT = 10.0 / b;
  345.         final double t = sqrtT * sqrtT;
  346.         double w = DELTA[DELTA.length - 1] * s[s.length - 1];
  347.         for (int i = DELTA.length - 2; i >= 0; i--) {
  348.             w = t * w + DELTA[i] * s[i];
  349.         }
  350.         return w * p / b;
  351.     }

  353.     /**
  354.      * Returns the value of Δ(p) + Δ(q) - Δ(p + q), with p, q ≥ 10. Based on
  355.      * the <em>NSWC Library of Mathematics Subroutines</em> double precision
  356.      * implementation, {@code DBCORR}. In
  357.      * {@code BetaTest.testSumDeltaMinusDeltaSum()}, this private method is
  358.      * accessed through reflection.
  359.      *
  360.      * @param p First argument.
  361.      * @param q Second argument.
  362.      * @return the value of {@code Delta(p) + Delta(q) - Delta(p + q)}.
  363.      * @throws MathIllegalArgumentException if {@code p < 10.0} or {@code q < 10.0}.
  364.      */
  365.     private static double sumDeltaMinusDeltaSum(final double p,
  366.                                                 final double q) {

  368.         if (p < 10.0) {
  369.             throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_SMALL,
  370.                                                    p, 10.0);
  371.         }
  372.         if (q < 10.0) {
  373.             throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_SMALL,
  374.                                                    q, 10.0);
  375.         }

  377.         final double a = FastMath.min(p, q);
  378.         final double b = FastMath.max(p, q);
  379.         final double sqrtT = 10.0 / a;
  380.         final double t = sqrtT * sqrtT;
  381.         double z = DELTA[DELTA.length - 1];
  382.         for (int i = DELTA.length - 2; i >= 0; i--) {
  383.             z = t * z + DELTA[i];
  384.         }
  385.         return z / a + deltaMinusDeltaSum(a, b);
  386.     }

  388.     /**
  389.      * Returns the value of log B(p, q) for 0 ≤ x ≤ 1 and p, q &gt; 0. Based on the
  390.      * <em>NSWC Library of Mathematics Subroutines</em> implementation,
  391.      * {@code DBETLN}.
  392.      *
  393.      * @param p First argument.
  394.      * @param q Second argument.
  395.      * @return the value of {@code log(Beta(p, q))}, {@code NaN} if
  396.      * {@code p <= 0} or {@code q <= 0}.
  397.      */
  398.     public static double logBeta(final double p, final double q) {
  399.         if (Double.isNaN(p) || Double.isNaN(q) || (p <= 0.0) || (q <= 0.0)) {
  400.             return Double.NaN;
  401.         }

  403.         final double a = FastMath.min(p, q);
  404.         final double b = FastMath.max(p, q);
  405.         if (a >= 10.0) {
  406.             final double w = sumDeltaMinusDeltaSum(a, b);
  407.             final double h = a / b;
  408.             final double c = h / (1.0 + h);
  409.             final double u = -(a - 0.5) * FastMath.log(c);
  410.             final double v = b * FastMath.log1p(h);
  411.             if (u <= v) {
  412.                 return (((-0.5 * FastMath.log(b) + HALF_LOG_TWO_PI) + w) - u) - v;
  413.             } else {
  414.                 return (((-0.5 * FastMath.log(b) + HALF_LOG_TWO_PI) + w) - v) - u;
  415.             }
  416.         } else if (a > 2.0) {
  417.             if (b > 1000.0) {
  418.                 final int n = (int) FastMath.floor(a - 1.0);
  419.                 double prod = 1.0;
  420.                 double ared = a;
  421.                 for (int i = 0; i < n; i++) {
  422.                     ared -= 1.0;
  423.                     prod *= ared / (1.0 + ared / b);
  424.                 }
  425.                 return (FastMath.log(prod) - n * FastMath.log(b)) +
  426.                         (Gamma.logGamma(ared) +
  427.                          logGammaMinusLogGammaSum(ared, b));
  428.             } else {
  429.                 double prod1 = 1.0;
  430.                 double ared = a;
  431.                 while (ared > 2.0) {
  432.                     ared -= 1.0;
  433.                     final double h = ared / b;
  434.                     prod1 *= h / (1.0 + h);
  435.                 }
  436.                 if (b < 10.0) {
  437.                     double prod2 = 1.0;
  438.                     double bred = b;
  439.                     while (bred > 2.0) {
  440.                         bred -= 1.0;
  441.                         prod2 *= bred / (ared + bred);
  442.                     }
  443.                     return FastMath.log(prod1) +
  444.                            FastMath.log(prod2) +
  445.                            (Gamma.logGamma(ared) +
  446.                            (Gamma.logGamma(bred) -
  447.                             logGammaSum(ared, bred)));
  448.                 } else {
  449.                     return FastMath.log(prod1) +
  450.                            Gamma.logGamma(ared) +
  451.                            logGammaMinusLogGammaSum(ared, b);
  452.                 }
  453.             }
  454.         } else if (a >= 1.0) {
  455.             if (b > 2.0) {
  456.                 if (b < 10.0) {
  457.                     double prod = 1.0;
  458.                     double bred = b;
  459.                     while (bred > 2.0) {
  460.                         bred -= 1.0;
  461.                         prod *= bred / (a + bred);
  462.                     }
  463.                     return FastMath.log(prod) +
  464.                            (Gamma.logGamma(a) +
  465.                             (Gamma.logGamma(bred) -
  466.                              logGammaSum(a, bred)));
  467.                 } else {
  468.                     return Gamma.logGamma(a) +
  469.                            logGammaMinusLogGammaSum(a, b);
  470.                 }
  471.             } else {
  472.                 return Gamma.logGamma(a) +
  473.                        Gamma.logGamma(b) -
  474.                        logGammaSum(a, b);
  475.             }
  476.         } else {
  477.             if (b >= 10.0) {
  478.                 return Gamma.logGamma(a) +
  479.                        logGammaMinusLogGammaSum(a, b);
  480.             } else {
  481.                 // The following command is the original NSWC implementation.
  482.                 // return Gamma.logGamma(a) +
  483.                 // (Gamma.logGamma(b) - Gamma.logGamma(a + b));
  484.                 // The following command turns out to be more accurate.
  485.                 return FastMath.log(Gamma.gamma(a) * Gamma.gamma(b) /
  486.                                     Gamma.gamma(a + b));
  487.             }
  488.         }
  489.     }
  490. }
Created with JaCoCo 0.8.12.202403310830