Skip to content

The normal distribution, live

The standard normal distribution, N(x), the chance a standard normal is below x, by Abramowitz and Stegun's 26.2.17: for z ≥ 0 the upper tail is φ(z) times a polynomial in t = 1 / (1 + p z), within 7.5 × 10⁻⁸ (lusine.finance.normal_cdf). Below zero, N(x) is that upper tail at −x itself, never one less the upper, so a small chance keeps its digits. Below, x is swept from −9 to 9 a tenth at a time on your GPU as you read, N at x and at −x; thirty-seven points a half apart are held to scipy's exact distribution, and seventeen in the lower tail from −1 to −5.

This page runs the program on your GPU through WebGPU, so it needs a browser with WebGPU (a recent Chrome, Edge, Safari or Firefox). Without it the program does not run and nothing is drawn below.

What it shows

  • The strip above the controls is green while every check holds, and turns red for good once one fails.
  • The distribution: the chart draws N(x) and 1 − N(−x), in percent, against x, a point a step, one over the other. Below it, x, N(x), the largest error of the thirty-seven points, absolute, in units of 10⁻⁸, and the largest in the lower tail, relative, in units of 10⁻⁵.
  • The checks: scipy's ndtr, N at each of the thirty-seven points within A&S's 7.5 × 10⁻⁸ and two binary32 epsilons (2⁻²³) of scipy's exact distribution at the x binary32 holds, the formula's error and the float's; N(−x) = 1 − N(x), rising, the two tails adding to one within an epsilon at every x swept, and N never falling as x rises; the lower tail, N's error relative to scipy's within A&S's own, 0.000052, to three standard deviations below nothing, and 0.0016 to five, each and 16 epsilons; far out, N(−40) nothing and N(40) one. Each is held at every step.

Every number here is computed by the program on the GPU; the page only reads the values back and draws them. The sweep's start and step and how far out is far are the program's settings, which the host may write; this page offers no control to change them yet.

The program

The points, the tail, the sweep and the checks, the panels the page draws, and the one procedure, main, run each step: x a tenth on, N at it, at −x, at the points and in the tail, and the checks.

// The normal distribution's live page: x swept from -9 to 9 a tenth at a time, the library's N at x and at -x, held to
// scipy's exact ndtr at thirty-seven points and in the lower tail at seventeen, with the page's checks computed here
// on the GPU and shown by the panels the program declares.

namespace normal_cdf;

lusine.version version = { { number = 1 } };

// The reader's settings: the sweep's first x and its step, and far out, 40.
host struct { float first = -9.0; float stride = 0.1; float far = 40.0; } settings;

