View Javadoc
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>&lt;</u> 0.5 that is used is:<br>
245      * &nbsp; &nbsp; erf(z) = (exp(-z<sup>2</sup>) / &radic;&pi;) &Sigma; [ 2z<sup>2n + 1</sup> / (2n + 1)!!]<br>
246      * The Taylor series for erf(z) for abs(z) &gt; 3.7 that is used is:<br>
247      * &nbsp; &nbsp; erf(z) = 1 - (exp(-z<sup>2</sup>) / &radic;&pi;) &Sigma; [ (-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 &lt; z &lt; 3.7 it approximates erf(z) using the following Taylor series:<br>
251      * &nbsp; &nbsp; erf(z) = (2/&radic;&pi;) (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>&lt;</u> 0.5 that is used is:<br>
274      * &nbsp; &nbsp; erf(z) = (exp(-z<sup>2</sup>) / &radic;&pi;) &Sigma; [ 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>&lt;</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      * &nbsp; &nbsp; erf(z) = 1 - (exp(-z<sup>2</sup>) / &radic;&pi;) &Sigma; [ (-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>&gt;</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      * &nbsp; &nbsp; erf(z) = (2/&radic;&pi;) (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 &lt; 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 &lt; 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) = &Gamma;(z) &Gamma;(w) / &Gamma;(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 &lt; 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 }