The normal distribution¶
The standard normal distribution, N(x), the chance a standard normal is below x, by Abramowitz & Stegun 26.2.17: for z >= 0, 1 - N(z) = phi(z) (b1 t + b2 t^2 + b3 t^3 + b4 t^4 + b5 t^5), t = 1 / (1 + p z), with |error| < 7.5e-8 (the formula's own bound, absolute; 7.45e-8 the most found against scipy's ndtr over [-12, 12]).
Below zero it is that upper tail at -x itself, never 1 - N(-x), so a small N(x) is not lost to the rounding of 1; its relative error is the formula's, 3e-7 within one standard deviation, 5e-5 within three, 1.6e-3 within five. Above zero, N(x) is 1 less the tail, good to binary32's rounding of 1 (2^-24 absolute). The six numbers are the approximation's, written as the source prints them. Source: Abramowitz & Stegun, Handbook of Mathematical Functions (1964), 26.2.17 (Zelen & Severo); its coefficients as Wikipedia's "Normal distribution" quotes them (A&S itself not opened). Checked against a reference in numpy, which has its own tests.
- Source
- Abramowitz & Stegun, Handbook of Mathematical Functions (1964), 26.2.17 (Zelen & Severo); its coefficients as Wikipedia's "Normal distribution" quotes them (A&S itself not opened).
- Proof
- reference, invariant, limit.
normal_cdf.trait
namespace lusine.finance;
float normal_cdf(float x)
{
let z = x < 0.0 ? -x : x;
let t = 1.0 / (1.0 + 0.2316419 * z);
let poly = t * (0.319381530 + t * (-0.356563782 + t * (1.781477937 + t * (-1.821255978 + t * 1.330274429))));
let upper = normal_pdf(x = z) * poly;
return x < 0.0 ? upper : 1.0 - upper;
}
Uses, without proving it here: finance.normal_pdf.
N(x) by Abramowitz & Stegun 26.2.17, |error| < 7.5e-8, its lower tail taken directly.
import os, sys
root = os.getcwd()
while not os.path.isfile(os.path.join(root, 'docsite.json')):
root = os.path.dirname(root)
sys.path.insert(0, os.path.join(root, 'tools'))
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import ndtr
from scipy.stats import norm
from lusine import trait_run
EPS = float(np.finfo(np.float32).eps) # Trait's float is binary32
YEAR = 31557600.0 # lusine's year row, the Julian year in seconds: the host speaks base units
The kernel¶
The law as the library states it, and as this page runs it.
print(trait_run.source('finance/normal_cdf.trait'))
// lusine.finance.normal_cdf: the standard normal distribution, N(x), the chance a standard normal is below x, by
// Abramowitz & Stegun 26.2.17: for z >= 0, 1 - N(z) = phi(z) (b1 t + b2 t^2 + b3 t^3 + b4 t^4 + b5 t^5), t = 1 / (1 +
// p z), with |error| < 7.5e-8 (the formula's own bound, absolute; 7.45e-8 the most found against scipy's ndtr over
// [-12, 12]). Below zero it is that upper tail at -x itself, never 1 - N(-x), so a small N(x) is not lost to the
// rounding of 1; its relative error is the formula's, 3e-7 within one standard deviation, 5e-5 within three, 1.6e-3
// within five. Above zero, N(x) is 1 less the tail, good to binary32's rounding of 1 (2^-24 absolute). The six numbers
// are the approximation's, written as the source prints them.
// Source: Abramowitz & Stegun, Handbook of Mathematical Functions (1964), 26.2.17 (Zelen & Severo); its coefficients as
// Wikipedia's "Normal distribution" quotes them (A&S itself not opened).
// Checked against a reference in numpy, which has its own tests.
namespace lusine.finance;
float normal_cdf(float x)
{
let z = x < 0.0 ? -x: x;
let t = 1.0 / (1.0 + 0.2316419 * z);
let poly = t * (0.319381530 + t * (-0.356563782 + t *));
let upper = normal_pdf(x = z) * poly;
return x < 0.0 ? upper: 1.0 - upper;
}
The program¶
The kernel over a lattice of items, one an x; the independent account is scipy's exact ndtr at the binary32 x the device holds. The coefficients are A&S 26.2.17's as Wikipedia's "Normal distribution" quotes them (A&S itself was not opened); the bound is measured here.
print(trait_run.law_program('finance/normal_cdf.trait', 'normal_cdf'))
def cdf(x):
return trait_run.law('finance/normal_cdf.trait', 'normal_cdf', x=x)
namespace probe;
lusine.version version = { { number = 1 } };
struct item {
host float x;
float result_0;
};
lattice<item, 1> items(spacing = 1 m, boundary = wall);
void main() { items.result_0 = lusine.finance.normal_cdf(items.x); }
Reference: scipy's distribution¶
Over [-9, 9] each N(x) is within A&S's 7.5e-8 and two EPS of binary32's rounding of scipy's: the formula's error and the float's.
x = np.linspace(-9.0, 9.0, 721)
xs = np.float32(x).astype(float)
y = cdf(x)
err = y - ndtr(xs)
assert np.max(np.abs(err)) <= 7.5e-8 + 2 * EPS, np.max(np.abs(err))
print(f'largest error {np.max(np.abs(err)):.2e}')
fig, ax = plt.subplots(figsize=(7, 3.5))
ax.plot(xs, err)
ax.set_xlabel('x')
ax.set_ylabel('N(x) - scipy')
plt.show()
largest error 1.64e-07
Invariant: N(-x) = 1 - N(x), rising¶
Off nothing, the two tails add to one within an EPS, and N never falls as x rises.
x = np.linspace(0.01, 9.0, 900)
hi, lo = cdf(x), cdf(-x)
assert np.max(np.abs(hi + lo - 1)) <= EPS, np.max(np.abs(hi + lo - 1))
assert np.all(np.diff(cdf(np.linspace(-9, 9, 721))) >= 0)
Limit: the lower tail keeps its digits¶
Below nothing N(x) is the tail itself, never 1 less the upper, so a small chance is not lost to the rounding of 1: within three standard deviations its relative error is A&S's own, 5.1e-5, and binary32's; at five, A&S's 1.6e-3. Far out N is 0 and 1.
x = -np.linspace(1.0, 5.0, 401)
xs = np.float32(x).astype(float)
rel = np.abs(cdf(x) / ndtr(xs) - 1)
assert np.max(rel[xs >= -3]) <= 5.2e-5 + 16 * EPS, np.max(rel[xs >= -3])
assert np.max(rel) <= 1.6e-3 + 16 * EPS, np.max(rel)
assert list(cdf([-40.0, 40.0])) == [0.0, 1.0]
print(f'to 3 sd {np.max(rel[xs >= -3]):.2e}, to 5 sd {np.max(rel):.2e}')
to 3 sd 5.12e-05, to 5 sd 1.58e-03