001 /* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the NOTICE file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017 package org.apache.commons.math.special; 018 019 import org.apache.commons.math.MathException; 020 import org.apache.commons.math.util.ContinuedFraction; 021 import org.apache.commons.math.util.FastMath; 022 023 /** 024 * This is a utility class that provides computation methods related to the 025 * Beta family of functions. 026 * 027 * @version $Revision: 990655 $ $Date: 2010-08-29 23:49:40 +0200 (dim. 29 ao??t 2010) $ 028 */ 029 public class Beta { 030 031 /** Maximum allowed numerical error. */ 032 private static final double DEFAULT_EPSILON = 10e-15; 033 034 /** 035 * Default constructor. Prohibit instantiation. 036 */ 037 private Beta() { 038 super(); 039 } 040 041 /** 042 * Returns the 043 * <a href="http://mathworld.wolfram.com/RegularizedBetaFunction.html"> 044 * regularized beta function</a> I(x, a, b). 045 * 046 * @param x the value. 047 * @param a the a parameter. 048 * @param b the b parameter. 049 * @return the regularized beta function I(x, a, b) 050 * @throws MathException if the algorithm fails to converge. 051 */ 052 public static double regularizedBeta(double x, double a, double b) 053 throws MathException 054 { 055 return regularizedBeta(x, a, b, DEFAULT_EPSILON, Integer.MAX_VALUE); 056 } 057 058 /** 059 * Returns the 060 * <a href="http://mathworld.wolfram.com/RegularizedBetaFunction.html"> 061 * regularized beta function</a> I(x, a, b). 062 * 063 * @param x the value. 064 * @param a the a parameter. 065 * @param b the b parameter. 066 * @param epsilon When the absolute value of the nth item in the 067 * series is less than epsilon the approximation ceases 068 * to calculate further elements in the series. 069 * @return the regularized beta function I(x, a, b) 070 * @throws MathException if the algorithm fails to converge. 071 */ 072 public static double regularizedBeta(double x, double a, double b, 073 double epsilon) throws MathException 074 { 075 return regularizedBeta(x, a, b, epsilon, Integer.MAX_VALUE); 076 } 077 078 /** 079 * Returns the regularized beta function I(x, a, b). 080 * 081 * @param x the value. 082 * @param a the a parameter. 083 * @param b the b parameter. 084 * @param maxIterations Maximum number of "iterations" to complete. 085 * @return the regularized beta function I(x, a, b) 086 * @throws MathException if the algorithm fails to converge. 087 */ 088 public static double regularizedBeta(double x, double a, double b, 089 int maxIterations) throws MathException 090 { 091 return regularizedBeta(x, a, b, DEFAULT_EPSILON, maxIterations); 092 } 093 094 /** 095 * Returns the regularized beta function I(x, a, b). 096 * 097 * The implementation of this method is based on: 098 * <ul> 099 * <li> 100 * <a href="http://mathworld.wolfram.com/RegularizedBetaFunction.html"> 101 * Regularized Beta Function</a>.</li> 102 * <li> 103 * <a href="http://functions.wolfram.com/06.21.10.0001.01"> 104 * Regularized Beta Function</a>.</li> 105 * </ul> 106 * 107 * @param x the value. 108 * @param a the a parameter. 109 * @param b the b parameter. 110 * @param epsilon When the absolute value of the nth item in the 111 * series is less than epsilon the approximation ceases 112 * to calculate further elements in the series. 113 * @param maxIterations Maximum number of "iterations" to complete. 114 * @return the regularized beta function I(x, a, b) 115 * @throws MathException if the algorithm fails to converge. 116 */ 117 public static double regularizedBeta(double x, final double a, 118 final double b, double epsilon, int maxIterations) throws MathException 119 { 120 double ret; 121 122 if (Double.isNaN(x) || Double.isNaN(a) || Double.isNaN(b) || (x < 0) || 123 (x > 1) || (a <= 0.0) || (b <= 0.0)) 124 { 125 ret = Double.NaN; 126 } else if (x > (a + 1.0) / (a + b + 2.0)) { 127 ret = 1.0 - regularizedBeta(1.0 - x, b, a, epsilon, maxIterations); 128 } else { 129 ContinuedFraction fraction = new ContinuedFraction() { 130 131 @Override 132 protected double getB(int n, double x) { 133 double ret; 134 double m; 135 if (n % 2 == 0) { // even 136 m = n / 2.0; 137 ret = (m * (b - m) * x) / 138 ((a + (2 * m) - 1) * (a + (2 * m))); 139 } else { 140 m = (n - 1.0) / 2.0; 141 ret = -((a + m) * (a + b + m) * x) / 142 ((a + (2 * m)) * (a + (2 * m) + 1.0)); 143 } 144 return ret; 145 } 146 147 @Override 148 protected double getA(int n, double x) { 149 return 1.0; 150 } 151 }; 152 ret = FastMath.exp((a * FastMath.log(x)) + (b * FastMath.log(1.0 - x)) - 153 FastMath.log(a) - logBeta(a, b, epsilon, maxIterations)) * 154 1.0 / fraction.evaluate(x, epsilon, maxIterations); 155 } 156 157 return ret; 158 } 159 160 /** 161 * Returns the natural logarithm of the beta function B(a, b). 162 * 163 * @param a the a parameter. 164 * @param b the b parameter. 165 * @return log(B(a, b)) 166 */ 167 public static double logBeta(double a, double b) { 168 return logBeta(a, b, DEFAULT_EPSILON, Integer.MAX_VALUE); 169 } 170 171 /** 172 * Returns the natural logarithm of the beta function B(a, b). 173 * 174 * The implementation of this method is based on: 175 * <ul> 176 * <li><a href="http://mathworld.wolfram.com/BetaFunction.html"> 177 * Beta Function</a>, equation (1).</li> 178 * </ul> 179 * 180 * @param a the a parameter. 181 * @param b the b parameter. 182 * @param epsilon When the absolute value of the nth item in the 183 * series is less than epsilon the approximation ceases 184 * to calculate further elements in the series. 185 * @param maxIterations Maximum number of "iterations" to complete. 186 * @return log(B(a, b)) 187 */ 188 public static double logBeta(double a, double b, double epsilon, 189 int maxIterations) { 190 191 double ret; 192 193 if (Double.isNaN(a) || Double.isNaN(b) || (a <= 0.0) || (b <= 0.0)) { 194 ret = Double.NaN; 195 } else { 196 ret = Gamma.logGamma(a) + Gamma.logGamma(b) - 197 Gamma.logGamma(a + b); 198 } 199 200 return ret; 201 } 202 }