root/trunk/mathextra/normaldist.d

/**
 * Macros:
 *
 *  TABLE_SV = <table border=1 cellpadding=4 cellspacing=0>
 *      <caption>Special Values</caption>
 *      $0</table>
 *
 *  NAN = $(RED NAN)
 *  SUP = <span style="vertical-align:super;font-size:smaller">$0</span>
 */
module mathextra.normaldist;
private import std.math;
debug import std.stdio;

unittest {
    debug writefln("--- UnitTest: " __FILE__ " ---");
}

const real SQRT2PI =   0x1.40d931ff62705966p+1;    // 2.5066282746310005024
const real EXP_2  = 0.13533528323661269189L; /* exp(-2) */

/***
Computes the normal distribution function.

The normal (or Gaussian, or bell-shaped) distribution is
defined as:

normalDist(x) = &int;$(SUP x)<sub>-$(INFINITY)</sub>
   exp(-t$(SUP 2)/2)/$(SQRT)2*&pi;
    = 0.5 + 0.5 * erf(x/sqrt(2))
    = 0.5 * erfc(- x/sqrt(2))

To maintain accuracy at high values of x, use
normalDistribution(x) = 1 - normalDistribution(-x).

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

    $(TABLE_SV
    $(TR $(TH x)               $(TH cos(x)) $(TH invalid?)  )
    $(TR $(TD $(NAN))          $(TD $(NAN)) $(TD yes)   )
    $(TR $(TD &plusmn;&infin;) $(TD $(NAN)) $(TD yes)   )
    )
References:
G. Marsaglia, "Evaluating the Normal Distribution",
Journal of Statistical Software <b>11</b>, (July 2004).
*/
real normalDistribution(real x)
{
  return 0.5 * erfc(-x * SQRT1_2);
}

private {
static real P0[] = [
   -0x1.758f4d969484bfdcp-7,    // -0.011400139698853582732
   0x1.53cee17a59259dd2p-3, // 0.16592193750979583221
   -0x1.ea01e4400a9427a2p-1,    // -0.95704568177942689081
   0x1.61f7504a0105341ap+1, // 2.7653599130008302859
   -0x1.09475a594d0399f6p+2,    // -4.1449800369337538286
   0x1.7c59e7a0df99e3e2p+1, // 2.971493676711545292
   -0x1.87a81da52edcdf14p-1,    // -0.76495449677843806914
   0x1.1fb149fd3f83600cp-7  // 0.0087796794200550691607
];

static real Q0[] = [
   -0x1.64b92ae791e64bb2p-7,    // -0.010886331510064192632
   0x1.7585c7d597298286p-3, // 0.1823840725000038842
   -0x1.40011be4f7591ce6p+0,    // -1.2500169214248199725
   0x1.1fc067d8430a425ep+2, // 4.4961185085232139506
   -0x1.21008ffb1e7ccdf2p+3,    // -9.0313186554593813887
   0x1.3d1581cf9bc12fccp+3, // 9.9088753752567182205
   -0x1.53723a89fd8f083cp+2,    // -5.3038469646037218604
   0x1p+0   // 1
];

static real P1[] = [
   0x1.20ceea49ea142f12p-13,    // 0.00013771451113809605662
   0x1.cbe8a7267aea80bp-7,  // 0.014035302749980729871
   0x1.79fea765aa787c48p-2, // 0.36913549001712241224
   0x1.d1f59faa1f4c4864p+1, // 3.6403083401370131097
   0x1.1c22e426a013bb96p+4, // 17.75851836288460008
   0x1.a8675a0c51ef3202p+5, // 53.050464721918523919
   0x1.75782c4f83614164p+6, // 93.367356531518738722
   0x1.7a2f3d90948f1666p+6, // 94.546133288447683183
   0x1.5cd116ee4c088c3ap+5, // 43.602094518370966827
   0x1.1361e3eb6e3cc20ap+2  // 4.3028497504355521807
];

static real Q1[] = [
   0x1.3a4ce1406cea98fap-13,    // 0.00014987006762866754669
   0x1.f45332623335cda2p-7, // 0.015268706895221911913
   0x1.98f28bbd4b98db1p-2,  // 0.39936273901812389627
   0x1.ec3b24f9c698091cp+1, // 3.8455549449546995474
   0x1.1cc56ecda7cf58e4p+4, // 17.79820137342627204
   0x1.92c6f7376bf8c058p+5, // 50.347151215536627131
   0x1.4154c25aa47519b4p+6, // 80.332772651946720635
   0x1.1b321d3b927849eap+6, // 70.798939638914882544
   0x1.403a5f5a4ce7b202p+4, // 20.014251091705301368
   0x1p+0   // 1
];

static real P2[] = [
   0x1.8c124a850116a6d8p-21,    // 7.3774056430545041787e-07
   0x1.534abda3c2fb90bap-13,    // 0.0001617870121822776094
   0x1.29a055ec93a4718cp-7, // 0.0090828342009931074419
   0x1.6468e98aad6dd474p-3, // 0.17402822927913678347
   0x1.3dab2ef4c67a601cp+0, // 1.2408933017345389353
   0x1.e1fb3a1e70c67464p+1, // 3.7654793404231444828
   0x1.b6cce8035ff57b02p+2, // 6.8562564881284157607
   0x1.9f4c9e749ff35f62p+1  // 3.2445257253129069325
];

static real Q2[] = [
   0x1.af03f4fc0655e006p-21,    // 8.0282885006885383316e-07
   0x1.713192048d11fb2p-13, // 0.00017604524340842589303
   0x1.4357e5bbf5fef536p-7, // 0.0098676559208996361084
   0x1.7fdac8749985d43cp-3, // 0.18742901426157036096
   0x1.4a080c813a2d8e84p+0, // 1.2891853156563028786
   0x1.c3a4b423cdb41bdap+1, // 3.528463857156936774
   0x1.8160694e24b5557ap+2, // 6.0215094817275106307
   0x1p+0   // 1
];

static real P3[] = [
   -0x1.55da447ae3806168p-34,   // -7.7728283809481633868e-11
   -0x1.145635641f8778a6p-24,   // -6.4339663876133447143e-08
   -0x1.abf46d6b48040128p-17,   // -1.2754046756102807876e-05
   -0x1.7da550945da790fcp-11,   // -0.00072793152007373443093
   -0x1.aa0b2a31157775fap-8,    // -0.0065009096152460679857
   0x1.b11d97522eed26bcp-3, // 0.21148222178987070632
   0x1.1106d22f9ae89238p+1, // 2.1330206615874130532
   0x1.029a358e1e630f64p+1  // 2.0203310913027725356
];

static real Q3[] = [
   -0x1.74022dd5523e6f84p-34,   // -8.4584942637876803775e-11
   -0x1.2cb60d61e29ee836p-24,   // -7.0014768675591937804e-08
   -0x1.d19e6ec03a85e556p-17,   // -1.3876523894802171788e-05
   -0x1.9ea2a7b4422f6502p-11,   // -0.00079085420887378582886
   -0x1.c54b1e852f107162p-8,    // -0.0069167088997199649828
   0x1.e05268dd3c07989ep-3, // 0.23453218388704381964
   0x1.239c6aff14afbf82p+1, // 2.2782109971534491995
   0x1p+0   // 1
];

}

