source: src/main/java/agents/anac/y2019/harddealer/math3/special/Beta.java

Last change on this file was 204, checked in by Katsuhide Fujita, 5 years ago

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