My goal is to implement a CDF (cumulative distribution function) for the Standard Normal Distribution and to do so with the best possible numeric accuracy.
This is my starting point: (1 + erf(x / sqrt(2))) / 2
After reading John Cook's article about why the erf() and erfc() functions are both necessary, I also tried this: (2 - erfc(x / sqrt(2))) / 2
The result of 1 + erf(-2) differs from 2 - erfc(-2). Presumably the second is more accurate.
Assuming that various math libraries provide high quality implementations of sqrt, erf, and erfc, does anyone know how to find the range of values where the second formula becomes more accurate than the first?
For negative inputs, erfc(-x) is superior 1 + erf(x) because it avoids the loss of precision at the addition step. The improvement quickly becomes massive as x becomes more negative.
For positive inputs, erfc(-x) is often better than 1 + erf(x) but the difference is only 1 ULP (unit in the last place) and not really of consequence.
Conclusion: It is reasonable to use erfc(-x) for all input values. So, a numerically favorable implementation of the CDF for the Standard Normal Distribution is 0.5 * erfc(-x / sqrt(2))
Here is a Python script to compare the two approaches.
For 10,000 random values in each range bin, it shows the maximum difference in between the two approaches as measured in ULPs. It also shows how often each technique beats or agrees with the other as compared to a reference implementation using mpmath:
import math
from mpmath import mp # Get this with: pip install mpmath
from collections import Counter
from random import uniform
from itertools import pairwise
mp.dps = 50
def best(x):
'Which is better? erf, erfc, or are they the same?'
ref = 1 + mp.erf(x)
e = 1 + math.erf(x)
c = math.erfc(-x)
de = abs(e - ref)
dc = abs(c - ref)
return 'erf' if de < dc else 'erfc' if dc < de else '='
def err(x):
'Diffence in ulp'
e = 1 + math.erf(x)
c = math.erfc(-x)
return abs(round((e - c) / math.ulp(c)))
def frange(lo, hi, steps=10):
step = (hi - lo) / steps
return [lo + i * step for i in range(steps+1)]
for lo, hi in pairwise(frange(-5, 5, steps=20)):
xarr = [uniform(lo, hi) for i in range(10_000)]
winners = Counter(map(best, xarr)).most_common()
max_err = max(map(err, xarr))
print(f'{lo:-5.2f} to {hi:-5.2f}: {max_err:12d} ulp', winners)
The results depend on the underlying math library implementations of erf and erfc . On macOS with clang-1600.0.26.6, this outputs:
-5.00 to -4.50: 409383042002 ulp [('erfc', 10000)]
-4.50 to -4.00: 3202234149 ulp [('erfc', 10000)]
-4.00 to -3.50: 25164371 ulp [('erfc', 10000)]
-3.50 to -3.00: 786399 ulp [('erfc', 10000)]
-3.00 to -2.50: 24489 ulp [('erfc', 9995), ('erf', 5)]
-2.50 to -2.00: 1532 ulp [('erfc', 9924), ('erf', 41), ('=', 35)]
-2.00 to -1.50: 95 ulp [('erfc', 9507), ('erf', 311), ('=', 182)]
-1.50 to -1.00: 11 ulp [('erfc', 7508), ('=', 1269), ('erf', 1223)]
-1.00 to -0.50: 2 ulp [('=', 4107), ('erfc', 4093), ('erf', 1800)]
-0.50 to 0.00: 1 ulp [('=', 9782), ('erfc', 144), ('erf', 74)]
0.00 to 0.50: 1 ulp [('=', 9211), ('erf', 727), ('erfc', 62)]
0.50 to 1.00: 1 ulp [('=', 9536), ('erfc', 377), ('erf', 87)]
1.00 to 1.50: 1 ulp [('=', 8761), ('erfc', 1113), ('erf', 126)]
1.50 to 2.00: 1 ulp [('=', 8749), ('erfc', 1217), ('erf', 34)]
2.00 to 2.50: 1 ulp [('=', 8773), ('erfc', 1224), ('erf', 3)]
2.50 to 3.00: 1 ulp [('=', 8701), ('erfc', 1299)]
3.00 to 3.50: 1 ulp [('=', 8772), ('erfc', 1228)]
3.50 to 4.00: 1 ulp [('=', 8739), ('erfc', 1261)]
4.00 to 4.50: 1 ulp [('=', 8760), ('erfc', 1240)]
4.50 to 5.00: 1 ulp [('=', 8743), ('erfc', 1257)]