Classes in this File | Line Coverage | Branch Coverage | Complexity | ||||||||
Beta |
|
| 1.4444444444444444;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 | } |