public class Erf extends Object
erf(x)
erfc(x) = 1 - erf(x)
phi(x) = 0.5 * erfc(-x/sqrt(2))
phic(x) = 0.5 * erfc(x/sqrt(2))
Note that phi(x) gives the probability for an observation smaller than x for a Gaussian probability distribution with zero mean and unit standard deviation, while phic(x) gives the probability for an observation larger than x.
The algorithms for erf(x) and erfc(x) are based on Schonfelder's work. See J.L. Schonfelder, "Chebyshev Expansions for the Error and Related Functions", Math. Comp. 32, 1232 (1978). The calculations of phi(x) and phic(x) are trivially calculated using the erfc(x) algorithm.
Schonfelder's algorithms are advertised to provide "30 digit" accurracy. Since this level of accuracy exceeds the machine precision for doubles, summation terms whose relative weight is below machine precision are dropped.
In this algorithm, we calculate
erf(x) = x* y(t) for |x| < 2
erf(x) = 1 - exp(-x*x) * y(t) / x for x >= 2
erfc(|x|) = exp(-x*x)*y(|x|)
The functions y(x) are expanded in terms of Chebyshev polynomials, where there is a different set of coefficients a[r] for each of the above 3 cases.
y(x) = Sum'( a[r] * T(r, t) )
The notation Sum' indicates that the r = 0 term is divided by 2.
The variable t is defined as
t = ( x*x - 2 ) / 2 for erf(x) with x < 2
t = ( 21 - 2*x*x ) / (5 + 2*x*x) for erf(x) with x >= 2
t = ( 4*|x| - 15 ) / ( 4*|x| + 15 ) for erfc(x)
The code and implementation are based on Alan Genz's FORTRAN source code that can be found at http://www.math.wsu.edu/faculty/genz/homepage.
Genz's code was a bit tricky to "reverse engineer", so we go through the way these calculations are performed in some detail. Rather than calculate y(x) directly, he calculates
bm = Sum( a[r] * U(r, t) ) r = 0 : N
bp = Sum( a[r] * U(r-2, t) ) r = 2 : N
where U(r, t) are Chebyshev polynomials of the second kind. The coefficients a[r] decrease with r, and the value of N is chosen where a[N] / a[0] is ~10^-16, reflecting the machine precision for doubles.
The Chebyshev polynomials of the second kind U(r, t) are calculated using the recursion relation:
U(r, t) = 2 * t * U(r-1, t) - U(r-2, t)
Genz uses the identity
T(r, t) = ( U(r, t) - U(r-2, t) ) / 2
to calculate y(x)
y(x) = ( bm - bp ) / 2.
Note that we get the correct contributions for the r = 0 and r = 1 terms by ignoring these terms in the bp sum, including getting the desired factor of 1/2 in the contribution from the r = 0 term.
| Modifier and Type | Field and Description |
|---|---|
private static double[] |
a1 |
private static double[] |
a2 |
private static double[] |
a3 |
private static double |
rtwo |
| Constructor and Description |
|---|
Erf() |
| Modifier and Type | Method and Description |
|---|---|
static double |
erf(double x)
Calculate the error function
|
static double |
erfc(double x)
Calculate the error function complement
|
static double |
phi(double x)
Calcualate the probability for an observation smaller than x for a
Gaussian probability distribution with zero mean and unit standard
deviation
|
static double |
phic(double x)
Calculate the probability for an observation larger than x for a
Gaussian probability distribution with zero mean and unit standard
deviation
|
private static double rtwo
private static double[] a1
private static double[] a2
private static double[] a3
public static double erf(double x)
x - argumentpublic static double erfc(double x)
x - argumentpublic static double phi(double x)
x - argumentpublic static double phic(double x)
x - argumentCopyright © 2016 Linear Collider Detector (LCD). All rights reserved.