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