// Thirty-seven x, -9 to 9 a half apart, each with N(x) by scipy's exact ndtr, as two binary32 numbers, its rounding and
// the rest; and the kernel's N(x) and its error, absolute.
struct point { float x; float reference; float low; float n; float error; };
host table<point> points(capacity = 37) = {
  { x = -9.0, reference = 1.1285884256611585e-19, low = -1.9707326600905217e-27 },
  { x = -8.5, reference = 9.479534487610775e-18, low = 3.3459246867680744e-25 },
  { x = -8.0, reference = 6.220960424517508e-16, low = 1.497542303579722e-23 },
  { x = -7.5, reference = 3.190891696919565e-14, low = -2.400867998999134e-22 },
  { x = -7.0, reference = 1.279812518982515e-12, low = 2.490331924137166e-20 },
  { x = -6.5, reference = 4.0160007264145037e-11, low = -1.425554175482303e-18 },
  { x = -6.0, reference = 9.865877004244794e-10, low = -5.5386786280521704e-17 },
  { x = -5.5, reference = 1.8989561922921894e-08, low = 5.429658010027805e-16 },
  { x = -5.0, reference = 2.8665158424701076e-07, low = -1.2367817748128e-14 },
  { x = -4.5, reference = 3.3976730264839716e-06, low = 9.82460843907372e-14 },
  { x = -4.0, reference = 3.16712430503685e-05, low = -1.2172486326192389e-12 },
  { x = -3.5, reference = 0.00023262908507604152, low = -6.0405166112109665e-12 },
  { x = -3.0, reference = 0.0013498980551958084, low = -2.356571535133689e-11 },
  { x = -2.5, reference = 0.006209665443748236, low = -1.1797210430763982e-10 },
  { x = -2.0, reference = 0.0227501317858696, low = 1.6230959631879927e-10 },
  { x = -1.5, reference = 0.06680720299482346, low = -1.725965370624749e-09 },
  { x = -1.0, reference = 0.15865525603294373, low = -2.1014865403401473e-09 },
  { x = -0.5, reference = 0.3085375428199768, low = -4.09399003586941e-09 }, { x = 0.0, reference = 0.5, low = 0.0 },
  { x = 0.5, reference = 0.6914624571800232, low = 4.09399003586941e-09 },
  { x = 1.0, reference = 0.8413447737693787, low = -2.7700835403265955e-08 },
  { x = 1.5, reference = 0.933192789554596, low = 9.176545745503972e-09 },
  { x = 2.0, reference = 0.9772498607635498, low = 7.288270875704939e-09 },
  { x = 2.5, reference = 0.993790328502655, low = 6.1715690335972795e-09 },
  { x = 3.0, reference = 0.998650074005127, low = 2.796324238829584e-08 },
  { x = 3.5, reference = 0.9997673630714417, low = 7.849522809522114e-09 },
  { x = 4.0, reference = 0.9999683499336243, low = -2.1175457831645872e-08 },
  { x = 4.5, reference = 0.9999966025352478, low = -2.083725414436799e-10 },
  { x = 5.0, reference = 0.9999997019767761, low = 1.137165206444024e-08 },
  { x = 5.5, reference = 1.0, low = -1.8989561922921894e-08 },
  { x = 6.0, reference = 1.0, low = -9.865877004244794e-10 },
  { x = 6.5, reference = 1.0, low = -4.0159986447463325e-11 },
  { x = 7.0, reference = 1.0, low = -1.2798651027878805e-12 },
  { x = 7.5, reference = 1.0, low = -3.186340080674199e-14 },
  { x = 8.0, reference = 1.0, low = -6.661338147750939e-16 }, { x = 8.5, reference = 1.0, low = 0.0 },
  { x = 9.0, reference = 1.0, low = 0.0 },
};

// The lower tail, x from -1 to -5 a quarter apart, each with N(x) by scipy's, and whether it is within three standard
// deviations; the kernel's N(x) and its error, relative.
struct tail { float x; float reference; float low; bool near; float n; float error; };
host table<tail> tails(capacity = 17) = {
  { x = -1.0, reference = 0.15865525603294373, low = -2.1014865403401473e-09, near = true },
  { x = -1.25, reference = 0.10564977675676346, low = -3.0899081071567025e-09, near = true },
  { x = -1.5, reference = 0.06680720299482346, low = -1.725965370624749e-09, near = true },
  { x = -1.75, reference = 0.040059156715869904, low = 1.4794718228294812e-10, near = true },
  { x = -2.0, reference = 0.0227501317858696, low = 1.6230959631879927e-10, near = true },
  { x = -2.25, reference = 0.012224473059177399, low = -4.0413269952743747e-10, near = true },
  { x = -2.5, reference = 0.006209665443748236, low = -1.1797210430763982e-10, near = true },
  { x = -2.75, reference = 0.0029797633178532124, low = -8.279865681970477e-11, near = true },
  { x = -3.0, reference = 0.0013498980551958084, low = -2.356571535133689e-11, near = true },
  { x = -3.25, reference = 0.0005770250572822988, low = -1.4891532451599687e-11, near = false },
  { x = -3.5, reference = 0.00023262908507604152, low = -6.0405166112109665e-12, near = false },
  { x = -3.75, reference = 8.84172841324471e-05, low = 1.068356666154846e-12, near = false },
  { x = -4.0, reference = 3.16712430503685e-05, low = -1.2172486326192389e-12, near = false },
  { x = -4.25, reference = 1.0688525435398333e-05, low = 3.3953607070022485e-13, near = false },
  { x = -4.5, reference = 3.3976730264839716e-06, low = 9.82460843907372e-14, near = false },
  { x = -4.75, reference = 1.017083263832319e-06, low = -2.1263615258208628e-14, near = false },
  { x = -5.0, reference = 2.8665158424701076e-07, low = -1.2367817748128e-14, near = false },
};

// The sweep: the step, from 0 to 180 and round again; x, N(x) and N(-x), and N at the step before; N(x) and 1 - N(-x)
// in percent, for the chart; the largest absolute error of the thirty-seven, in units of 1e-8, and the largest
// relative error in the tail, in units of 1e-5.
struct {
  integer step = -1; float x; float n; float mirrored; float before; float percent; float complement; float worst;
  float tail;
} sweep;

