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.util.ContinuedFraction;
020 import org.apache.commons.math.util.FastMath;
021
022 /**
023 * This is a utility class that provides computation methods related to the
024 * Beta family of functions.
025 *
026 * @version $Id: Beta.java 1131229 2011-06-03 20:49:25Z luc $
027 */
028 public class Beta {
029 /** Maximum allowed numerical error. */
030 private static final double DEFAULT_EPSILON = 10e-15;
031
032 /**
033 * Default constructor. Prohibit instantiation.
034 */
035 private Beta() {}
036
037 /**
038 * Returns the
039 * <a href="http://mathworld.wolfram.com/RegularizedBetaFunction.html">
040 * regularized beta function</a> I(x, a, b).
041 *
042 * @param x Value.
043 * @param a Parameter {@code a}.
044 * @param b Parameter {@code b}.
045 * @return the regularized beta function I(x, a, b).
046 * @throws org.apache.commons.math.exception.MaxCountExceededException
047 * if the algorithm fails to converge.
048 */
049 public static double regularizedBeta(double x, double a, double b) {
050 return regularizedBeta(x, a, b, DEFAULT_EPSILON, Integer.MAX_VALUE);
051 }
052
053 /**
054 * Returns the
055 * <a href="http://mathworld.wolfram.com/RegularizedBetaFunction.html">
056 * regularized beta function</a> I(x, a, b).
057 *
058 * @param x Value.
059 * @param a Parameter {@code a}.
060 * @param b Parameter {@code b}.
061 * @param epsilon When the absolute value of the nth item in the
062 * series is less than epsilon the approximation ceases to calculate
063 * further elements in the series.
064 * @return the regularized beta function I(x, a, b)
065 * @throws org.apache.commons.math.exception.MaxCountExceededException
066 * if the algorithm fails to converge.
067 */
068 public static double regularizedBeta(double x,
069 double a, double b,
070 double epsilon) {
071 return regularizedBeta(x, a, b, epsilon, Integer.MAX_VALUE);
072 }
073
074 /**
075 * Returns the regularized beta function I(x, a, b).
076 *
077 * @param x the value.
078 * @param a Parameter {@code a}.
079 * @param b Parameter {@code b}.
080 * @param maxIterations Maximum number of "iterations" to complete.
081 * @return the regularized beta function I(x, a, b)
082 * @throws org.apache.commons.math.exception.MaxCountExceededException
083 * if the algorithm fails to converge.
084 */
085 public static double regularizedBeta(double x,
086 double a, double b,
087 int maxIterations) {
088 return regularizedBeta(x, a, b, DEFAULT_EPSILON, maxIterations);
089 }
090
091 /**
092 * Returns the regularized beta function I(x, a, b).
093 *
094 * The implementation of this method is based on:
095 * <ul>
096 * <li>
097 * <a href="http://mathworld.wolfram.com/RegularizedBetaFunction.html">
098 * Regularized Beta Function</a>.</li>
099 * <li>
100 * <a href="http://functions.wolfram.com/06.21.10.0001.01">
101 * Regularized Beta Function</a>.</li>
102 * </ul>
103 *
104 * @param x the value.
105 * @param a Parameter {@code a}.
106 * @param b Parameter {@code b}.
107 * @param epsilon When the absolute value of the nth item in the
108 * series is less than epsilon the approximation ceases to calculate
109 * further elements in the series.
110 * @param maxIterations Maximum number of "iterations" to complete.
111 * @return the regularized beta function I(x, a, b)
112 * @throws org.apache.commons.math.exception.MaxCountExceededException
113 * if the algorithm fails to converge.
114 */
115 public static double regularizedBeta(double x,
116 final double a, final double b,
117 double epsilon, int maxIterations) {
118 double ret;
119
120 if (Double.isNaN(x) ||
121 Double.isNaN(a) ||
122 Double.isNaN(b) ||
123 x < 0 ||
124 x > 1 ||
125 a <= 0.0 ||
126 b <= 0.0) {
127 ret = Double.NaN;
128 } else if (x > (a + 1.0) / (a + b + 2.0)) {
129 ret = 1.0 - regularizedBeta(1.0 - x, b, a, epsilon, maxIterations);
130 } else {
131 ContinuedFraction fraction = new ContinuedFraction() {
132
133 @Override
134 protected double getB(int n, double x) {
135 double ret;
136 double m;
137 if (n % 2 == 0) { // even
138 m = n / 2.0;
139 ret = (m * (b - m) * x) /
140 ((a + (2 * m) - 1) * (a + (2 * m)));
141 } else {
142 m = (n - 1.0) / 2.0;
143 ret = -((a + m) * (a + b + m) * x) /
144 ((a + (2 * m)) * (a + (2 * m) + 1.0));
145 }
146 return ret;
147 }
148
149 @Override
150 protected double getA(int n, double x) {
151 return 1.0;
152 }
153 };
154 ret = FastMath.exp((a * FastMath.log(x)) + (b * FastMath.log(1.0 - x)) -
155 FastMath.log(a) - logBeta(a, b, epsilon, maxIterations)) *
156 1.0 / fraction.evaluate(x, epsilon, maxIterations);
157 }
158
159 return ret;
160 }
161
162 /**
163 * Returns the natural logarithm of the beta function B(a, b).
164 *
165 * @param a Parameter {@code a}.
166 * @param b Parameter {@code b}.
167 * @return log(B(a, b)).
168 */
169 public static double logBeta(double a, double b) {
170 return logBeta(a, b, DEFAULT_EPSILON, Integer.MAX_VALUE);
171 }
172
173 /**
174 * Returns the natural logarithm of the beta function B(a, b).
175 *
176 * The implementation of this method is based on:
177 * <ul>
178 * <li><a href="http://mathworld.wolfram.com/BetaFunction.html">
179 * Beta Function</a>, equation (1).</li>
180 * </ul>
181 *
182 * @param a Parameter {@code a}.
183 * @param b Parameter {@code b}.
184 * @param epsilon When the absolute value of the nth item in the
185 * series is less than epsilon the approximation ceases to calculate
186 * further elements in the series.
187 * @param maxIterations Maximum number of "iterations" to complete.
188 * @return log(B(a, b)).
189 */
190 public static double logBeta(double a, double b,
191 double epsilon,
192 int maxIterations) {
193 double ret;
194
195 if (Double.isNaN(a) ||
196 Double.isNaN(b) ||
197 a <= 0.0 ||
198 b <= 0.0) {
199 ret = Double.NaN;
200 } else {
201 ret = Gamma.logGamma(a) + Gamma.logGamma(b) -
202 Gamma.logGamma(a + b);
203 }
204
205 return ret;
206 }
207 }