1 /*
2 * Copyright (C) 2015 Alberto Irurueta Carro (alberto@irurueta.com)
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
17 package com.irurueta.statistics;
18
19 /**
20 * Defines the gamma function, which is a function that extends the concept of
21 * factorials from natural to real and complex numbers (except zero and negative
22 * integer values).
23 * If argument of gamma function is natural and positive, then it relates to
24 * factorials as: Gamma(n) = (n - 1)!
25 */
26 public class Gamma extends GaussLegendreQuadrature {
27
28 /**
29 * Maximum number of logarithm of factorials cached.
30 */
31 protected static final int MAX_CACHED_LOG_FACTORIALS = 2000;
32
33 /**
34 * Defines when to switch to quadrature method.
35 */
36 private static final int ASWITCH = 100;
37
38 /**
39 * Epsilon for double. It is related to machine precision.
40 */
41 private static final double EPS = Math.ulp(1.0);
42
43 /**
44 * Constant related to machine precision.
45 */
46 private static final double FPMIN = Double.MIN_VALUE / EPS;
47
48 /**
49 * Maximum value supported to estimate factorials for.
50 */
51 private static final int MAX_FACTORIALS = 170;
52
53 /**
54 * Default number of maximum iterations to compute incomplete gamma functions.
55 */
56 private static final int DEFAULT_MAX_ITERATIONS = 100;
57
58 /**
59 * Coefficients for computation of logarithm of gamma function.
60 */
61 private static final double[] COF = {57.1562356658629235, -59.5979603554754912,
62 14.1360979747417471, -0.491913816097620199, .339946499848118887e-4,
63 .465236289270485756e-4, -.983744753048795646e-4, .158088703224912494e-3,
64 -.210264441724104883e-3, .217439618115212643e-3, -.164318106536763890e-3,
65 .844182239838527433e-4, -.261908384015814087e-4, .368991826595316234e-5};
66
67 /**
68 * Indicates whether factorials have been initialized and stored in a table
69 * for future faster access.
70 */
71 private static boolean factorialsInitialized;
72
73 /**
74 * Table where factorials are cached for faster future access.
75 */
76 private static double[] factorialsTable;
77
78 /**
79 * Indicates whether logarithm of factorials has been initialized and
80 * stored in a table for future faster access.
81 */
82 private static boolean logarithmOfFactorialsInitialized;
83
84 /**
85 * Table where logarithms of factorials are cached for faster future access.
86 */
87 private static double[] logarithmOfFactorialsTable;
88
89 /**
90 * Logarithm of gamma function.
91 */
92 private double gln;
93
94 /**
95 * Returns logarithm of gamma function.
96 *
97 * @return logarithm of gamma function.
98 */
99 public double getGln() {
100 return gln;
101 }
102
103 /**
104 * Returns the value ln(gamma(xx)) for xx > 0.
105 *
106 * @param xx a value.
107 * @return the logarithm of gamma function.
108 * @throws IllegalArgumentException if value is negative.
109 */
110 public static double gammln(final double xx) {
111 if (xx <= 0.0) {
112 throw new IllegalArgumentException("bad arg in gammln");
113 }
114
115 int j;
116 double x;
117 double tmp;
118 double y;
119 double ser;
120
121 y = x = xx;
122 tmp = x + 5.24218750000000000;
123 tmp = (x + 0.5) * Math.log(tmp) - tmp;
124 ser = 0.999999999999997092;
125 for (j = 0; j < 14; j++) {
126 ser += COF[j] / ++y;
127 }
128 return tmp + Math.log(2.5066282746310005 * ser / x);
129 }
130
131 /**
132 * Returns the factorial of n (n!) as a floating-point number.
133 * Factorials up to 22! have exact double precision representations.
134 * Factorials from 23! to 170! are approximate due to IEEE double
135 * representation.
136 * Factorials equal or greater than 170! are out of range.
137 *
138 * @param n value to compute factorial for.
139 * @return factorial of n, (i.e. n!).
140 * @throws IllegalArgumentException if provided value generates a factorial
141 * that cannot be represented using double precision.
142 */
143 public static double factrl(final int n) {
144 if (n < 0 || n > MAX_FACTORIALS) {
145 throw new IllegalArgumentException("factrl out of range");
146 }
147
148 if (!factorialsInitialized) {
149 factorialsInitialized = true;
150 factorialsTable = new double[171];
151 factorialsTable[0] = 1.;
152 for (int i = 1; i < (MAX_FACTORIALS + 1); i++) {
153 factorialsTable[i] = i * factorialsTable[i - 1];
154 }
155 }
156
157 return factorialsTable[n];
158 }
159
160 /**
161 * Returns logarithm of n!
162 *
163 * @param n value to compute logarithm of factorial for.
164 * @return logarithm of factorial.
165 * @throws IllegalArgumentException if provided value is negative.
166 */
167 public static double factln(final int n) {
168 if (n < 0) {
169 throw new IllegalArgumentException("negative arg in factln");
170 }
171
172 if (!logarithmOfFactorialsInitialized) {
173 logarithmOfFactorialsInitialized = true;
174 logarithmOfFactorialsTable = new double[MAX_CACHED_LOG_FACTORIALS];
175 for (int i = 0; i < MAX_CACHED_LOG_FACTORIALS; i++) {
176 logarithmOfFactorialsTable[i] = gammln(i + 1.);
177 }
178 }
179 if (n < MAX_CACHED_LOG_FACTORIALS) {
180 return logarithmOfFactorialsTable[n];
181 }
182 return gammln(n + 1.);
183 }
184
185 /**
186 * Returns the binomial coefficient (n k) as a floating-point number, which
187 * indicates the discrete probability distribution of getting exactly k
188 * successes in n trials.
189 *
190 * @param n number of trials.
191 * @param k number of successes.
192 * @return binomial value.
193 * @throws IllegalArgumentException if either n or k are negative or if k is greater than n.
194 */
195 public static double bico(final int n, final int k) {
196 if (k < 0 || k > n) {
197 throw new IllegalArgumentException("bad args in bico");
198 }
199
200 if (n < MAX_FACTORIALS + 1) {
201 return Math.floor(0.5 + factrl(n) / (factrl(k) * factrl(n - k)));
202 }
203
204 // The floor function cleans up round-off error for smaller values of n
205 // and k.
206 return Math.floor(0.5 + Math.exp(factln(n) - factln(k) - factln(n - k)));
207 }
208
209 /**
210 * Returns the value of the beta function B(z, w) (aka Euler integral).
211 *
212 * @param z a parameter.
213 * @param w a parameter.
214 * @return value of beta function.
215 * @throws IllegalArgumentException if either z or w are negative.
216 */
217 public static double beta(final double z, final double w) {
218 return Math.exp(gammln(z) + gammln(w) - gammln(z + w));
219 }
220
221 /**
222 * Returns the incomplete gamma function P(a,x).
223 *
224 * @param a a parameter.
225 * @param x x parameter.
226 * @return value of incomplete gamma function.
227 * @throws IllegalArgumentException if provided values are invalid.
228 * @throws MaxIterationsExceededException if convergence cannot be reached.
229 */
230 public double gammp(final double a, final double x) throws MaxIterationsExceededException {
231 if (x < 0.0 || a <= 0.0) {
232 throw new IllegalArgumentException("bad args in gammp");
233 }
234
235 if (x == 0.0) {
236 return 0.0;
237 } else if ((int) a >= ASWITCH) {
238 return gammpapprox(a, x, 1);
239 } else if (x < a + 1.0) {
240 return gser(a, x);
241 } else {
242 return 1.0 - gcf(a, x);
243 }
244 }
245
246 /**
247 * Returns the incomplete gamma function Q(a, x) = 1 - P(a, x).
248 *
249 * @param a a parameter.
250 * @param x x parameter.
251 * @return value of incomplete gamma function.
252 * @throws IllegalArgumentException if provided values are invalid.
253 * @throws MaxIterationsExceededException if convergence cannot be reached.
254 */
255 public double gammq(final double a, final double x) throws MaxIterationsExceededException {
256 if (x < 0.0 || a <= 0.0) {
257 throw new IllegalArgumentException("bad args in gammq");
258 }
259
260 if (x == 0.0) {
261 return 1.0;
262 } else if ((int) a >= ASWITCH) {
263 return gammpapprox(a, x, 0);
264 } else if (x < a + 1.0) {
265 return 1.0 - gser(a, x);
266 } else {
267 return gcf(a, x);
268 }
269 }
270
271 /**
272 * Returns the incomplete gamma function P(a, x) evaluated by its series
273 * representation.
274 *
275 * @param a a parameter to obtain its gamma logarithm.
276 * @param x x parameter.
277 * @return incomplete gamma function.
278 * @throws MaxIterationsExceededException if maximum number of iterations is
279 * exceeded.
280 */
281 private double gser(final double a, final double x) throws MaxIterationsExceededException {
282 return gser(a, x, DEFAULT_MAX_ITERATIONS);
283 }
284
285 /**
286 * Returns the incomplete gamma function P(a, x) evaluated by its series
287 * representation.
288 *
289 * @param a a parameter to obtain its gamma logarithm.
290 * @param x x parameter.
291 * @param maxIterations maximum number of iterations.
292 * @return incomplete gamma function.
293 * @throws MaxIterationsExceededException if maximum number of iterations is
294 * exceeded.
295 */
296 @SuppressWarnings("SameParameterValue")
297 private double gser(final double a, final double x, final int maxIterations) throws MaxIterationsExceededException {
298 double sum;
299 double del;
300 double ap;
301
302 gln = gammln(a);
303 ap = a;
304 del = sum = 1.0 / a;
305 for (; ; ) {
306 ++ap;
307 del *= x / ap;
308 sum += del;
309 if (Math.abs(del) < Math.abs(sum) * EPS) {
310 return sum * Math.exp(-x + a * Math.log(x) - gln);
311 }
312
313 if (ap >= maxIterations) {
314 throw new MaxIterationsExceededException();
315 }
316 }
317 }
318
319 /**
320 * Returns the incomplete gamma function Q(a, x) evaluated by its continued
321 * fraction representation.
322 *
323 * @param a a parameter to obtain its gamma logarithm.
324 * @param x x parameter.
325 * @return incomplete gamma function.
326 * @throws MaxIterationsExceededException if maximum number of iterations is
327 * exceeded.
328 */
329 private double gcf(final double a, final double x) throws MaxIterationsExceededException {
330 return gcf(a, x, DEFAULT_MAX_ITERATIONS);
331 }
332
333 /**
334 * Returns the incomplete gamma function Q(a, x) evaluated by its continued
335 * fraction representation.
336 *
337 * @param a a parameter to obtain its gamma logarithm.
338 * @param x x parameter.
339 * @param maxIterations maximum number of iterations.
340 * @return incomplete gamma function.
341 * @throws MaxIterationsExceededException if maximum number of iterations is
342 * exceeded.
343 */
344 @SuppressWarnings("SameParameterValue")
345 private double gcf(final double a, final double x, final int maxIterations) throws MaxIterationsExceededException {
346 int i;
347 double an;
348 double b;
349 double c;
350 double d;
351 double del;
352 double h;
353 gln = gammln(a);
354 b = x + 1.0 - a;
355 c = 1.0 / FPMIN;
356 d = 1.0 / b;
357 h = d;
358 for (i = 1; ; i++) {
359 an = -i * (i - a);
360 b += 2.0;
361 d = an * d + b;
362 if (Math.abs(d) < FPMIN) {
363 d = FPMIN;
364 }
365 c = b + an / c;
366 if (Math.abs(c) < FPMIN) {
367 c = FPMIN;
368 }
369 d = 1.0 / d;
370 del = d * c;
371 h *= del;
372 if (Math.abs(del - 1.0) <= EPS) {
373 break;
374 }
375 if (i >= maxIterations) {
376 throw new MaxIterationsExceededException();
377 }
378 }
379 return Math.exp(-x + a * Math.log(x) - gln) * h;
380 }
381
382 /**
383 * Incomplete gamma by quadrature. Returns P(a, x) or Q(a, x), when psig is
384 * 1 or 0 respectively.
385 *
386 * @param a a parameter.
387 * @param x x parameter.
388 * @param psig a flag.
389 * @return incomplete gamma by quadrature.
390 */
391 private double gammpapprox(final double a, final double x, final int psig) {
392 int j;
393 final double xu;
394 double t;
395 double sum;
396 final double ans;
397 final double a1 = a - 1.0;
398 final double lna1 = Math.log(a1);
399 final double sqrta1 = Math.sqrt(a1);
400 gln = gammln(a);
401 if (x > a1) {
402 xu = Math.max(a1 + 11.5 * sqrta1, x + 6.0 * sqrta1);
403 } else {
404 xu = Math.max(0., Math.min(a1 - 7.5 * sqrta1, x - 5.0 * sqrta1));
405 }
406 sum = 0;
407 for (j = 0; j < N_GAU; j++) {
408 t = x + (xu - x) * Y[j];
409 sum += W[j] * Math.exp(-(t - a1) + a1 * (Math.log(t) - lna1));
410 }
411 ans = sum * (xu - x) * Math.exp(a1 * (lna1 - 1.) - gln);
412
413 if (psig != 0) {
414 return ans > 0.0 ? 1.0 - ans : -ans;
415 } else {
416 return ans >= 0.0 ? ans : 1.0 + ans;
417 }
418 }
419
420 /**
421 * Inverse function on x of P(a, x).
422 * Returns x such that P(a,x) = p for an argument p between 0 and 1.
423 *
424 * @param p argument p.
425 * @param a a parameter.
426 * @return inverse value.
427 * @throws IllegalArgumentException if arguments are invalid.
428 * @throws MaxIterationsExceededException if maximum number of iterations is
429 * exceeded.
430 */
431 public double invgammp(final double p, final double a) throws MaxIterationsExceededException {
432 int j;
433 double x;
434 double err;
435 double t;
436 double u;
437 final double pp;
438 double lna1 = 0.0;
439 double afac = 0.0;
440 final double a1 = a - 1.0;
441 gln = gammln(a);
442 if (a <= 0.) {
443 throw new IllegalArgumentException("a must be pos in invgammap");
444 }
445 if (p >= 1.) {
446 return Math.max(100., a + 100. * Math.sqrt(a));
447 }
448 if (p <= 0.) {
449 return 0.0;
450 }
451 if (a > 1.) {
452 lna1 = Math.log(a1);
453 afac = Math.exp(a1 * (lna1 - 1.) - gln);
454 pp = p < 0.5 ? p : 1. - p;
455 t = Math.sqrt(-2. * Math.log(pp));
456 x = (2.30753 + t * 0.27061) / (1. + t * (0.99229 + t * 0.04481)) - t;
457 if (p < 0.5) {
458 x = -x;
459 }
460 x = Math.max(1.e-3, a * Math.pow(1. - 1. / (9. * a) - x / (3. * Math.sqrt(a)), 3));
461 } else {
462 t = 1.0 - a * (0.253 + a * 0.12);
463 if (p < t) {
464 x = Math.pow(p / t, 1. / a);
465 } else {
466 x = 1. - Math.log(1. - (p - t) / (1. - t));
467 }
468 }
469 for (j = 0; j < 12; j++) {
470 if (x <= 0.0) {
471 return 0.0;
472 }
473 err = gammp(a, x) - p;
474 if (a > 1.) {
475 t = afac * Math.exp(-(x - a1) + a1 * (Math.log(x) - lna1));
476 } else {
477 t = Math.exp(-x + a1 * Math.log(x) - gln);
478 }
479 u = err / t;
480 x -= (t = u / (1. - 0.5 * Math.min(1., u * ((a - 1.) / x - 1))));
481 if (x <= 0.) {
482 x = 0.5 * (x + t);
483 }
484 if (Math.abs(t) < EPS * x) {
485 break;
486 }
487 }
488 return x;
489 }
490 }