// The checks: every N(x) within A&S's 7.5e-8 and two binary32 epsilons (2^-23) of scipy's, the formula's error and the
// float's; N(-x) + N(x) one within an epsilon, and N never falling as x rises; the lower tail its own, never 1 less the
// upper, so its relative error within three standard deviations A&S's 5.2e-5 and to five 1.6e-3, each and 16 epsilons;
// and far out N 0 and 1.
struct { bool reference = true; bool symmetric = true; bool tail = true; bool limits = true; } checks;

// A strip of lamps the page draws, lit green while every check holds and red once one fails.
struct lamp { vector<float, 4> colour; };
lattice<lamp, 2> lamps(spacing = 1 m, boundary = wall);

vector<float, 4> lamp_colour(bool holds)
{
  return holds ? vector<float, 4>(0.15, 0.6, 0.3, 1.0) : vector<float, 4>(0.75, 0.15, 0.15, 1.0);
}

float eps() { return 1.0 / 8388608.0; }

float absolute(float n, float reference, float low) { return abs(n - reference - low); }
float relative(float n, float reference, float low) { return abs(n - reference - low) / reference; }

integer outside(table<float> errors) { return count(errors) where errors > 0.000000075 + 2.0 * eps(); }

integer outside_tail(table<float> errors, table<bool> nears)
{
  return count(errors) where errors > (nears ? 0.000052 : 0.0016) + 16.0 * eps();
}

output table<lusine.widget> panels = {
  lusine.series { { panel = "the distribution", chart = "N(x) and 1 - N(-x), %", quantity = sweep.percent,
                    steps = 181, against = sweep.x, label = "N(x)" } },
  lusine.series { { panel = "the distribution", chart = "N(x) and 1 - N(-x), %", quantity = sweep.complement,
                    steps = 181, against = sweep.x, label = "1 - N(-x)" } },
  lusine.value { { panel = "the distribution", quantity = sweep.x, label = "x" } },
  lusine.value { { panel = "the distribution", quantity = sweep.percent, label = "N(x), %" } },
  lusine.value { { panel = "the distribution", quantity = sweep.worst, label = "largest error, 1e-8" } },
  lusine.value { { panel = "the distribution", quantity = sweep.tail, label = "largest in the tail, 1e-5" } },
  lusine.value { { panel = "the checks", quantity = checks.reference, label = "scipy's ndtr" } },
  lusine.value { { panel = "the checks", quantity = checks.symmetric, label = "N(-x) = 1 - N(x), rising" } },
  lusine.value { { panel = "the checks", quantity = checks.tail, label = "the lower tail" } },
  lusine.value { { panel = "the checks", quantity = checks.limits, label = "far out" } },
};

void main()
{
  sweep.step = sweep.step >= 180 ? 0 : sweep.step + 1;
  sweep.x = settings.first + float(sweep.step) * settings.stride;
  sweep.before = sweep.step == 0 ? 0.0 : sweep.n;
  sweep.n = lusine.finance.normal_cdf(sweep.x);
  sweep.mirrored = lusine.finance.normal_cdf(x = -sweep.x);
  checks.symmetric = checks.symmetric && abs(sweep.n + sweep.mirrored - 1.0) <= eps() && sweep.n >= sweep.before;
  sweep.percent = 100.0 * sweep.n;
  sweep.complement = 100.0 * (1.0 - sweep.mirrored);
  points.n = lusine.finance.normal_cdf(points.x);
  points.error = absolute(points.n, points.reference, points.low);
  tails.n = lusine.finance.normal_cdf(tails.x);
  tails.error = relative(tails.n, tails.reference, tails.low);
  sweep.worst = 100000000.0 * max(points.error);
  sweep.tail = 100000.0 * max(tails.error);
  checks.reference = checks.reference && outside(errors = points.error) == 0;
  checks.tail = checks.tail && outside_tail(errors = tails.error, nears = tails.near) == 0;
  checks.limits = checks.limits && lusine.finance.normal_cdf(x = -settings.far) == 0.0
                  && lusine.finance.normal_cdf(x = settings.far) == 1.0;
  lamps.colour = lamp_colour(holds = checks.reference && checks.symmetric && checks.tail && checks.limits);
}

The kernel's other proof (721 points against scipy's, and 401 in the lower tail) is on its proof page.