Skip to content

The normal density, live

The standard normal density, φ(x) = e^(−x²/2) / √(2π), dimensionless in and out: a number of standard deviations, and a density per standard deviation (lusine.finance.normal_pdf). Below, x is swept from −9 to 9 a tenth at a time on your GPU as you read, the library's density at each, at −x beside it; thirty-seven points a half apart are held to a reference in float64; and the density is taken over a grid a quarter apart from −10 to 10 and summed.

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 density: the chart draws φ(x) in thousandths against x, a point a step, the bell filling in from the left. Below it, x, the density there, the largest error of the thirty-seven points over its bound, and the grid's sum.
  • The checks: the float64 reference, the density at each of the thirty-seven points within (4 + x²/2) binary32 epsilons (2⁻²³) of e^(−x²/2) / √(2π) in float64, relative: the exponent's argument is rounded once, and its rounding, relative to x²/2, exp carries into the result; even, φ(−x) the same bits as φ(x) at every x swept; a density, h Σ φ(k h) over the grid one within 10⁻⁶, the trapezoid rule exact for the bell to far below binary32's rounding at h = ¼; the peak and the tails, φ(0) 1/√(2π) as binary32 holds it, bit for bit, and φ(40) nothing. 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, the grid's spacing and how far out the tail is taken are the program's settings, which the host may write; this page offers no control to change them yet.

The program

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

// The normal density's live page: x swept from -9 to 9 a tenth at a time, the library's density at each beside a
// reference in float64 at thirty-seven points, the density summed over a grid a quarter apart, its evenness, its peak
// and its tails, with the page's checks computed here on the GPU and shown by the panels the program declares.

namespace normal_pdf;

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

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

// Thirty-seven x, -9 to 9 a half apart, each with the density there in float64, e^(-x^2 / 2) / sqrt(2 pi), as two
// binary32 numbers, its rounding and the rest; and the kernel's error relative to it, in binary32 epsilons (2^-23),
// against the bound, 4 + x^2 / 2.
struct point { float x; float reference; float low; float density; float error; float bound; };
host table<point> points(capacity = 37) = {
  { x = -9.0, reference = 1.0279773620817797e-18, low = -4.914888091972684e-27 },
  { x = -8.5, reference = 8.166235380055849e-17, low = 2.516136938227098e-24 },
  { x = -8.0, reference = 5.05227116246636e-15, low = -7.892946525011385e-23 },
  { x = -7.5, reference = 2.43432064870458e-13, low = -1.1567557026644364e-20 },
  { x = -7.0, reference = 9.134720757586479e-12, low = -3.4922187689281045e-19 },
  { x = -6.5, reference = 2.669556520284999e-10, low = 9.447785641339765e-18 },
  { x = -6.0, reference = 6.075882907907726e-09, low = -5.808443732501639e-17 },
  { x = -5.5, reference = 1.0769760194762057e-07, low = -1.522187768796771e-15 },
  { x = -5.0, reference = 1.486719497734157e-06, low = 1.7000140986773993e-14 },
  { x = -4.5, reference = 1.598374183231499e-05, low = -7.254095327896559e-13 },
  { x = -4.0, reference = 0.00013383022451307625, low = 1.2518090947502558e-12 },
  { x = -3.5, reference = 0.0008726827218197286, low = -2.677396845895963e-11 },
  { x = -3.0, reference = 0.004431848414242268, low = -2.3042600813538705e-12 },
  { x = -2.5, reference = 0.01752830110490322, low = -6.113347050806794e-10 },
  { x = -2.0, reference = 0.05399096757173538, low = -1.0585473608770712e-09 },
  { x = -1.5, reference = 0.1295175999403, low = -4.274408382798356e-09 },
  { x = -1.0, reference = 0.2419707179069519, low = 6.6121916830752525e-09 },
  { x = -0.5, reference = 0.3520653247833252, low = 1.9809742735077407e-09 },
  { x = 0.0, reference = 0.3989422917366028, low = -1.1335170135851058e-08 },
  { x = 0.5, reference = 0.3520653247833252, low = 1.9809742735077407e-09 },
  { x = 1.0, reference = 0.2419707179069519, low = 6.6121916830752525e-09 },
  { x = 1.5, reference = 0.1295175999403, low = -4.274408382798356e-09 },
  { x = 2.0, reference = 0.05399096757173538, low = -1.0585473608770712e-09 },
  { x = 2.5, reference = 0.01752830110490322, low = -6.113347050806794e-10 },
  { x = 3.0, reference = 0.004431848414242268, low = -2.3042600813538705e-12 },
  { x = 3.5, reference = 0.0008726827218197286, low = -2.677396845895963e-11 },
  { x = 4.0, reference = 0.00013383022451307625, low = 1.2518090947502558e-12 },
  { x = 4.5, reference = 1.598374183231499e-05, low = -7.254095327896559e-13 },
  { x = 5.0, reference = 1.486719497734157e-06, low = 1.7000140986773993e-14 },
  { x = 5.5, reference = 1.0769760194762057e-07, low = -1.522187768796771e-15 },
  { x = 6.0, reference = 6.075882907907726e-09, low = -5.808443732501639e-17 },
  { x = 6.5, reference = 2.669556520284999e-10, low = 9.447785641339765e-18 },
  { x = 7.0, reference = 9.134720757586479e-12, low = -3.4922187689281045e-19 },
  { x = 7.5, reference = 2.43432064870458e-13, low = -1.1567557026644364e-20 },
  { x = 8.0, reference = 5.05227116246636e-15, low = -7.892946525011385e-23 },
  { x = 8.5, reference = 8.166235380055849e-17, low = 2.516136938227098e-24 },
  { x = 9.0, reference = 1.0279773620817797e-18, low = -4.914888091972684e-27 },
};

