1 package nl.tudelft.simulation.jstats.math;
2
3 import org.djutils.exceptions.Throw;
4
5 /**
6 * The ProbMath class defines some very basic probabilistic mathematical functions.
7 * <p>
8 * Copyright (c) 2004-2024 Delft University of Technology, Jaffalaan 5, 2628 BX Delft, the Netherlands. All rights reserved. See
9 * for project information <a href="https://simulation.tudelft.nl/" target="_blank"> https://simulation.tudelft.nl</a>. The DSOL
10 * project is distributed under a three-clause BSD-style license, which can be found at
11 * <a href="https://https://simulation.tudelft.nl/dsol/docs/latest/license.html" target="_blank">
12 * https://https://simulation.tudelft.nl/dsol/docs/latest/license.html</a>.
13 * </p>
14 * @author <a href="https://www.tudelft.nl/averbraeck" target="_blank"> Alexander Verbraeck</a>
15 */
16 public final class ProbMath
17 {
18 /** stored values of n! as a double value. */
19 public static final double[] FACTORIAL_DOUBLE;
20
21 /** stored values of n! as a long value. */
22 public static final long[] FACTORIAL_LONG;
23
24 /** stored values of a^n as a double value. */
25 public static final double[] POW2_DOUBLE;
26
27 /** stored values of 2^n as a long value. */
28 public static final long[] POW2_LONG;
29
30 /** calculate n! and 2^n. */
31 static
32 {
33 FACTORIAL_DOUBLE = new double[171];
34 FACTORIAL_DOUBLE[0] = 1.0;
35 double d = 1.0;
36 for (int i = 1; i <= 170; i++)
37 {
38 d = d * i;
39 FACTORIAL_DOUBLE[i] = d;
40 }
41
42 FACTORIAL_LONG = new long[21];
43 FACTORIAL_LONG[0] = 1L;
44 long n = 1L;
45 for (int i = 1; i <= 20; i++)
46 {
47 n = n * i;
48 FACTORIAL_LONG[i] = n;
49 }
50
51 POW2_DOUBLE = new double[1024];
52 POW2_DOUBLE[0] = 1.0;
53 for (int i = 1; i < 1024; i++)
54 {
55 POW2_DOUBLE[i] = 2.0 * POW2_DOUBLE[i - 1];
56 }
57
58 POW2_LONG = new long[63];
59 POW2_LONG[0] = 1L;
60 for (int i = 1; i < 63; i++)
61 {
62 POW2_LONG[i] = 2L * POW2_LONG[i - 1];
63 }
64 }
65
66 /**
67 * Utility class.
68 */
69 private ProbMath()
70 {
71 // unreachable code for the utility class
72 }
73
74 /**
75 * Compute n factorial as a double. fac(n) = n * (n-1) * (n-2) * ... 2 * 1.
76 * @param n int; the value to calculate n! for
77 * @return double; n factorial
78 */
79 public static double factorial(final int n)
80 {
81 Throw.when(n < 0, IllegalArgumentException.class, "n! with n<0 is invalid");
82 Throw.when(n > 170, IllegalArgumentException.class, "n! with n>170 is larger than Double.MAX_VALUE");
83 return FACTORIAL_DOUBLE[n];
84 }
85
86 /**
87 * Compute n factorial as a double. fac(n) = n * (n-1) * (n-2) * ... 2 * 1.
88 * @param n long; the value to calculate n! for
89 * @return double; n factorial
90 */
91 public static double factorial(final long n)
92 {
93 Throw.when(n < 0, IllegalArgumentException.class, "n! with n<0 is invalid");
94 Throw.when(n > 170, IllegalArgumentException.class, "n! with n>170 is larger than Double.MAX_VALUE");
95 return FACTORIAL_DOUBLE[(int) n];
96 }
97
98 /**
99 * Compute n factorial as a long. fac(n) = n * (n-1) * (n-2) * ... 2 * 1.
100 * @param n int; the value to calculate n! for
101 * @return double; n factorial
102 */
103 public static long fac(final int n)
104 {
105 Throw.when(n < 0, IllegalArgumentException.class, "n! with n<0 is invalid");
106 Throw.when(n > 20, IllegalArgumentException.class, "n! with n>20 is more than 64 bits long");
107 return FACTORIAL_LONG[n];
108 }
109
110 /**
111 * Compute n factorial as a long. fac(n) = n * (n-1) * (n-2) * ... 2 * 1.
112 * @param n long; the value to calculate n! for
113 * @return double; n factorial
114 */
115 public static long fac(final long n)
116 {
117 Throw.when(n < 0, IllegalArgumentException.class, "n! with n<0 is invalid");
118 Throw.when(n > 20, IllegalArgumentException.class, "n! with n>20 is more than 64 bits long");
119 return FACTORIAL_LONG[(int) n];
120 }
121
122 /**
123 * Compute the k-permutations of n.
124 * @param n int; the first parameter
125 * @param k int; the second parameter
126 * @return the k-permutations of n
127 */
128 public static double permutations(final int n, final int k)
129 {
130 if (k > n)
131 {
132 throw new IllegalArgumentException("permutations of (n,k) with k>n");
133 }
134 return factorial(n) / factorial(n - k);
135 }
136
137 /**
138 * Compute the k-permutations of n.
139 * @param n long; the first parameter
140 * @param k long; the second parameter
141 * @return the k-permutations of n
142 */
143 public static double permutations(final long n, final long k)
144 {
145 if (k > n)
146 {
147 throw new IllegalArgumentException("permutations of (n,k) with k>n");
148 }
149 return factorial(n) / factorial(n - k);
150 }
151
152 /**
153 * Computes the k-permutations of n as a long.
154 * @param n int; the first parameter
155 * @param k int; the second parameter
156 * @return the k-permutations of n
157 */
158 public static long perm(final int n, final int k)
159 {
160 if (k > n)
161 {
162 throw new IllegalArgumentException("permutations of (n,k) with k>n");
163 }
164 return fac(n) / fac(n - k);
165 }
166
167 /**
168 * Computes the k-permutations of n as a long.
169 * @param n long; the first parameter
170 * @param k long; the second parameter
171 * @return the k-permutations of n
172 */
173 public static long perm(final long n, final long k)
174 {
175 if (k > n)
176 {
177 throw new IllegalArgumentException("permutations of (n,k) with k>n");
178 }
179 return fac(n) / fac(n - k);
180 }
181
182 /**
183 * computes the combinations of n over k.
184 * @param n int; the first parameter
185 * @param k int; the second parameter
186 * @return the combinations of n over k
187 */
188 public static double combinations(final int n, final int k)
189 {
190 if (k > n)
191 {
192 throw new IllegalArgumentException("combinations of (n,k) with k>n");
193 }
194 return factorial(n) / (factorial(k) * factorial(n - k));
195 }
196
197 /**
198 * computes the combinations of n over k.
199 * @param n long; the first parameter
200 * @param k long; the second parameter
201 * @return the combinations of n over k
202 */
203 public static double combinations(final long n, final long k)
204 {
205 if (k > n)
206 {
207 throw new IllegalArgumentException("combinations of (n,k) with k>n");
208 }
209 return factorial(n) / (factorial(k) * factorial(n - k));
210 }
211
212 /**
213 * computes the combinations of n over k as a long.
214 * @param n int; the first parameter
215 * @param k int; the second parameter
216 * @return the combinations of n over k
217 */
218 public static long comb(final int n, final int k)
219 {
220 if (k > n)
221 {
222 throw new IllegalArgumentException("combinations of (n,k) with k>n");
223 }
224 return fac(n) / (fac(k) * fac(n - k));
225 }
226
227 /**
228 * computes the combinations of n over k as a long.
229 * @param n long; the first parameter
230 * @param k long; the second parameter
231 * @return the combinations of n over k
232 */
233 public static long comb(final long n, final long k)
234 {
235 if (k > n)
236 {
237 throw new IllegalArgumentException("combinations of (n,k) with k>n");
238 }
239 return fac(n) / (fac(k) * fac(n - k));
240 }
241
242 /**
243 * Approximates erf(z) using a Taylor series.<br>
244 * The Taylor series for erf(z) for abs(z) <u><</u> 0.5 that is used is:<br>
245 * erf(z) = (exp(-z<sup>2</sup>) / √π) Σ [ 2z<sup>2n + 1</sup> / (2n + 1)!!]<br>
246 * The Taylor series for erf(z) for abs(z) > 3.7 that is used is:<br>
247 * erf(z) = 1 - (exp(-z<sup>2</sup>) / √π) Σ [ (-1)<sup>n</sup> (2n - 1)!! z<sup>-(2n +
248 * 1)</sup> / 2<sup>n</sup>]<br>
249 * See <a href="https://mathworld.wolfram.com/Erf.html">https://mathworld.wolfram.com/Erf.html</a>. <br>
250 * For 0.5 < z < 3.7 it approximates erf(z) using the following Taylor series:<br>
251 * erf(z) = (2/√π) (z - z<sup>3</sup>/3 + z<sup>5</sup>/10 - z<sup>7</sup>/42 + z<sup>9</sup>/216 -
252 * ...)<br>
253 * The factors are given by <a href="https://oeis.org/A007680">https://oeis.org/A007680</a>, which evaluates to a(n) =
254 * (2n+1)n!. See <a href="https://en.wikipedia.org/wiki/Error_function">https://en.wikipedia.org/wiki/Error_function</a>.
255 * @param z double; the value to calculate the error function for
256 * @return erf(z)
257 */
258 public static double erf(final double z)
259 {
260 double zpos = Math.abs(z);
261 if (zpos < 0.5)
262 {
263 return erfSmall(z);
264 }
265 if (zpos > 3.8)
266 {
267 return erfBig(z);
268 }
269 return erfTaylor(z);
270 }
271
272 /**
273 * The Taylor series for erf(z) for abs(z) <u><</u> 0.5 that is used is:<br>
274 * erf(z) = (exp(-z<sup>2</sup>) / √π) Σ [ 2z<sup>2n + 1</sup> / (2n + 1)!!]<br>
275 * where the !! operator is the 'double factorial' operator which is (n).(n-2)...8.4.2 for even n, and (n).(n-2)...3.5.1 for
276 * odd n. See <a href="https://mathworld.wolfram.com/Erf.html">https://mathworld.wolfram.com/Erf.html</a> formula (9) and
277 * (10). This function would work well for z <u><</u> 0.5.
278 * @param z double; the parameter
279 * @return double; erf(x)
280 */
281 private static double erfSmall(final double z)
282 {
283 double zpos = Math.abs(z);
284 // @formatter:off
285 double sum = zpos
286 + 2.0 * Math.pow(zpos, 3) / 3.0
287 + 4.0 * Math.pow(zpos, 5) / 15.0
288 + 8.0 * Math.pow(zpos, 7) / 105.0
289 + 16.0 * Math.pow(zpos, 9) / 945.0
290 + 32.0 * Math.pow(zpos, 11) / 10395.0
291 + 64.0 * Math.pow(zpos, 13) / 135135.0
292 + 128.0 * Math.pow(zpos, 15) / 2027025.0
293 + 256.0 * Math.pow(zpos, 17) / 34459425.0
294 + 512.0 * Math.pow(zpos, 19) / 654729075.0
295 + 1024.0 * Math.pow(zpos, 21) / 13749310575.0;
296 // @formatter:on
297 return Math.signum(z) * sum * 2.0 * Math.exp(-zpos * zpos) / Math.sqrt(Math.PI);
298 }
299
300 /**
301 * Calculate erf(z) for large values using the Taylor series:<br>
302 * erf(z) = 1 - (exp(-z<sup>2</sup>) / √π) Σ [ (-1)<sup>n</sup> (2n - 1)!! z<sup>-(2n +
303 * 1)</sup> / 2<sup>n</sup>]<br>
304 * where the !! operator is the 'double factorial' operator which is (n).(n-2)...8.4.2 for even n, and (n).(n-2)...3.5.1 for
305 * odd n. See <a href="https://mathworld.wolfram.com/Erf.html">https://mathworld.wolfram.com/Erf.html</a> formula (18) to
306 * (20). This function would work well for z <u>></u> 3.7.
307 * @param z double; the argument
308 * @return double; erf(z)
309 */
310 private static double erfBig(final double z)
311 {
312 double zpos = Math.abs(z);
313 // @formatter:off
314 double sum = 1.0 / zpos
315 - (1.0 / 2.0) * Math.pow(zpos, -3)
316 + (3.0 / 4.0) * Math.pow(zpos, -5)
317 - (15.0 / 8.0) * Math.pow(zpos, -7)
318 + (105.0 / 16.0) * Math.pow(zpos, -9)
319 - (945.0 / 32.0) * Math.pow(zpos, -11)
320 + (10395.0 / 64.0) * Math.pow(zpos, -13)
321 - (135135.0 / 128.0) * Math.pow(zpos, -15)
322 + (2027025.0 / 256.0) * Math.pow(zpos, -17);
323 // @formatter:on
324 return Math.signum(z) * (1.0 - sum * Math.exp(-zpos * zpos) / Math.sqrt(Math.PI));
325 }
326
327 /**
328 * Calculate erf(z) using the Taylor series:<br>
329 * erf(z) = (2/√π) (z - z<sup>3</sup>/3 + z<sup>5</sup>/10 - z<sup>7</sup>/42 + z<sup>9</sup>/216 -
330 * ...)<br>
331 * The factors are given by <a href="https://oeis.org/A007680">https://oeis.org/A007680</a>, which evaluates to a(n) =
332 * (2n+1)n!. See <a href="https://en.wikipedia.org/wiki/Error_function">https://en.wikipedia.org/wiki/Error_function</a>.
333 * This works pretty well on the interval [0.5,3.7].
334 * @param z double; the argument
335 * @return double; erf(z)
336 */
337 private static double erfTaylor(final double z)
338 {
339 double zpos = Math.abs(z);
340 double d = zpos;
341 double zpow = zpos;
342 for (int i = 1; i < 64; i++) // max 64 steps
343 {
344 // calculate Math.pow(zpos, 2 * i + 1) / ((2 * i + 1) * factorial(i));
345 zpow *= zpos * zpos;
346 double term = zpow / ((2.0 * i + 1.0) * ProbMath.factorial(i));
347 d += term * ((i & 1) == 0 ? 1 : -1);
348 if (term < 1E-16)
349 {
350 break;
351 }
352 }
353 return Math.signum(z) * d * 2.0 / Math.sqrt(Math.PI);
354 }
355
356 /**
357 * Approximates erf<sup>-1</sup>(p) based on
358 * <a href="http://www.naic.edu/~jeffh/inverse_cerf.c">http://www.naic.edu/~jeffh/inverse_cerf.c</a> code.
359 * @param y double; the cumulative probability to calculate the inverse error function for
360 * @return erf<sup>-1</sup>(p)
361 */
362 public static double erfInv(final double y)
363 {
364 double ax, t, ret;
365 ax = Math.abs(y);
366
367 /*
368 * This approximation, taken from Table 10 of Blair et al., is valid for |x|<=0.75 and has a maximum relative error of
369 * 4.47 x 10^-8.
370 */
371 if (ax <= 0.75)
372 {
373
374 double[] p = new double[] {-13.0959967422, 26.785225760, -9.289057635};
375 double[] q = new double[] {-12.0749426297, 30.960614529, -17.149977991, 1.00000000};
376
377 t = ax * ax - 0.75 * 0.75;
378 ret = ax * (p[0] + t * (p[1] + t * p[2])) / (q[0] + t * (q[1] + t * (q[2] + t * q[3])));
379
380 }
381 else if (ax >= 0.75 && ax <= 0.9375)
382 {
383 double[] p = new double[] {-.12402565221, 1.0688059574, -1.9594556078, .4230581357};
384 double[] q = new double[] {-.08827697997, .8900743359, -2.1757031196, 1.0000000000};
385
386 /*
387 * This approximation, taken from Table 29 of Blair et al., is valid for .75<=|x|<=.9375 and has a maximum relative
388 * error of 4.17 x 10^-8.
389 */
390 t = ax * ax - 0.9375 * 0.9375;
391 ret = ax * (p[0] + t * (p[1] + t * (p[2] + t * p[3]))) / (q[0] + t * (q[1] + t * (q[2] + t * q[3])));
392
393 }
394 else if (ax >= 0.9375 && ax <= (1.0 - 1.0e-9))
395 {
396 double[] p =
397 new double[] {.1550470003116, 1.382719649631, .690969348887, -1.128081391617, .680544246825, -.16444156791};
398 double[] q = new double[] {.155024849822, 1.385228141995, 1.000000000000};
399
400 /*
401 * This approximation, taken from Table 50 of Blair et al., is valid for .9375<=|x|<=1-10^-100 and has a maximum
402 * relative error of 2.45 x 10^-8.
403 */
404 t = 1.0 / Math.sqrt(-Math.log(1.0 - ax));
405 ret = (p[0] / t + p[1] + t * (p[2] + t * (p[3] + t * (p[4] + t * p[5])))) / (q[0] + t * (q[1] + t * (q[2])));
406 }
407 else
408 {
409 ret = Double.POSITIVE_INFINITY;
410 }
411
412 return Math.signum(y) * ret;
413 }
414
415 /** Coefficients for the ln(gamma(x)) function. */
416 private static final double[] GAMMALN_COF = {76.18009172947146, -86.50532032941677, 24.01409824083091, -1.231739572450155,
417 0.1208650973866179e-2, -0.5395239384953e-5};
418
419 /**
420 * Calculates ln(gamma(x)). Java version of gammln function in Numerical Recipes in C, p.214.
421 * @param xx double; the value to calculate the gamma function for
422 * @return double; gamma(x)
423 * @throws IllegalArgumentException when x is < 0
424 */
425 public static double gammaln(final double xx)
426 {
427 Throw.when(xx < 0, IllegalArgumentException.class, "gamma function not defined for real values < 0");
428 double x, y, tmp, ser;
429 x = xx;
430 y = x;
431 tmp = x + 5.5;
432 tmp -= (x + 0.5) * Math.log(tmp);
433 ser = 1.000000000190015;
434 for (int j = 0; j <= 5; j++)
435 {
436 ser += GAMMALN_COF[j] / ++y;
437 }
438 return -tmp + Math.log(2.5066282746310005 * ser / x);
439 }
440
441 /**
442 * Calculates gamma(x). Based on the gammln function in Numerical Recipes in C, p.214.
443 * @param x double; the value to calculate the gamma function for
444 * @return double; gamma(x)
445 * @throws IllegalArgumentException when x is < 0
446 */
447 public static double gamma(final double x)
448 {
449 return Math.exp(gammaln(x));
450 }
451
452 /**
453 * Calculates Beta(z, w) where Beta(z, w) = Γ(z) Γ(w) / Γ(z + w).
454 * @param z double; beta function parameter 1
455 * @param w ; beta function parameter 2
456 * @return double; beta(z, w)
457 * @throws IllegalArgumentException when one of the parameters is < 0
458 */
459 public static double beta(final double z, final double w)
460 {
461 Throw.when(z < 0 || w < 0, IllegalArgumentException.class, "beta function not defined for negative arguments");
462 return Math.exp(gammaln(z) + gammaln(w) - gammaln(z + w));
463 }
464
465 }