1 /*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * https://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17
18 /*
19 * This is not the original file distributed by the Apache Software Foundation
20 * It has been modified by the Hipparchus project
21 */
22 package org.hipparchus.special;
23
24 import org.hipparchus.exception.LocalizedCoreFormats;
25 import org.hipparchus.exception.MathIllegalArgumentException;
26 import org.hipparchus.util.ContinuedFraction;
27 import org.hipparchus.util.FastMath;
28 import org.hipparchus.util.MathUtils;
29
30 /**
31 * <p>
32 * This is a utility class that provides computation methods related to the
33 * Beta family of functions.
34 * </p>
35 * <p>
36 * Implementation of {@link #logBeta(double, double)} is based on the
37 * algorithms described in
38 * </p>
39 * <ul>
40 * <li><a href="http://dx.doi.org/10.1145/22721.23109">Didonato and Morris
41 * (1986)</a>, <em>Computation of the Incomplete Gamma Function Ratios
42 * and their Inverse</em>, TOMS 12(4), 377-393,</li>
43 * <li><a href="http://dx.doi.org/10.1145/131766.131776">Didonato and Morris
44 * (1992)</a>, <em>Algorithm 708: Significant Digit Computation of the
45 * Incomplete Beta Function Ratios</em>, TOMS 18(3), 360-373,</li>
46 * </ul>
47 * <p>
48 * and implemented in the
49 * <a href="http://www.dtic.mil/docs/citations/ADA476840">NSWC Library of Mathematical Functions</a>,
50 * available
51 * <a href="http://www.ualberta.ca/CNS/RESEARCH/Software/NumericalNSWC/site.html">here</a>.
52 * This library is "approved for public release", and the
53 * <a href="http://www.dtic.mil/dtic/pdf/announcements/CopyrightGuidance.pdf">Copyright guidance</a>
54 * indicates that unless otherwise stated in the code, all FORTRAN functions in
55 * this library are license free. Since no such notice appears in the code these
56 * functions can safely be ported to Hipparchus.
57 * </p>
58 */
59 public class Beta {
60 /** Maximum allowed numerical error. */
61 private static final double DEFAULT_EPSILON = 1E-14;
62
63 /** The constant value of ½log 2π. */
64 private static final double HALF_LOG_TWO_PI = .9189385332046727;
65
66 /**
67 * <p>
68 * The coefficients of the series expansion of the Δ function. This function
69 * is defined as follows
70 * </p>
71 * <center>Δ(x) = log Γ(x) - (x - 0.5) log a + a - 0.5 log 2π,</center>
72 * <p>
73 * see equation (23) in Didonato and Morris (1992). The series expansion,
74 * which applies for x ≥ 10, reads
75 * </p>
76 * <pre>
77 * 14
78 * ====
79 * 1 \ 2 n
80 * Δ(x) = --- > d (10 / x)
81 * x / n
82 * ====
83 * n = 0
84 * <pre>
85 */
86 private static final double[] DELTA = {
87 .833333333333333333333333333333E-01,
88 -.277777777777777777777777752282E-04,
89 .793650793650793650791732130419E-07,
90 -.595238095238095232389839236182E-09,
91 .841750841750832853294451671990E-11,
92 -.191752691751854612334149171243E-12,
93 .641025640510325475730918472625E-14,
94 -.295506514125338232839867823991E-15,
95 .179643716359402238723287696452E-16,
96 -.139228964661627791231203060395E-17,
97 .133802855014020915603275339093E-18,
98 -.154246009867966094273710216533E-19,
99 .197701992980957427278370133333E-20,
100 -.234065664793997056856992426667E-21,
101 .171348014966398575409015466667E-22
102 };
103
104 /**
105 * Default constructor. Prohibit instantiation.
106 */
107 private Beta() {}
108
109 /**
110 * Returns the
111 * <a href="http://mathworld.wolfram.com/RegularizedBetaFunction.html">
112 * regularized beta function</a> I(x, a, b).
113 *
114 * @param x Value.
115 * @param a Parameter {@code a}.
116 * @param b Parameter {@code b}.
117 * @return the regularized beta function I(x, a, b).
118 * @throws org.hipparchus.exception.MathIllegalStateException
119 * if the algorithm fails to converge.
120 */
121 public static double regularizedBeta(double x, double a, double b) {
122 return regularizedBeta(x, a, b, DEFAULT_EPSILON, Integer.MAX_VALUE);
123 }
124
125 /**
126 * Returns the
127 * <a href="http://mathworld.wolfram.com/RegularizedBetaFunction.html">
128 * regularized beta function</a> I(x, a, b).
129 *
130 * @param x Value.
131 * @param a Parameter {@code a}.
132 * @param b Parameter {@code b}.
133 * @param epsilon When the absolute value of the nth item in the
134 * series is less than epsilon the approximation ceases to calculate
135 * further elements in the series.
136 * @return the regularized beta function I(x, a, b)
137 * @throws org.hipparchus.exception.MathIllegalStateException
138 * if the algorithm fails to converge.
139 */
140 public static double regularizedBeta(double x,
141 double a, double b,
142 double epsilon) {
143 return regularizedBeta(x, a, b, epsilon, Integer.MAX_VALUE);
144 }
145
146 /**
147 * Returns the regularized beta function I(x, a, b).
148 *
149 * @param x the value.
150 * @param a Parameter {@code a}.
151 * @param b Parameter {@code b}.
152 * @param maxIterations Maximum number of "iterations" to complete.
153 * @return the regularized beta function I(x, a, b)
154 * @throws org.hipparchus.exception.MathIllegalStateException
155 * if the algorithm fails to converge.
156 */
157 public static double regularizedBeta(double x,
158 double a, double b,
159 int maxIterations) {
160 return regularizedBeta(x, a, b, DEFAULT_EPSILON, maxIterations);
161 }
162
163 /**
164 * Returns the regularized beta function I(x, a, b).
165 *
166 * The implementation of this method is based on:
167 * <ul>
168 * <li>
169 * <a href="http://mathworld.wolfram.com/RegularizedBetaFunction.html">
170 * Regularized Beta Function</a>.</li>
171 * <li>
172 * <a href="http://functions.wolfram.com/06.21.10.0001.01">
173 * Regularized Beta Function</a>.</li>
174 * </ul>
175 *
176 * @param x the value.
177 * @param a Parameter {@code a}.
178 * @param b Parameter {@code b}.
179 * @param epsilon When the absolute value of the nth item in the
180 * series is less than epsilon the approximation ceases to calculate
181 * further elements in the series.
182 * @param maxIterations Maximum number of "iterations" to complete.
183 * @return the regularized beta function I(x, a, b)
184 * @throws org.hipparchus.exception.MathIllegalStateException
185 * if the algorithm fails to converge.
186 */
187 public static double regularizedBeta(double x,
188 final double a, final double b,
189 double epsilon, int maxIterations) {
190 double ret;
191
192 if (Double.isNaN(x) ||
193 Double.isNaN(a) ||
194 Double.isNaN(b) ||
195 x < 0 ||
196 x > 1 ||
197 a <= 0 ||
198 b <= 0) {
199 ret = Double.NaN;
200 } else if (x > (a + 1) / (2 + b + a) &&
201 1 - x <= (b + 1) / (2 + b + a)) {
202 ret = 1 - regularizedBeta(1 - x, b, a, epsilon, maxIterations);
203 } else {
204 ContinuedFraction fraction = new ContinuedFraction() {
205
206 /** {@inheritDoc} */
207 @Override
208 protected double getB(int n, double x) {
209 double ret;
210 double m;
211 if (n % 2 == 0) { // even
212 m = n / 2.0;
213 ret = (m * (b - m) * x) /
214 ((a + (2 * m) - 1) * (a + (2 * m)));
215 } else {
216 m = (n - 1.0) / 2.0;
217 ret = -((a + m) * (a + b + m) * x) /
218 ((a + (2 * m)) * (a + (2 * m) + 1.0));
219 }
220 return ret;
221 }
222
223 /** {@inheritDoc} */
224 @Override
225 protected double getA(int n, double x) {
226 return 1.0;
227 }
228 };
229 ret = FastMath.exp((a * FastMath.log(x)) + (b * FastMath.log1p(-x)) -
230 FastMath.log(a) - logBeta(a, b)) *
231 1.0 / fraction.evaluate(x, epsilon, maxIterations);
232 }
233
234 return ret;
235 }
236
237 /**
238 * Returns the value of log Γ(a + b) for 1 ≤ a, b ≤ 2. Based on the
239 * <em>NSWC Library of Mathematics Subroutines</em> double precision
240 * implementation, {@code DGSMLN}. In {@code BetaTest.testLogGammaSum()},
241 * this private method is accessed through reflection.
242 *
243 * @param a First argument.
244 * @param b Second argument.
245 * @return the value of {@code log(Gamma(a + b))}.
246 * @throws MathIllegalArgumentException if {@code a} or {@code b} is lower than
247 * {@code 1.0} or greater than {@code 2.0}.
248 */
249 private static double logGammaSum(final double a, final double b)
250 throws MathIllegalArgumentException {
251
252 MathUtils.checkRangeInclusive(a, 1, 2);
253 MathUtils.checkRangeInclusive(b, 1, 2);
254
255 final double x = (a - 1.0) + (b - 1.0);
256 if (x <= 0.5) {
257 return Gamma.logGamma1p(1.0 + x);
258 } else if (x <= 1.5) {
259 return Gamma.logGamma1p(x) + FastMath.log1p(x);
260 } else {
261 return Gamma.logGamma1p(x - 1.0) + FastMath.log(x * (1.0 + x));
262 }
263 }
264
265 /**
266 * Returns the value of log[Γ(b) / Γ(a + b)] for a ≥ 0 and b ≥ 10. Based on
267 * the <em>NSWC Library of Mathematics Subroutines</em> double precision
268 * implementation, {@code DLGDIV}. In
269 * {@code BetaTest.testLogGammaMinusLogGammaSum()}, this private method is
270 * accessed through reflection.
271 *
272 * @param a First argument.
273 * @param b Second argument.
274 * @return the value of {@code log(Gamma(b) / Gamma(a + b))}.
275 * @throws MathIllegalArgumentException if {@code a < 0.0} or {@code b < 10.0}.
276 */
277 private static double logGammaMinusLogGammaSum(final double a,
278 final double b)
279 throws MathIllegalArgumentException {
280
281 if (a < 0.0) {
282 throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_SMALL,
283 a, 0.0);
284 }
285 if (b < 10.0) {
286 throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_SMALL,
287 b, 10.0);
288 }
289
290 /*
291 * d = a + b - 0.5
292 */
293 final double d;
294 final double w;
295 if (a <= b) {
296 d = b + (a - 0.5);
297 w = deltaMinusDeltaSum(a, b);
298 } else {
299 d = a + (b - 0.5);
300 w = deltaMinusDeltaSum(b, a);
301 }
302
303 final double u = d * FastMath.log1p(a / b);
304 final double v = a * (FastMath.log(b) - 1.0);
305
306 return u <= v ? (w - u) - v : (w - v) - u;
307 }
308
309 /**
310 * Returns the value of Δ(b) - Δ(a + b), with 0 ≤ a ≤ b and b ≥ 10. Based
311 * on equations (26), (27) and (28) in Didonato and Morris (1992).
312 *
313 * @param a First argument.
314 * @param b Second argument.
315 * @return the value of {@code Delta(b) - Delta(a + b)}
316 * @throws MathIllegalArgumentException if {@code a < 0} or {@code a > b}
317 * @throws MathIllegalArgumentException if {@code b < 10}
318 */
319 private static double deltaMinusDeltaSum(final double a,
320 final double b)
321 throws MathIllegalArgumentException {
322
323 MathUtils.checkRangeInclusive(a, 0, b);
324 if (b < 10) {
325 throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_SMALL,
326 b, 10);
327 }
328
329 final double h = a / b;
330 final double p = h / (1.0 + h);
331 final double q = 1.0 / (1.0 + h);
332 final double q2 = q * q;
333 /*
334 * s[i] = 1 + q + ... - q**(2 * i)
335 */
336 final double[] s = new double[DELTA.length];
337 s[0] = 1.0;
338 for (int i = 1; i < s.length; i++) {
339 s[i] = 1.0 + (q + q2 * s[i - 1]);
340 }
341 /*
342 * w = Delta(b) - Delta(a + b)
343 */
344 final double sqrtT = 10.0 / b;
345 final double t = sqrtT * sqrtT;
346 double w = DELTA[DELTA.length - 1] * s[s.length - 1];
347 for (int i = DELTA.length - 2; i >= 0; i--) {
348 w = t * w + DELTA[i] * s[i];
349 }
350 return w * p / b;
351 }
352
353 /**
354 * Returns the value of Δ(p) + Δ(q) - Δ(p + q), with p, q ≥ 10. Based on
355 * the <em>NSWC Library of Mathematics Subroutines</em> double precision
356 * implementation, {@code DBCORR}. In
357 * {@code BetaTest.testSumDeltaMinusDeltaSum()}, this private method is
358 * accessed through reflection.
359 *
360 * @param p First argument.
361 * @param q Second argument.
362 * @return the value of {@code Delta(p) + Delta(q) - Delta(p + q)}.
363 * @throws MathIllegalArgumentException if {@code p < 10.0} or {@code q < 10.0}.
364 */
365 private static double sumDeltaMinusDeltaSum(final double p,
366 final double q) {
367
368 if (p < 10.0) {
369 throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_SMALL,
370 p, 10.0);
371 }
372 if (q < 10.0) {
373 throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_SMALL,
374 q, 10.0);
375 }
376
377 final double a = FastMath.min(p, q);
378 final double b = FastMath.max(p, q);
379 final double sqrtT = 10.0 / a;
380 final double t = sqrtT * sqrtT;
381 double z = DELTA[DELTA.length - 1];
382 for (int i = DELTA.length - 2; i >= 0; i--) {
383 z = t * z + DELTA[i];
384 }
385 return z / a + deltaMinusDeltaSum(a, b);
386 }
387
388 /**
389 * Returns the value of log B(p, q) for 0 ≤ x ≤ 1 and p, q > 0. Based on the
390 * <em>NSWC Library of Mathematics Subroutines</em> implementation,
391 * {@code DBETLN}.
392 *
393 * @param p First argument.
394 * @param q Second argument.
395 * @return the value of {@code log(Beta(p, q))}, {@code NaN} if
396 * {@code p <= 0} or {@code q <= 0}.
397 */
398 public static double logBeta(final double p, final double q) {
399 if (Double.isNaN(p) || Double.isNaN(q) || (p <= 0.0) || (q <= 0.0)) {
400 return Double.NaN;
401 }
402
403 final double a = FastMath.min(p, q);
404 final double b = FastMath.max(p, q);
405 if (a >= 10.0) {
406 final double w = sumDeltaMinusDeltaSum(a, b);
407 final double h = a / b;
408 final double c = h / (1.0 + h);
409 final double u = -(a - 0.5) * FastMath.log(c);
410 final double v = b * FastMath.log1p(h);
411 if (u <= v) {
412 return (((-0.5 * FastMath.log(b) + HALF_LOG_TWO_PI) + w) - u) - v;
413 } else {
414 return (((-0.5 * FastMath.log(b) + HALF_LOG_TWO_PI) + w) - v) - u;
415 }
416 } else if (a > 2.0) {
417 if (b > 1000.0) {
418 final int n = (int) FastMath.floor(a - 1.0);
419 double prod = 1.0;
420 double ared = a;
421 for (int i = 0; i < n; i++) {
422 ared -= 1.0;
423 prod *= ared / (1.0 + ared / b);
424 }
425 return (FastMath.log(prod) - n * FastMath.log(b)) +
426 (Gamma.logGamma(ared) +
427 logGammaMinusLogGammaSum(ared, b));
428 } else {
429 double prod1 = 1.0;
430 double ared = a;
431 while (ared > 2.0) {
432 ared -= 1.0;
433 final double h = ared / b;
434 prod1 *= h / (1.0 + h);
435 }
436 if (b < 10.0) {
437 double prod2 = 1.0;
438 double bred = b;
439 while (bred > 2.0) {
440 bred -= 1.0;
441 prod2 *= bred / (ared + bred);
442 }
443 return FastMath.log(prod1) +
444 FastMath.log(prod2) +
445 (Gamma.logGamma(ared) +
446 (Gamma.logGamma(bred) -
447 logGammaSum(ared, bred)));
448 } else {
449 return FastMath.log(prod1) +
450 Gamma.logGamma(ared) +
451 logGammaMinusLogGammaSum(ared, b);
452 }
453 }
454 } else if (a >= 1.0) {
455 if (b > 2.0) {
456 if (b < 10.0) {
457 double prod = 1.0;
458 double bred = b;
459 while (bred > 2.0) {
460 bred -= 1.0;
461 prod *= bred / (a + bred);
462 }
463 return FastMath.log(prod) +
464 (Gamma.logGamma(a) +
465 (Gamma.logGamma(bred) -
466 logGammaSum(a, bred)));
467 } else {
468 return Gamma.logGamma(a) +
469 logGammaMinusLogGammaSum(a, b);
470 }
471 } else {
472 return Gamma.logGamma(a) +
473 Gamma.logGamma(b) -
474 logGammaSum(a, b);
475 }
476 } else {
477 if (b >= 10.0) {
478 return Gamma.logGamma(a) +
479 logGammaMinusLogGammaSum(a, b);
480 } else {
481 // The following command is the original NSWC implementation.
482 // return Gamma.logGamma(a) +
483 // (Gamma.logGamma(b) - Gamma.logGamma(a + b));
484 // The following command turns out to be more accurate.
485 return FastMath.log(Gamma.gamma(a) * Gamma.gamma(b) /
486 Gamma.gamma(a + b));
487 }
488 }
489 }
490 }