// The grid, x = k h for k from -40 to 40: each one's density.
struct node { integer k; float density; };
host table<node> grid(capacity = 81) = {
  { k = -40 }, { k = -39 }, { k = -38 }, { k = -37 }, { k = -36 }, { k = -35 }, { k = -34 }, { k = -33 }, { k = -32 },
  { k = -31 }, { k = -30 }, { k = -29 }, { k = -28 }, { k = -27 }, { k = -26 }, { k = -25 }, { k = -24 }, { k = -23 },
  { k = -22 }, { k = -21 }, { k = -20 }, { k = -19 }, { k = -18 }, { k = -17 }, { k = -16 }, { k = -15 }, { k = -14 },
  { k = -13 }, { k = -12 }, { k = -11 }, { k = -10 }, { k = -9 }, { k = -8 }, { k = -7 }, { k = -6 }, { k = -5 },
  { k = -4 }, { k = -3 }, { k = -2 }, { k = -1 }, { k = 0 }, { k = 1 }, { k = 2 }, { k = 3 }, { k = 4 }, { k = 5 },
  { k = 6 }, { k = 7 }, { k = 8 }, { k = 9 }, { k = 10 }, { k = 11 }, { k = 12 }, { k = 13 }, { k = 14 }, { k = 15 },
  { k = 16 }, { k = 17 }, { k = 18 }, { k = 19 }, { k = 20 }, { k = 21 }, { k = 22 }, { k = 23 }, { k = 24 },
  { k = 25 }, { k = 26 }, { k = 27 }, { k = 28 }, { k = 29 }, { k = 30 }, { k = 31 }, { k = 32 }, { k = 33 },
  { k = 34 }, { k = 35 }, { k = 36 }, { k = 37 }, { k = 38 }, { k = 39 }, { k = 40 },
};

// The sweep: the step, from 0 to 180 and round again; x, the density there and at -x, and in thousandths for the chart;
// the largest error of the points, in epsilons over its bound; and the grid's sum.
struct {
  integer step = -1; float x; float density; float mirrored; float thousandths; float worst; float area;
} sweep;

// The checks: every point's density within (4 + x^2 / 2) epsilons of the float64 reference, relative, the exponent's
// argument rounded once and exp carrying its rounding, relative to x^2 / 2, into the result; the density even, phi(-x)
// phi(x) bit for bit; the grid's h sum phi(k h) one within 1e-6, the trapezoid rule exact for the bell to far below
// binary32; and the peak 1 / sqrt(2 pi) as binary32 holds it, bit for bit, and 40 out nothing.
struct { bool reference = true; bool even = true; bool area = 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 error_of(float density, float reference, float low) { return abs(density - reference - low) / reference / eps(); }
float bound_of(float x) { return 4.0 + x * x / 2.0; }

integer outside(table<float> errors, table<float> bounds) { return count(errors) where errors > bounds; }

float worst(table<float> errors, table<float> bounds) { return max(errors / bounds); }

float at(integer k, float spacing) { return float(k) * spacing; }

output table<lusine.widget> panels = {
  lusine.series { { panel = "the density", chart = "phi(x), thousandths", quantity = sweep.thousandths, steps = 181,
                    against = sweep.x, label = "the kernel" } },
  lusine.value { { panel = "the density", quantity = sweep.x, label = "x" } },
  lusine.value { { panel = "the density", quantity = sweep.density, label = "phi(x)" } },
  lusine.value { { panel = "the density", quantity = sweep.worst, label = "largest error over its bound" } },
  lusine.value { { panel = "the density", quantity = sweep.area, label = "h sum phi(k h)" } },
  lusine.value { { panel = "the checks", quantity = checks.reference, label = "the float64 reference" } },
  lusine.value { { panel = "the checks", quantity = checks.even, label = "even" } },
  lusine.value { { panel = "the checks", quantity = checks.area, label = "a density" } },
  lusine.value { { panel = "the checks", quantity = checks.limits, label = "the peak and the tails" } },
};

void main()
{
  sweep.step = sweep.step >= 180 ? 0 : sweep.step + 1;
  sweep.x = settings.first + float(sweep.step) * settings.stride;
  sweep.density = lusine.finance.normal_pdf(sweep.x);
  sweep.mirrored = lusine.finance.normal_pdf(x = -sweep.x);
  sweep.thousandths = 1000.0 * sweep.density;
  points.density = lusine.finance.normal_pdf(points.x);
  points.error = error_of(points.density, points.reference, points.low);
  points.bound = bound_of(points.x);
  sweep.worst = worst(errors = points.error, bounds = points.bound);
  grid.density = lusine.finance.normal_pdf(x = at(grid.k, settings.spacing));
  sweep.area = settings.spacing * sum(grid.density);
  checks.reference = checks.reference && outside(errors = points.error, bounds = points.bound) == 0;
  checks.even = checks.even && sweep.mirrored == sweep.density;
  checks.area = checks.area && abs(sweep.area - 1.0) <= 0.000001;
  checks.limits = checks.limits && lusine.finance.normal_pdf(x = 0.0 * settings.far) == 0.3989422917366028
                  && lusine.finance.normal_pdf(x = settings.far) == 0.0;
  lamps.colour = lamp_colour(holds = checks.reference && checks.even && checks.area && checks.limits);
}

The kernel's other proof (721 points against scipy's density, and the area over a finer grid) is on its proof page.