/******************************
 * Inverse of Normal distribution function
 *
 * Returns the argument, x, for which the area under the
 * Normal probability density function (integrated from
 * minus infinity to x) is equal to p.
 *
 * For small arguments 0 < p < exp(-2), the program computes
 * z = sqrt( -2 log(p) );  then the approximation is
 * x = z - log(z)/z  - (1/z) P(1/z) / Q(1/z) .
 * For larger arguments,  x/sqrt(2 pi) = w + w^3 R(w^2)/S(w^2)) ,
 * where w = p - 0.5 .
 */
real normalDistributionInv(real p)
in {
  assert(p>=0.0L && p<=1.0L); // domain error
}
body
{
    if (p == 0.0L) {
        return -real.infinity;
    }
    if( p == 1.0L ) {
        return real.infinity;
    }
    real x, z, y2, x0, x1;
    int code = 1;
    real y = p;
    if( y > (1.0L - EXP_2) ) {
        y = 1.0L - y;
        code = 0;
    }

    if ( y > EXP_2 ) {
        y = y - 0.5L;
        y2 = y * y;
        x = y + y * (y2 * poly( y2, P0)/poly( y2, Q0));
        x = x * SQRT2PI;
        return x;
    }

    x = sqrt( -2.0L * log(y) );
    x0 = x - log(x)/x;
    z = 1.0L/x;
    if( x < 8.0L ) {
        x1 = z * poly( z, P1)/poly( z, Q1);
    } else if( x < 32.0L ) {
        x1 = z * poly( z, P2)/poly( z, Q2);
    } else {
//  assert(0);
        x1 = z * poly( z, P3)/poly( z, Q3);
    }
    x = x0 - x1;
    if( code != 0 ) {
        x = -x;
    }
    return x;
}


unittest {
// The values below are from Excel 2003.
assert(fabs(normalDistributionInv(0.001) - (-3.09023230616779))< 0.00000000000005);
assert(fabs(normalDistributionInv(1e-50) - (-14.9333375347885))< 0.00000000000005);
assert(feqrel(normalDistributionInv(0.999), -normalDistributionInv(0.001))>real.mant_dig-6);

// Excel 2003 gets all the following values wrong!
assert(normalDistributionInv(0.0)==-real.infinity);
assert(normalDistributionInv(1.0)==real.infinity);
assert(normalDistributionInv(0.5)==0);

// I don't know the correct result for low values
// (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 = normalDistributionInv(1e-250L);
assert( fabs(unknown1 -(-33.79958617269L) ) < 0.00000005);
}

unittest {

//writefln("%.20g", normalDistribution(normalDistributionInv(1e-90)));
//normalDistribution(normalDistributionInv(1e-50))
assert(fabs(normalDistribution(1L) - (0.841344746068543))< 0.0000000000000005);
assert(isnan(normalDistribution(real.nan)));

//writefln("%.20f", normalDistributionInv(1e-250L));
//-3.09023230616779000000

//assert(normalDistribution(0)==0.5);
/+
  for (int i=0; i<=10; ++i) {
   writefln(i/10.0, "  ", normalDistributionInv(i/10.0));
  }
  real qq= 1-real.epsilon;
   writefln(qq, " %a  ", qq,  normalDistributionInv(qq));
   qq = real.min;
   writefln(qq, " %a  ", qq,  normalDistributionInv(qq));
   writefln(normalDistribution(9));
   writefln(normalDistribution(-150));
+/
  }