The normal density¶
The standard normal density, phi(x) = e^(-x^2 / 2) / sqrt(2 pi).
The 1 / sqrt(2 pi) is the law's own number, written to binary32's precision and beyond. Dimensionless in and out: a number of standard deviations, and a density per standard deviation. Source: Abramowitz & Stegun, Handbook of Mathematical Functions (1964), ch. 26 (Zelen & Severo). Checked against a reference in numpy, which has its own tests.
- Source
- Abramowitz & Stegun, Handbook of Mathematical Functions (1964), ch. 26 (Zelen & Severo).
- Proof
- reference, invariant, limit.
normal_pdf.trait
namespace lusine.finance;
float normal_pdf(float x)
{
return exp(-0.5 * x * x) * 0.3989422804014327;
}
Phi(x) = e^(-x^2 / 2) / sqrt(2 pi).
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_pdf.trait'))
// lusine.finance.normal_pdf: the standard normal density, phi(x) = e^(-x^2 / 2) / sqrt(2 pi). The 1 / sqrt(2 pi) is the
// law's own number, written to binary32's precision and beyond. Dimensionless in and out: a number of standard
// deviations, and a density per standard deviation.
// Source: Abramowitz & Stegun, Handbook of Mathematical Functions (1964), ch. 26 (Zelen & Severo).
// Checked against a reference in numpy, which has its own tests.
namespace lusine.finance;
float normal_pdf(float x)
{
return exp(-0.5 * x * x) * 0.3989422804014327;
}
The program¶
The kernel over a lattice of items, one an x. x is dimensionless, so it is stored as given; the independent account is scipy's density at the binary32 x the device holds.
print(trait_run.law_program('finance/normal_pdf.trait', 'normal_pdf'))
def pdf(x):
return trait_run.law('finance/normal_pdf.trait', 'normal_pdf', 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_pdf(items.x); }
Reference: scipy's density¶
Over [-9, 9] each value is within (4 + x^2 / 2) EPS of scipy's, relative: the exponent's argument is rounded once, and its rounding is relative to x^2 / 2, which exp carries into the result.
x = np.linspace(-9.0, 9.0, 721)
xs = np.float32(x).astype(float)
y = pdf(x)
rel = np.abs(y / norm.pdf(xs) - 1)
assert np.all(rel <= (4 + xs * xs / 2) * EPS), np.max(rel / EPS)
fig, ax = plt.subplots(figsize=(7, 3.5))
ax.plot(xs, rel / EPS)
ax.set_xlabel('x')
ax.set_ylabel('relative error (EPS)')
plt.show()
Invariant: even, and a density¶
phi(-x) is phi(x) bit for bit, and the density sums to one over a fine grid.
x = np.linspace(0.0, 10.0, 10001)
up, down = pdf(x), pdf(-x)
assert np.array_equal(up, down)
area = 2 * np.sum(up[1:]) * 1e-3 + up[0] * 1e-3
assert abs(area - 1) <= 1e-6, area
print('area', area)
area 1.0000000295804095
Limit: the peak and the tails¶
At nothing the density is 1 / sqrt(2 pi) as binary32 holds it; far out it is nothing.
peak, far = pdf([0.0, 40.0])
assert peak == np.float32(1 / np.sqrt(2 * np.pi)), peak
assert far == 0.0
print(peak, far)
0.3989422917366028 0.0