View Javadoc

1   /*
2    *  Class Erf
3    * 
4    */
5   package org.lcsim.util.probability;
6   
7   /**
8    *
9    * Calculates the following probability integrals:
10   *<p>
11   *   erf(x) <br>
12   *   erfc(x) = 1 - erf(x) <br>
13   *   phi(x) = 0.5 * erfc(-x/sqrt(2)) <br>
14   *   phic(x) = 0.5 * erfc(x/sqrt(2))
15   *<p>
16   * Note that phi(x) gives the probability for an observation smaller than x for
17   * a Gaussian probability distribution with zero mean and unit standard
18   * deviation, while phic(x) gives the probability for an observation larger
19   * than x.
20   *<p>
21   * The algorithms for erf(x) and erfc(x) are based on Schonfelder's work.
22   * See J.L. Schonfelder, "Chebyshev Expansions for the Error and Related
23   * Functions", Math. Comp. 32, 1232 (1978).  The calculations of phi(x)
24   * and phic(x) are trivially calculated using the erfc(x) algorithm.
25   *<p>
26   * Schonfelder's algorithms are advertised to provide "30 digit" accurracy.
27   * Since this level of accuracy exceeds the machine precision for doubles,
28   * summation terms whose relative weight is below machine precision are
29   * dropped.
30   *<p>
31   * In this algorithm, we calculate
32   *<p>
33   *   erf(x) = x* y(t) for |x| < 2 <br>
34   *   erf(x) = 1 - exp(-x*x) * y(t) / x for x >= 2 <br>
35   *   erfc(|x|) = exp(-x*x)*y(|x|) <br>
36   *<p>
37   * The functions y(x) are expanded in terms of Chebyshev polynomials, where
38   * there is a different set of coefficients a[r] for each of the above 3 cases.
39   *<p>
40   *   y(x) = Sum'( a[r] * T(r, t) )
41   *<p>
42   * The notation Sum' indicates that the r = 0 term is divided by 2.
43   *<p>
44   * The variable t is defined as
45   *<p>
46   *   t = ( x*x - 2 ) / 2   for erf(x) with x < 2 <br>
47   *   t = ( 21 - 2*x*x ) / (5 + 2*x*x)   for erf(x) with x >= 2 <br>
48   *   t = ( 4*|x| - 15 ) / ( 4*|x| + 15 )   for erfc(x)
49   *<p>
50   * The code and implementation are based on Alan Genz's FORTRAN source code
51   * that can be found at http://www.math.wsu.edu/faculty/genz/homepage.
52   *<p>
53   * Genz's code was a bit tricky to "reverse engineer", so we go through the
54   * way these calculations are performed in some detail.  Rather than calculate
55   * y(x) directly,  he calculates
56   *<p>
57   *   bm = Sum( a[r] * U(r, t) )    r = 0 : N <br>
58   *   bp = Sum( a[r] * U(r-2, t) )  r = 2 : N
59   *<p>
60   * where U(r, t) are Chebyshev polynomials of the second kind.  The coefficients
61   * a[r] decrease with r, and the value of N is chosen where a[N] / a[0] is
62   * ~10^-16, reflecting the machine precision for doubles.
63   *<p>
64   * The Chebyshev polynomials of the second kind U(r, t) are calculated using the
65   * recursion relation:
66   *<p>
67   *   U(r, t) = 2 * t * U(r-1, t) - U(r-2, t)
68   *<p>
69   * Genz uses the identity
70   *<p>
71   *   T(r, t) = ( U(r, t) - U(r-2, t) ) / 2
72   *<p>
73   * to calculate y(x)
74   *<p>
75   *   y(x) = ( bm - bp ) / 2.
76   *<p>
77   * Note that we get the correct contributions for the r = 0 and r = 1 terms by
78   * ignoring these terms in the bp sum, including getting the desired factor
79   * of 1/2 in the contribution from the r = 0 term.
80   *
81   * @author Richard Partridge
82   */
83  public class Erf {
84  
85      private static double rtwo = 1.414213562373095048801688724209e0;
86  
87      //  Coefficients for the erf(x) calculation with |x| < 2
88      private static double[] a1 = {
89          1.483110564084803581889448079057e0,
90          -3.01071073386594942470731046311e-1,
91          6.8994830689831566246603180718e-2,
92          -1.3916271264722187682546525687e-2,
93          2.420799522433463662891678239e-3,
94          -3.65863968584808644649382577e-4,
95          4.8620984432319048282887568e-5,
96          -5.749256558035684835054215e-6,
97          6.11324357843476469706758e-7,
98          -5.8991015312958434390846e-8,
99          5.207009092068648240455e-9,
100         -4.23297587996554326810e-10,
101         3.1881135066491749748e-11,
102         -2.236155018832684273e-12,
103         1.46732984799108492e-13,
104         -9.044001985381747e-15,
105         5.25481371547092e-16};
106 
107     //  Coefficients for the err(x) calculation with x > 2
108     private static double[] a2 = {
109       1.077977852072383151168335910348e0,
110       -2.6559890409148673372146500904e-2,
111       -1.487073146698099509605046333e-3,
112       -1.38040145414143859607708920e-4,
113       -1.1280303332287491498507366e-5,
114       -1.172869842743725224053739e-6,
115       -1.03476150393304615537382e-7,
116       -1.1899114085892438254447e-8,
117       -1.016222544989498640476e-9,
118       -1.37895716146965692169e-10,
119       -9.369613033737303335e-12,
120       -1.918809583959525349e-12,
121       -3.7573017201993707e-14,
122       -3.7053726026983357e-14,
123       2.627565423490371e-15,
124       -1.121322876437933e-15,
125       1.84136028922538e-16};
126 
127     //  Coefficients for the erfc(x) calculation
128     private static double[] a3 = {
129         6.10143081923200417926465815756e-1,
130         -4.34841272712577471828182820888e-1,
131         1.76351193643605501125840298123e-1,
132         -6.0710795609249414860051215825e-2,
133         1.7712068995694114486147141191e-2,
134         -4.321119385567293818599864968e-3,
135         8.54216676887098678819832055e-4,
136         -1.27155090609162742628893940e-4,
137         1.1248167243671189468847072e-5,
138         3.13063885421820972630152e-7,
139         -2.70988068537762022009086e-7,
140         3.0737622701407688440959e-8,
141         2.515620384817622937314e-9,
142         -1.028929921320319127590e-9,
143         2.9944052119949939363e-11,
144         2.6051789687266936290e-11,
145         -2.634839924171969386e-12,
146         -6.43404509890636443e-13,
147         1.12457401801663447e-13,
148         1.7281533389986098e-14,
149         -4.264101694942375e-15,
150         -5.45371977880191e-16,
151         1.58697607761671e-16,
152         2.0899837844334e-17,
153         -5.900526869409e-18};
154 
155     /**
156      * Calculate the error function
157      * 
158      * @param x argument
159      * @return error function
160      */
161     public static double erf(double x) {
162 
163         //  Initialize
164         double xa = Math.abs(x);
165         double erf;
166 
167         //  Case 1: |x| < 2
168         if (xa < 2.) {
169 
170             //  First calculate 2*t
171             double tt = x*x - 2.;
172 
173             //  Initialize the recursion variables.
174             double bm = 0.;
175             double b = 0.;
176             double bp = 0.;
177 
178             //  Calculate bm and bp as defined above
179             for (int i = 16; i >= 0; i--) {
180                 bp = b;
181                 b = bm;
182                 bm = tt * b - bp + a1[i];
183             }
184 
185             //  Finally, calculate erfc using the Chebyshev polynomial identity
186             erf = x * (bm - bp) / 2.;
187 
188         } else {
189 
190             //  Case 2: |x| >= 2
191 
192             //  First calculate 2*t
193             double tt = (42. - 4 * xa*xa) / (5. + 2 * xa*xa);
194 
195             //  Initialize the recursion variables.
196             double bm = 0.;
197             double b = 0.;
198             double bp = 0.;
199 
200             //  Calculate bm and bp as defined above
201             for (int i = 16; i >= 0; i--) {
202                 bp = b;
203                 b = bm;
204                 bm = tt * b - bp + a2[i];
205             }
206 
207             //  Finally, calculate erfc using the Chebyshev polynomial identity
208             erf = 1. - Math.exp(-x * x) * (bm - bp) / (2. * xa);
209 
210             //  Take care of negative argument for case 2
211             if (x < 0.) erf = -erf;
212         }
213 
214         //  Finished both cases!
215         return erf;
216     }
217 
218     /**
219      * Calculate the error function complement
220      * @param x argument
221      * @return error function complement
222      */
223     public static double erfc(double x) {
224 
225         //  Initialize
226         double xa = Math.abs(x);
227         double erfc;
228 
229         //  Set phi to 0 when the argument is too big
230         if (xa > 100.) {
231             erfc = 0.;
232         } else {
233 
234             //  First calculate 2*t
235             double tt = (8. * xa - 30.) / (4. * xa + 15.);
236 
237             //  Initialize the recursion variables.
238             double bm = 0.;
239             double b = 0.;
240             double bp = 0.;
241 
242             //  Calculate bm and bp as defined above
243             for (int i = 24; i >= 0; i--) {
244                 bp = b;
245                 b = bm;
246                 bm = tt * b - bp + a3[i];
247             }
248 
249             //  Finally, calculate erfc using the Chebyshev polynomial identity
250             erfc = Math.exp(-x * x) * (bm - bp) / 2.;
251         }
252 
253         //  Cacluate erfc for negative arguments
254         if (x < 0.) erfc = 2. - erfc;
255 
256         return erfc;
257     }
258 
259     /**
260      * Calcualate the probability for an observation smaller than x for a
261      * Gaussian probability distribution with zero mean and unit standard
262      * deviation
263      * 
264      * @param x argument
265      * @return probability integral
266      */
267     public static double phi(double x) {
268         return 0.5 * erfc( -x / rtwo);
269     }
270 
271     /**
272      * Calculate the probability for an observation larger than x for a
273      * Gaussian probability distribution with zero mean and unit standard
274      * deviation
275      *
276      * @param x argument
277      * @return probability integral
278      */
279     public static double phic(double x) {
280         return 0.5 * erfc(x / rtwo);
281     }
282 }
283