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.math3.special;
018
019 import org.apache.commons.math3.exception.NumberIsTooSmallException;
020 import org.apache.commons.math3.exception.OutOfRangeException;
021 import org.apache.commons.math3.util.ContinuedFraction;
022 import org.apache.commons.math3.util.FastMath;
023
024 /**
025 * <p>
026 * This is a utility class that provides computation methods related to the
027 * Beta family of functions.
028 * </p>
029 * <p>
030 * Implementation of {@link #logBeta(double, double)} is based on the
031 * algorithms described in
032 * <ul>
033 * <li><a href="http://dx.doi.org/10.1145/22721.23109">Didonato and Morris
034 * (1986)</a>, <em>Computation of the Incomplete Gamma Function Ratios
035 * and their Inverse</em>, TOMS 12(4), 377-393,</li>
036 * <li><a href="http://dx.doi.org/10.1145/131766.131776">Didonato and Morris
037 * (1992)</a>, <em>Algorithm 708: Significant Digit Computation of the
038 * Incomplete Beta Function Ratios</em>, TOMS 18(3), 360-373,</li>
039 * </ul>
040 * and implemented in the
041 * <a href="http://www.dtic.mil/docs/citations/ADA476840">NSWC Library of Mathematical Functions</a>,
042 * available
043 * <a href="http://www.ualberta.ca/CNS/RESEARCH/Software/NumericalNSWC/site.html">here</a>.
044 * This library is "approved for public release", and the
045 * <a href="http://www.dtic.mil/dtic/pdf/announcements/CopyrightGuidance.pdf">Copyright guidance</a>
046 * indicates that unless otherwise stated in the code, all FORTRAN functions in
047 * this library are license free. Since no such notice appears in the code these
048 * functions can safely be ported to Commons-Math.
049 * </p>
050 *
051 *
052 * @version $Id: Beta.java 1420669 2012-12-12 13:40:35Z erans $
053 */
054 public class Beta {
055 /** Maximum allowed numerical error. */
056 private static final double DEFAULT_EPSILON = 1E-14;
057
058 /** The constant value of ½log 2π. */
059 private static final double HALF_LOG_TWO_PI = .9189385332046727;
060
061 /**
062 * <p>
063 * The coefficients of the series expansion of the Δ function. This function
064 * is defined as follows
065 * </p>
066 * <center>Δ(x) = log Γ(x) - (x - 0.5) log a + a - 0.5 log 2π,</center>
067 * <p>
068 * see equation (23) in Didonato and Morris (1992). The series expansion,
069 * which applies for x ≥ 10, reads
070 * </p>
071 * <pre>
072 * 14
073 * ====
074 * 1 \ 2 n
075 * Δ(x) = --- > d (10 / x)
076 * x / n
077 * ====
078 * n = 0
079 * <pre>
080 */
081 private static final double[] DELTA = {
082 .833333333333333333333333333333E-01,
083 -.277777777777777777777777752282E-04,
084 .793650793650793650791732130419E-07,
085 -.595238095238095232389839236182E-09,
086 .841750841750832853294451671990E-11,
087 -.191752691751854612334149171243E-12,
088 .641025640510325475730918472625E-14,
089 -.295506514125338232839867823991E-15,
090 .179643716359402238723287696452E-16,
091 -.139228964661627791231203060395E-17,
092 .133802855014020915603275339093E-18,
093 -.154246009867966094273710216533E-19,
094 .197701992980957427278370133333E-20,
095 -.234065664793997056856992426667E-21,
096 .171348014966398575409015466667E-22
097 };
098
099 /**
100 * Default constructor. Prohibit instantiation.
101 */
102 private Beta() {}
103
104 /**
105 * Returns the
106 * <a href="http://mathworld.wolfram.com/RegularizedBetaFunction.html">
107 * regularized beta function</a> I(x, a, b).
108 *
109 * @param x Value.
110 * @param a Parameter {@code a}.
111 * @param b Parameter {@code b}.
112 * @return the regularized beta function I(x, a, b).
113 * @throws org.apache.commons.math3.exception.MaxCountExceededException
114 * if the algorithm fails to converge.
115 */
116 public static double regularizedBeta(double x, double a, double b) {
117 return regularizedBeta(x, a, b, DEFAULT_EPSILON, Integer.MAX_VALUE);
118 }
119
120 /**
121 * Returns the
122 * <a href="http://mathworld.wolfram.com/RegularizedBetaFunction.html">
123 * regularized beta function</a> I(x, a, b).
124 *
125 * @param x Value.
126 * @param a Parameter {@code a}.
127 * @param b Parameter {@code b}.
128 * @param epsilon When the absolute value of the nth item in the
129 * series is less than epsilon the approximation ceases to calculate
130 * further elements in the series.
131 * @return the regularized beta function I(x, a, b)
132 * @throws org.apache.commons.math3.exception.MaxCountExceededException
133 * if the algorithm fails to converge.
134 */
135 public static double regularizedBeta(double x,
136 double a, double b,
137 double epsilon) {
138 return regularizedBeta(x, a, b, epsilon, Integer.MAX_VALUE);
139 }
140
141 /**
142 * Returns the regularized beta function I(x, a, b).
143 *
144 * @param x the value.
145 * @param a Parameter {@code a}.
146 * @param b Parameter {@code b}.
147 * @param maxIterations Maximum number of "iterations" to complete.
148 * @return the regularized beta function I(x, a, b)
149 * @throws org.apache.commons.math3.exception.MaxCountExceededException
150 * if the algorithm fails to converge.
151 */
152 public static double regularizedBeta(double x,
153 double a, double b,
154 int maxIterations) {
155 return regularizedBeta(x, a, b, DEFAULT_EPSILON, maxIterations);
156 }
157
158 /**
159 * Returns the regularized beta function I(x, a, b).
160 *
161 * The implementation of this method is based on:
162 * <ul>
163 * <li>
164 * <a href="http://mathworld.wolfram.com/RegularizedBetaFunction.html">
165 * Regularized Beta Function</a>.</li>
166 * <li>
167 * <a href="http://functions.wolfram.com/06.21.10.0001.01">
168 * Regularized Beta Function</a>.</li>
169 * </ul>
170 *
171 * @param x the value.
172 * @param a Parameter {@code a}.
173 * @param b Parameter {@code b}.
174 * @param epsilon When the absolute value of the nth item in the
175 * series is less than epsilon the approximation ceases to calculate
176 * further elements in the series.
177 * @param maxIterations Maximum number of "iterations" to complete.
178 * @return the regularized beta function I(x, a, b)
179 * @throws org.apache.commons.math3.exception.MaxCountExceededException
180 * if the algorithm fails to converge.
181 */
182 public static double regularizedBeta(double x,
183 final double a, final double b,
184 double epsilon, int maxIterations) {
185 double ret;
186
187 if (Double.isNaN(x) ||
188 Double.isNaN(a) ||
189 Double.isNaN(b) ||
190 x < 0 ||
191 x > 1 ||
192 a <= 0.0 ||
193 b <= 0.0) {
194 ret = Double.NaN;
195 } else if (x > (a + 1.0) / (a + b + 2.0)) {
196 ret = 1.0 - regularizedBeta(1.0 - x, b, a, epsilon, maxIterations);
197 } else {
198 ContinuedFraction fraction = new ContinuedFraction() {
199
200 @Override
201 protected double getB(int n, double x) {
202 double ret;
203 double m;
204 if (n % 2 == 0) { // even
205 m = n / 2.0;
206 ret = (m * (b - m) * x) /
207 ((a + (2 * m) - 1) * (a + (2 * m)));
208 } else {
209 m = (n - 1.0) / 2.0;
210 ret = -((a + m) * (a + b + m) * x) /
211 ((a + (2 * m)) * (a + (2 * m) + 1.0));
212 }
213 return ret;
214 }
215
216 @Override
217 protected double getA(int n, double x) {
218 return 1.0;
219 }
220 };
221 ret = FastMath.exp((a * FastMath.log(x)) + (b * FastMath.log(1.0 - x)) -
222 FastMath.log(a) - logBeta(a, b)) *
223 1.0 / fraction.evaluate(x, epsilon, maxIterations);
224 }
225
226 return ret;
227 }
228
229 /**
230 * Returns the natural logarithm of the beta function B(a, b).
231 *
232 * The implementation of this method is based on:
233 * <ul>
234 * <li><a href="http://mathworld.wolfram.com/BetaFunction.html">
235 * Beta Function</a>, equation (1).</li>
236 * </ul>
237 *
238 * @param a Parameter {@code a}.
239 * @param b Parameter {@code b}.
240 * @param epsilon This parameter is ignored.
241 * @param maxIterations This parameter is ignored.
242 * @return log(B(a, b)).
243 * @deprecated as of version 3.1, this method is deprecated as the
244 * computation of the beta function is no longer iterative; it will be
245 * removed in version 4.0. Current implementation of this method
246 * internally calls {@link #logBeta(double, double)}.
247 */
248 @Deprecated
249 public static double logBeta(double a, double b,
250 double epsilon,
251 int maxIterations) {
252
253 return logBeta(a, b);
254 }
255
256
257 /**
258 * Returns the value of log Γ(a + b) for 1 ≤ a, b ≤ 2. Based on the
259 * <em>NSWC Library of Mathematics Subroutines</em> double precision
260 * implementation, {@code DGSMLN}. In {@link BetaTest#testLogGammaSum()},
261 * this private method is accessed through reflection.
262 *
263 * @param a First argument.
264 * @param b Second argument.
265 * @return the value of {@code log(Gamma(a + b))}.
266 * @throws OutOfRangeException if {@code a} or {@code b} is lower than
267 * {@code 1.0} or greater than {@code 2.0}.
268 */
269 private static double logGammaSum(final double a, final double b)
270 throws OutOfRangeException {
271
272 if ((a < 1.0) || (a > 2.0)) {
273 throw new OutOfRangeException(a, 1.0, 2.0);
274 }
275 if ((b < 1.0) || (b > 2.0)) {
276 throw new OutOfRangeException(b, 1.0, 2.0);
277 }
278
279 final double x = (a - 1.0) + (b - 1.0);
280 if (x <= 0.5) {
281 return Gamma.logGamma1p(1.0 + x);
282 } else if (x <= 1.5) {
283 return Gamma.logGamma1p(x) + FastMath.log1p(x);
284 } else {
285 return Gamma.logGamma1p(x - 1.0) + FastMath.log(x * (1.0 + x));
286 }
287 }
288
289 /**
290 * Returns the value of log[Γ(b) / Γ(a + b)] for a ≥ 0 and b ≥ 10. Based on
291 * the <em>NSWC Library of Mathematics Subroutines</em> double precision
292 * implementation, {@code DLGDIV}. In
293 * {@link BetaTest#testLogGammaMinusLogGammaSum()}, this private method is
294 * accessed through reflection.
295 *
296 * @param a First argument.
297 * @param b Second argument.
298 * @return the value of {@code log(Gamma(b) / Gamma(a + b))}.
299 * @throws NumberIsTooSmallException if {@code a < 0.0} or {@code b < 10.0}.
300 */
301 private static double logGammaMinusLogGammaSum(final double a,
302 final double b)
303 throws NumberIsTooSmallException {
304
305 if (a < 0.0) {
306 throw new NumberIsTooSmallException(a, 0.0, true);
307 }
308 if (b < 10.0) {
309 throw new NumberIsTooSmallException(b, 10.0, true);
310 }
311
312 /*
313 * d = a + b - 0.5
314 */
315 final double d;
316 final double w;
317 if (a <= b) {
318 d = b + (a - 0.5);
319 w = deltaMinusDeltaSum(a, b);
320 } else {
321 d = a + (b - 0.5);
322 w = deltaMinusDeltaSum(b, a);
323 }
324
325 final double u = d * FastMath.log1p(a / b);
326 final double v = a * (FastMath.log(b) - 1.0);
327
328 return u <= v ? (w - u) - v : (w - v) - u;
329 }
330
331 /**
332 * Returns the value of Δ(b) - Δ(a + b), with 0 ≤ a ≤ b and b ≥ 10. Based
333 * on equations (26), (27) and (28) in Didonato and Morris (1992).
334 *
335 * @param a First argument.
336 * @param b Second argument.
337 * @return the value of {@code Delta(b) - Delta(a + b)}
338 * @throws OutOfRangeException if {@code a < 0} or {@code a > b}
339 * @throws NumberIsTooSmallException if {@code b < 10}
340 */
341 private static double deltaMinusDeltaSum(final double a,
342 final double b)
343 throws OutOfRangeException, NumberIsTooSmallException {
344
345 if ((a < 0) || (a > b)) {
346 throw new OutOfRangeException(a, 0, b);
347 }
348 if (b < 10) {
349 throw new NumberIsTooSmallException(b, 10, true);
350 }
351
352 final double h = a / b;
353 final double p = h / (1.0 + h);
354 final double q = 1.0 / (1.0 + h);
355 final double q2 = q * q;
356 /*
357 * s[i] = 1 + q + ... - q**(2 * i)
358 */
359 final double[] s = new double[DELTA.length];
360 s[0] = 1.0;
361 for (int i = 1; i < s.length; i++) {
362 s[i] = 1.0 + (q + q2 * s[i - 1]);
363 }
364 /*
365 * w = Delta(b) - Delta(a + b)
366 */
367 final double sqrtT = 10.0 / b;
368 final double t = sqrtT * sqrtT;
369 double w = DELTA[DELTA.length - 1] * s[s.length - 1];
370 for (int i = DELTA.length - 2; i >= 0; i--) {
371 w = t * w + DELTA[i] * s[i];
372 }
373 return w * p / b;
374 }
375
376 /**
377 * Returns the value of Δ(p) + Δ(q) - Δ(p + q), with p, q ≥ 10. Based on
378 * the <em>NSWC Library of Mathematics Subroutines</em> double precision
379 * implementation, {@code DBCORR}. In
380 * {@link BetaTest#testSumDeltaMinusDeltaSum()}, this private method is
381 * accessed through reflection.
382 *
383 * @param p First argument.
384 * @param q Second argument.
385 * @return the value of {@code Delta(p) + Delta(q) - Delta(p + q)}.
386 * @throws NumberIsTooSmallException if {@code p < 10.0} or {@code q < 10.0}.
387 */
388 private static double sumDeltaMinusDeltaSum(final double p,
389 final double q) {
390
391 if (p < 10.0) {
392 throw new NumberIsTooSmallException(p, 10.0, true);
393 }
394 if (q < 10.0) {
395 throw new NumberIsTooSmallException(q, 10.0, true);
396 }
397
398 final double a = FastMath.min(p, q);
399 final double b = FastMath.max(p, q);
400 final double sqrtT = 10.0 / a;
401 final double t = sqrtT * sqrtT;
402 double z = DELTA[DELTA.length - 1];
403 for (int i = DELTA.length - 2; i >= 0; i--) {
404 z = t * z + DELTA[i];
405 }
406 return z / a + deltaMinusDeltaSum(a, b);
407 }
408
409 /**
410 * Returns the value of log B(p, q) for 0 ≤ x ≤ 1 and p, q > 0. Based on the
411 * <em>NSWC Library of Mathematics Subroutines</em> implementation,
412 * {@code DBETLN}.
413 *
414 * @param p First argument.
415 * @param q Second argument.
416 * @return the value of {@code log(Beta(p, q))}, {@code NaN} if
417 * {@code p <= 0} or {@code q <= 0}.
418 */
419 public static double logBeta(final double p, final double q) {
420 if (Double.isNaN(p) || Double.isNaN(q) || (p <= 0.0) || (q <= 0.0)) {
421 return Double.NaN;
422 }
423
424 final double a = FastMath.min(p, q);
425 final double b = FastMath.max(p, q);
426 if (a >= 10.0) {
427 final double w = sumDeltaMinusDeltaSum(a, b);
428 final double h = a / b;
429 final double c = h / (1.0 + h);
430 final double u = -(a - 0.5) * FastMath.log(c);
431 final double v = b * FastMath.log1p(h);
432 if (u <= v) {
433 return (((-0.5 * FastMath.log(b) + HALF_LOG_TWO_PI) + w) - u) - v;
434 } else {
435 return (((-0.5 * FastMath.log(b) + HALF_LOG_TWO_PI) + w) - v) - u;
436 }
437 } else if (a > 2.0) {
438 if (b > 1000.0) {
439 final int n = (int) FastMath.floor(a - 1.0);
440 double prod = 1.0;
441 double ared = a;
442 for (int i = 0; i < n; i++) {
443 ared -= 1.0;
444 prod *= ared / (1.0 + ared / b);
445 }
446 return (FastMath.log(prod) - n * FastMath.log(b)) +
447 (Gamma.logGamma(ared) +
448 logGammaMinusLogGammaSum(ared, b));
449 } else {
450 double prod1 = 1.0;
451 double ared = a;
452 while (ared > 2.0) {
453 ared -= 1.0;
454 final double h = ared / b;
455 prod1 *= h / (1.0 + h);
456 }
457 if (b < 10.0) {
458 double prod2 = 1.0;
459 double bred = b;
460 while (bred > 2.0) {
461 bred -= 1.0;
462 prod2 *= bred / (ared + bred);
463 }
464 return FastMath.log(prod1) +
465 FastMath.log(prod2) +
466 (Gamma.logGamma(ared) +
467 (Gamma.logGamma(bred) -
468 logGammaSum(ared, bred)));
469 } else {
470 return FastMath.log(prod1) +
471 Gamma.logGamma(ared) +
472 logGammaMinusLogGammaSum(ared, b);
473 }
474 }
475 } else if (a >= 1.0) {
476 if (b > 2.0) {
477 if (b < 10.0) {
478 double prod = 1.0;
479 double bred = b;
480 while (bred > 2.0) {
481 bred -= 1.0;
482 prod *= bred / (a + bred);
483 }
484 return FastMath.log(prod) +
485 (Gamma.logGamma(a) +
486 (Gamma.logGamma(bred) -
487 logGammaSum(a, bred)));
488 } else {
489 return Gamma.logGamma(a) +
490 logGammaMinusLogGammaSum(a, b);
491 }
492 } else {
493 return Gamma.logGamma(a) +
494 Gamma.logGamma(b) -
495 logGammaSum(a, b);
496 }
497 } else {
498 if (b >= 10.0) {
499 return Gamma.logGamma(a) +
500 logGammaMinusLogGammaSum(a, b);
501 } else {
502 // The following command is the original NSWC implementation.
503 // return Gamma.logGamma(a) +
504 // (Gamma.logGamma(b) - Gamma.logGamma(a + b));
505 // The following command turns out to be more accurate.
506 return FastMath.log(Gamma.gamma(a) * Gamma.gamma(b) /
507 Gamma.gamma(a + b));
508 }
509 }
510 }
511 }