Coverage Report - org.apache.commons.math.special.Beta

Classes in this File Line Coverage Branch Coverage Complexity
Beta
84% 
100% 
1.444

 1  
 /*
 2  
  * Copyright 2003-2004 The Apache Software Foundation.
 3  
  *
 4  
  * Licensed under the Apache License, Version 2.0 (the "License");
 5  
  * you may not use this file except in compliance with the License.
 6  
  * You may obtain a copy of the License at
 7  
  *
 8  
  *      http://www.apache.org/licenses/LICENSE-2.0
 9  
  *
 10  
  * Unless required by applicable law or agreed to in writing, software
 11  
  * distributed under the License is distributed on an "AS IS" BASIS,
 12  
  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 13  
  * See the License for the specific language governing permissions and
 14  
  * limitations under the License.
 15  
  */
 16  
 package org.apache.commons.math.special;
 17  
 
 18  
 import java.io.Serializable;
 19  
 
 20  
 import org.apache.commons.math.MathException;
 21  
 import org.apache.commons.math.util.ContinuedFraction;
 22  
 
 23  
 /**
 24  
  * This is a utility class that provides computation methods related to the
 25  
  * Beta family of functions.
 26  
  *
 27  
  * @version $Revision$ $Date: 2005-02-26 05:11:52 -0800 (Sat, 26 Feb 2005) $
 28  
  */
 29  
 public class Beta implements Serializable {
 30  
     /** Maximum allowed numerical error. */
 31  
     private static final double DEFAULT_EPSILON = 10e-9;
 32  
 
 33  
     /**
 34  
      * Default constructor.  Prohibit instantiation.
 35  
      */
 36  
     private Beta() {
 37  0
         super();
 38  0
     }
 39  
 
 40  
     /**
 41  
      * Returns the
 42  
      * <a href="http://mathworld.wolfram.com/RegularizedBetaFunction.html">
 43  
      * regularized beta function</a> I(x, a, b).
 44  
      * 
 45  
      * @param x the value.
 46  
      * @param a the a parameter.
 47  
      * @param b the b parameter.
 48  
      * @return the regularized beta function I(x, a, b)
 49  
      * @throws MathException if the algorithm fails to converge.
 50  
      */
 51  
     public static double regularizedBeta(double x, double a, double b)
 52  
         throws MathException
 53  
     {
 54  4698
         return regularizedBeta(x, a, b, DEFAULT_EPSILON, Integer.MAX_VALUE);
 55  
     }
 56  
 
 57  
     /**
 58  
      * Returns the
 59  
      * <a href="http://mathworld.wolfram.com/RegularizedBetaFunction.html">
 60  
      * regularized beta function</a> I(x, a, b).
 61  
      * 
 62  
      * @param x the value.
 63  
      * @param a the a parameter.
 64  
      * @param b the b parameter.
 65  
      * @param epsilon When the absolute value of the nth item in the
 66  
      *                series is less than epsilon the approximation ceases
 67  
      *                to calculate further elements in the series.
 68  
      * @return the regularized beta function I(x, a, b)
 69  
      * @throws MathException if the algorithm fails to converge.
 70  
      */
 71  
     public static double regularizedBeta(double x, double a, double b,
 72  
         double epsilon) throws MathException
 73  
     {
 74  0
         return regularizedBeta(x, a, b, epsilon, Integer.MAX_VALUE);
 75  
     }
 76  
 
 77  
     /**
 78  
      * Returns the regularized beta function I(x, a, b).
 79  
      * 
 80  
      * @param x the value.
 81  
      * @param a the a parameter.
 82  
      * @param b the b parameter.
 83  
      * @param maxIterations Maximum number of "iterations" to complete. 
 84  
      * @return the regularized beta function I(x, a, b)
 85  
      * @throws MathException if the algorithm fails to converge.
 86  
      */
 87  
     public static double regularizedBeta(double x, double a, double b,
 88  
         int maxIterations) throws MathException
 89  
     {
 90  0
         return regularizedBeta(x, a, b, DEFAULT_EPSILON, maxIterations);
 91  
     }
 92  
     
 93  
     /**
 94  
      * Returns the regularized beta function I(x, a, b).
 95  
      * 
 96  
      * The implementation of this method is based on:
 97  
      * <ul>
 98  
      * <li>
 99  
      * <a href="http://mathworld.wolfram.com/RegularizedBetaFunction.html">
 100  
      * Regularized Beta Function</a>.</li>
 101  
      * <li>
 102  
      * <a href="http://functions.wolfram.com/06.21.10.0001.01">
 103  
      * Regularized Beta Function</a>.</li>
 104  
      * </ul>
 105  
      * 
 106  
      * @param x the value.
 107  
      * @param a the a parameter.
 108  
      * @param b the b parameter.
 109  
      * @param epsilon When the absolute value of the nth item in the
 110  
      *                series is less than epsilon the approximation ceases
 111  
      *                to calculate further elements in the series.
 112  
      * @param maxIterations Maximum number of "iterations" to complete. 
 113  
      * @return the regularized beta function I(x, a, b)
 114  
      * @throws MathException if the algorithm fails to converge.
 115  
      */
 116  
     public static double regularizedBeta(double x, final double a,
 117  
         final double b, double epsilon, int maxIterations) throws MathException
 118  
     {
 119  
         double ret;
 120  
 
 121  5282
         if (Double.isNaN(x) || Double.isNaN(a) || Double.isNaN(b) || (x < 0) ||
 122  
             (x > 1) || (a <= 0.0) || (b <= 0.0))
 123  
         {
 124  16
             ret = Double.NaN;
 125  5266
         } else if (x > (a + 1.0) / (a + b + 2.0)) {
 126  584
             ret = 1.0 - regularizedBeta(1.0 - x, b, a, epsilon, maxIterations);
 127  
         } else {
 128  4682
             ContinuedFraction fraction = new ContinuedFraction() {
 129  
                 protected double getB(int n, double x) {
 130  
                     double ret;
 131  
                     double m;
 132  22142
                     if (n % 2 == 0) { // even
 133  10498
                         m = n / 2.0;
 134  10498
                         ret = (m * (b - m) * x) /
 135  
                             ((a + (2 * m) - 1) * (a + (2 * m)));
 136  
                     } else {
 137  11644
                         m = (n - 1.0) / 2.0;
 138  11644
                         ret = -((a + m) * (a + b + m) * x) /
 139  
                                 ((a + (2 * m)) * (a + (2 * m) + 1.0));
 140  
                     }
 141  22142
                     return ret;
 142  
                 }
 143  
 
 144  4682
                 protected double getA(int n, double x) {
 145  26824
                     return 1.0;
 146  
                 }
 147  
             };
 148  4682
             ret = Math.exp((a * Math.log(x)) + (b * Math.log(1.0 - x)) -
 149  
                 Math.log(a) - logBeta(a, b, epsilon, maxIterations)) *
 150  
                 1.0 / fraction.evaluate(x, epsilon, maxIterations);
 151  
         }
 152  
 
 153  5282
         return ret;
 154  
     }
 155  
 
 156  
     /**
 157  
      * Returns the natural logarithm of the beta function B(a, b).
 158  
      * 
 159  
      * @param a the a parameter.
 160  
      * @param b the b parameter.
 161  
      * @return log(B(a, b))
 162  
      */
 163  
     public static double logBeta(double a, double b) {
 164  14
         return logBeta(a, b, DEFAULT_EPSILON, Integer.MAX_VALUE);
 165  
     }
 166  
     
 167  
     /**
 168  
      * Returns the natural logarithm of the beta function B(a, b).
 169  
      *
 170  
      * The implementation of this method is based on:
 171  
      * <ul>
 172  
      * <li><a href="http://mathworld.wolfram.com/BetaFunction.html">
 173  
      * Beta Function</a>, equation (1).</li>
 174  
      * </ul>
 175  
      * 
 176  
      * @param a the a parameter.
 177  
      * @param b the b parameter.
 178  
      * @param epsilon When the absolute value of the nth item in the
 179  
      *                series is less than epsilon the approximation ceases
 180  
      *                to calculate further elements in the series.
 181  
      * @param maxIterations Maximum number of "iterations" to complete. 
 182  
      * @return log(B(a, b))
 183  
      */
 184  
     public static double logBeta(double a, double b, double epsilon,
 185  
         int maxIterations) {
 186  
             
 187  
         double ret;
 188  
 
 189  4696
         if (Double.isNaN(a) || Double.isNaN(b) || (a <= 0.0) || (b <= 0.0)) {
 190  12
             ret = Double.NaN;
 191  
         } else {
 192  4684
             ret = Gamma.logGamma(a) + Gamma.logGamma(b) -
 193  
                 Gamma.logGamma(a + b);
 194  
         }
 195  
 
 196  4696
         return ret;
 197  
     }
 198  
 }