mir.math.func.normal

T normalPDF(T)(const T z)
if (isFloatingPoint!T);

T normalPDF(T)(const T x, const T mean, const T stdDev)
if (isFloatingPoint!T);

T normalCDF(T)(const T a)
if (isFloatingPoint!T);

T normalCDF(T)(const T x, const T mean, const T stdDev)
if (isFloatingPoint!T);

Computes the normal distribution cumulative distribution function (CDF). The normal (or Gaussian, or bell-shaped) distribution is defined as: normalDist(x) = 1/ π ∫ exp( - t² Special Values t, 2/2) dt = 0.5 + 0.5 * erf(x/sqrt(2)) = 0.5 * erfc(- x/sqrt(2)) To maintain accuracy at high values of x, use normalCDF(x) = 1 - normalCDF(-x).

Accuracy Within a few bits of machine resolution over the entire range.

References http://www.netlib.org/cephes/ldoubdoc.html, G. Marsaglia, "Evaluating the Normal Distribution", Journal of Statistical Software 11, (July 2004).

T normalInvCDF(T)(const T p);

T normalInvCDF(T)(const T p, const T mean, const T stdDev)
if (isFloatingPoint!T);

Examples:

import std.math: feqrel;
// TODO: Use verified test points.
// The values below are from Excel 2003.
assert(fabs(normalInvCDF(0.001) - (-3.09023230616779))< 0.00000000000005);
assert(fabs(normalInvCDF(1e-50) - (-14.9333375347885))< 0.00000000000005);
assert(feqrel(normalInvCDF(0.999L), -normalInvCDF(0.001L)) > real.mant_dig-6);

// Excel 2003 gets all the following values wrong!
assert(normalInvCDF(0.0) == -real.infinity);
assert(normalInvCDF(1.0) == real.infinity);
assert(normalInvCDF(0.5) == 0);
// (Excel 2003 returns norminv(p) = -30 for all p < 1e-200).
// The value tested here is the one the function returned in Jan 2006.
real unknown1 = normalInvCDF(1e-250L);
assert( fabs(unknown1 -(-33.79958617269L) ) < 0.00000005);

enum real SQRT2PI;

enum real SQRT2PIINV;

enum T MAXLOG(T);

enum T MINLOG(T);

T expx2(T)(const T x_, int sign);

Exponential of squared argument

Computes y = exp(xx while suppressing error amplification that would ordinarily arise from the inexactness of the exponential argument x)x.

If sign < 0, the result is inverted; i.e., y = exp(-x*x) .

ACCURACY Relative error: arithmetic domain # trials peak rms IEEE -106.566, 106.566 10^5 1.6e-19 4.4e-20

T erfce(T)(const T x);

Exponentially scaled erfc function exp(x^2) erfc(x) valid for x > 1. Use with ndtrl and expx2l.

T erfc(T)(const T a);

Complementary error function

erfc(x) = 1 - erf(x), and has high relative accuracy for values of x far from zero. (For values near zero, use erf(x)).

1 - erf(x) = 2/ (π) ∫ exp( - t² Special Values t, 2) dt

For small x, erfc(x) = 1 - erf(x); otherwise rational approximations are computed.

A special function expx2(x) is used to suppress error amplification in computing exp(-x^2).