mir.math.func.normal source code

1 /**
2 License: $(HTTP www.apache.org/licenses/LICENSE-2.0, Apache-2.0)
3 */
4 /**
5  * Error Functions and Normal Distribution.
6  *
7  * Copyright: Based on the CEPHES math library, which is
8  *            Copyright (C) 1994 Stephen L. Moshier (moshier@world.std.com).
9  * Authors:   Stephen L. Moshier, ported to D by Don Clugston and David Nadlinger. Adopted to Mir by Ilia Ki.
10  */
11 /**
12  * Macros:
13  *  NAN = $(RED NAN)
14  *  INTEGRAL = &#8747;
15  *  POWER = $1<sup>$2</sup>
16  *      <caption>Special Values</caption>
17  *      $0</table>
18  */
19 
20 module mir.math.func.normal;
21 
22 import std.traits: isFloatingPoint;
23 import mir.math.common;
24 
25 @safe pure nothrow @nogc:
26 
27 ///
28 T normalPDF(T)(const T z)
29  if (isFloatingPoint!T)
30 {
31     import mir.math.common: sqrt, exp;
32     return T(SQRT2PIINV) * exp(z * z * -0.5);
33 }
34 
35 /// ditto
36 T normalPDF(T)(const T x, const T mean, const T stdDev)
37 if (isFloatingPoint!T)
38 {
39     return normalPDF((x - mean) / stdDev) / stdDev;
40 }
41 
42 /++
43 Computes the normal distribution cumulative distribution function (CDF).
44 The normal (or Gaussian, or bell-shaped) distribution is
45 defined as:
46 normalDist(x) = 1/$(SQRT) &pi; $(INTEGRAL -$(INFINITY), x) exp( - $(POWER t, 2)/2) dt
47     = 0.5 + 0.5 * erf(x/sqrt(2))
48     = 0.5 * erfc(- x/sqrt(2))
49 To maintain accuracy at high values of x, use
50 normalCDF(x) = 1 - normalCDF(-x).
51 Accuracy:
52 Within a few bits of machine resolution over the entire
53 range.
54 References:
55 $(LINK http://www.netlib.org/cephes/ldoubdoc.html),
56 G. Marsaglia, "Evaluating the Normal Distribution",
57 Journal of Statistical Software <b>11</b>, (July 2004).
58 +/
59 T normalCDF(T)(const T a)
60     if (isFloatingPoint!T)
61 {
62     pragma(inline, false);
63     import mir.math.constant: SQRT1_2;
64 
65     T x = a * T(SQRT1_2);
66     T z = fabs(x);
67 
68     if (z < 1)
69     {
70         return 0.5f + 0.5f * erf(x);
71     }
72     else
73     {
74         T y = 0.5f * erfce(z);
75         /* Multiply by exp(-x^2 / 2)  */
76         z = expx2(a, -1);
77         y = y * sqrt(z);
78         if (y != y)
79             y = 0;
80         if (x > 0)
81             y = 1 - y;
82         return y;
83     }
84 }
85 
86 /// ditto
87 T normalCDF(T)(const T x, const T mean, const T stdDev)
88     if (isFloatingPoint!T)
89 {
90     return normalCDF((x - mean) / stdDev);
91 }
92 
93 version(mir_test)
94 @safe
95 unittest
96 {
97     assert(fabs(normalCDF(1.0L) - (0.841344746068543))< 0.0000000000000005);
98     assert(normalCDF(double.infinity) == 1);
99     assert(normalCDF(100.0) == 1);
100     assert(normalCDF(8.0) < 1);
101     assert(normalCDF(-double.infinity) == 0);
102     assert(normalCDF(-100.0) == 0);
103     assert(normalCDF(-8.0) > 0);
104 }
105 
106 ///
107 T normalInvCDF(T)(const T p)
108 in {
109   assert(p >= 0 && p <= 1, "Domain error");
110 }
111 do
112 {
113     pragma(inline, false);
114     if (p <= 0 || p >= 1)
115     {
116         if (p == 0)
117             return -T.infinity;
118         if ( p == 1 )
119             return T.infinity;
120         return T.nan; // domain error
121     }
122     int code = 1;
123     T y = p;
124     if ( y > (1 - T(EXP_2)) )
125     {
126         y = 1 - y;
127         code = 0;
128     }
129 
130     T x, z, y2, x0, x1;
131 
132     if ( y > T(EXP_2) )
133     {
134         y = y - 0.5L;
135         y2 = y * y;
136         x = y + y * (y2 * rationalPoly!(P0, Q0)(y2));
137         return x * double(SQRT2PI);
138     }
139 
140     x = sqrt( -2 * log(y) );
141     x0 = x - log(x)/x;
142     z = 1/x;
143     if ( x < 8 )
144     {
145         x1 = z * rationalPoly!(P1, Q1)(z);
146     }
147     else if ( x < 32 )
148     {
149         x1 = z * rationalPoly!(P2, Q2)(z);
150     }
151     else
152     {
153         x1 = z * rationalPoly!(P3, Q3)(z);
154     }
155     x = x0 - x1;
156     if ( code != 0 )
157     {
158         x = -x;
159     }
160     return x;
161 }
162 
163 /// ditto
164 T normalInvCDF(T)(const T p, const T mean, const T stdDev)
165     if (isFloatingPoint!T)
166 {
167     return normalInvCDF(p) * stdDev + mean;
168 }
169 
170 ///
171 version(mir_test)
172 @safe unittest
173 {
174     import std.math: feqrel;
175     // TODO: Use verified test points.
176     // The values below are from Excel 2003.
177     assert(fabs(normalInvCDF(0.001) - (-3.09023230616779))< 0.00000000000005);
178     assert(fabs(normalInvCDF(1e-50) - (-14.9333375347885))< 0.00000000000005);
179     assert(feqrel(normalInvCDF(0.999L), -normalInvCDF(0.001L)) > real.mant_dig-6);
180 
181     // Excel 2003 gets all the following values wrong!
182     assert(normalInvCDF(0.0) == -real.infinity);
183     assert(normalInvCDF(1.0) == real.infinity);
184     assert(normalInvCDF(0.5) == 0);
185     // (Excel 2003 returns norminv(p) = -30 for all p < 1e-200).
186     // The value tested here is the one the function returned in Jan 2006.
187     real unknown1 = normalInvCDF(1e-250L);
188     assert( fabs(unknown1 -(-33.79958617269L) ) < 0.00000005);
189 }
190 
191 ///
192 enum real SQRT2PI = 2.50662827463100050241576528481104525L; // sqrt(2pi)
193 ///
194 enum real SQRT2PIINV = 1 / SQRT2PI; // 1 / sqrt(2pi)
195 
196 package(mir) {
197 
198 enum real EXP_2  = 0.135335283236612691893999494972484403L; /* exp(-2) */
199 
200 ///
201 enum T MAXLOG(T) = log(T.max);
202 ///
203 enum T MINLOG(T) = log(T.min_normal * T.epsilon); // log(smallest denormal);
204 }
205 
206 /**
207  *  Exponential of squared argument
208  *
209  * Computes y = exp(x*x) while suppressing error amplification
210  * that would ordinarily arise from the inexactness of the
211  * exponential argument x*x.
212  *
213  * If sign < 0, the result is inverted; i.e., y = exp(-x*x) .
214  *
215  * ACCURACY:
216  *                      Relative error:
217  * arithmetic      domain        # trials      peak         rms
218  *   IEEE     -106.566, 106.566    10^5       1.6e-19     4.4e-20
219  */
220 package(mir)
221 T expx2(T)(const T x_, int sign)
222 {
223     /*
224     Cephes Math Library Release 2.9:  June, 2000
225     Copyright 2000 by Stephen L. Moshier
226     */
227     const T M = 32_768.0;
228     const T MINV = 3.0517578125e-5L;
229 
230     T x = fabs(x_);
231     if (sign < 0)
232         x = -x;
233 
234   /* Represent x as an exact multiple of M plus a residual.
235      M is a power of 2 chosen so that exp(m * m) does not overflow
236      or underflow and so that |x - m| is small.  */
237     T m = MINV * floor(M * x + 0.5f);
238     T f = x - m;
239 
240     /* x^2 = m^2 + 2mf + f^2 */
241     T u = m * m;
242     T u1 = 2 * m * f  +  f * f;
243 
244     if (sign < 0)
245     {
246         u = -u;
247         u1 = -u1;
248     }
249 
250     if (u + u1 > MAXLOG!T)
251         return T.infinity;
252 
253     /* u is exact, u1 is small.  */
254     return exp(u) * exp(u1);
255 }
256 
257 /**
258    Exponentially scaled erfc function
259    exp(x^2) erfc(x)
260    valid for x > 1.
261    Use with ndtrl and expx2l.  */
262 package(mir)
263 T erfce(T)(const T x)
264 {
265     T y = 1 / x;
266 
267     if (x < 8)
268     {
269         return rationalPoly!(P, Q)(y);
270     }
271     else
272     {
273         return y * rationalPoly!(R, S)(y * y);
274     }
275 }
276 
277 private T rationalPoly(alias numerator, alias denominator, T)(const T x) pure nothrow
278 {
279     return x.poly!numerator / x.poly!denominator;
280 }
281 
282 private T poly(alias vec, T)(const T x)
283 {
284     import mir.internal.utility: Iota;
285     T y = T(vec[$ - 1]);
286     foreach_reverse(i; Iota!(vec.length - 1))
287     {
288         y *= x;
289         y += T(vec[i]);
290     }
291     return y;
292 }
293 
294 private {
295 
296 /* erfc(x) = exp(-x^2) P(1/x)/Q(1/x)
297    1/8 <= 1/x <= 1
298    Peak relative error 5.8e-21  */
299 immutable real[10] P = [ -0x1.30dfa809b3cc6676p-17, 0x1.38637cd0913c0288p+18,
300    0x1.2f015e047b4476bp+22, 0x1.24726f46aa9ab08p+25, 0x1.64b13c6395dc9c26p+27,
301    0x1.294c93046ad55b5p+29, 0x1.5962a82f92576dap+30, 0x1.11a709299faba04ap+31,
302    0x1.11028065b087be46p+31, 0x1.0d8ef40735b097ep+30
303 ];
304 
305 immutable real[11] Q = [ 0x1.14d8e2a72dec49f4p+19, 0x1.0c880ff467626e1p+23,
306    0x1.04417ef060b58996p+26, 0x1.404e61ba86df4ebap+28, 0x1.0f81887bc82b873ap+30,
307    0x1.4552a5e39fb49322p+31, 0x1.11779a0ceb2a01cep+32, 0x1.3544dd691b5b1d5cp+32,
308    0x1.a91781f12251f02ep+31, 0x1.0d8ef3da605a1c86p+30, 1.0
309 ];
310 
311 /* erfc(x) = exp(-x^2) 1/x R(1/x^2) / S(1/x^2)
312     1/128 <= 1/x < 1/8
313     Peak relative error 1.9e-21  */
314 immutable real[5] R = [ 0x1.b9f6d8b78e22459ep-6, 0x1.1b84686b0a4ea43ap-1,
315     0x1.b8f6aebe96000c2ap+1, 0x1.cb1dbedac27c8ec2p+2, 0x1.cf885f8f572a4c14p+1
316 ];
317 
318 immutable real[6] S = [
319     0x1.87ae3cae5f65eb5ep-5, 0x1.01616f266f306d08p+0, 0x1.a4abe0411eed6c22p+2,
320     0x1.eac9ce3da600abaap+3, 0x1.5752a9ac2faebbccp+3, 1.0
321 ];
322 
323 /* erf(x)  = x P(x^2)/Q(x^2)
324     0 <= x <= 1
325     Peak relative error 7.6e-23  */
326 immutable real[7] T = [ 0x1.0da01654d757888cp+20, 0x1.2eb7497bc8b4f4acp+17,
327     0x1.79078c19530f72a8p+15, 0x1.4eaf2126c0b2c23p+11, 0x1.1f2ea81c9d272a2ep+8,
328     0x1.59ca6e2d866e625p+2, 0x1.c188e0b67435faf4p-4
329 ];
330 
331 immutable real[7] U = [ 0x1.dde6025c395ae34ep+19, 0x1.c4bc8b6235df35aap+18,
332     0x1.8465900e88b6903ap+16, 0x1.855877093959ffdp+13, 0x1.e5c44395625ee358p+9,
333     0x1.6a0fed103f1c68a6p+5, 1.0
334 ];
335 
336 }
337 
338 package(mir)
339 F erf(F)(const F x)
340  if(isFloatingPoint!F)
341 {
342     if (x == 0)
343         return x; // deal with negative zero
344     if (x == -F.infinity)
345         return -1;
346     if (x == F.infinity)
347         return 1;
348     immutable ax = fabs(x);
349     if (ax > 1)
350         return 1 - erfc(x);
351 
352     return x * rationalPoly!(T, U)(x * x);
353 }
354 
355 /**
356  *  Complementary error function
357  *
358  * erfc(x) = 1 - erf(x), and has high relative accuracy for
359  * values of x far from zero. (For values near zero, use erf(x)).
360  *
361  *  1 - erf(x) =  2/ $(SQRT)(&pi;)
362  *     $(INTEGRAL x, $(INFINITY)) exp( - $(POWER t, 2)) dt
363  *
364  *
365  * For small x, erfc(x) = 1 - erf(x); otherwise rational
366  * approximations are computed.
367  *
368  * A special function expx2(x) is used to suppress error amplification
369  * in computing exp(-x^2).
370  */
371 package(mir)
372 T erfc(T)(const T a)
373 {
374     if (a == T.infinity)
375         return 0;
376     if (a == -T.infinity)
377         return 2;
378 
379     immutable x = (a < 0) ? -a : a;
380 
381     if (x < 1)
382         return 1 - erf(a);
383 
384     if (-a * a < -MAXLOG!T)
385     {
386         // underflow
387         if (a < 0) return 2;
388         else return 0;
389     }
390 
391     T y;
392     immutable z = expx2(a, -1);
393 
394     y = 1 / x;
395     if (x < 8)
396         y = z * rationalPoly!(P, Q)(y);
397     else
398         y = z * y * rationalPoly!(R, S)(y * y);
399 
400     if (a < 0)
401         y = 2 - y;
402 
403     if (y == 0)
404     {
405         // underflow
406         if (a < 0) return 2;
407         else return 0;
408     }
409 
410     return y;
411 }
412 
413 private:
414 
415 static immutable real[8] P0 =
416 [ -0x1.758f4d969484bfdcp-7, 0x1.53cee17a59259dd2p-3,
417 -0x1.ea01e4400a9427a2p-1,  0x1.61f7504a0105341ap+1, -0x1.09475a594d0399f6p+2,
418 0x1.7c59e7a0df99e3e2p+1, -0x1.87a81da52edcdf14p-1,  0x1.1fb149fd3f83600cp-7
419 ];
420 
421 static immutable real[8] Q0 =
422 [ -0x1.64b92ae791e64bb2p-7, 0x1.7585c7d597298286p-3,
423 -0x1.40011be4f7591ce6p+0, 0x1.1fc067d8430a425ep+2, -0x1.21008ffb1e7ccdf2p+3,
424 0x1.3d1581cf9bc12fccp+3, -0x1.53723a89fd8f083cp+2, 1.0
425 ];
426 
427 static immutable real[10] P1 =
428 [ 0x1.20ceea49ea142f12p-13, 0x1.cbe8a7267aea80bp-7,
429 0x1.79fea765aa787c48p-2, 0x1.d1f59faa1f4c4864p+1, 0x1.1c22e426a013bb96p+4,
430 0x1.a8675a0c51ef3202p+5, 0x1.75782c4f83614164p+6, 0x1.7a2f3d90948f1666p+6,
431 0x1.5cd116ee4c088c3ap+5, 0x1.1361e3eb6e3cc20ap+2
432 ];
433 
434 static immutable real[10] Q1 =
435 [ 0x1.3a4ce1406cea98fap-13, 0x1.f45332623335cda2p-7,
436 0x1.98f28bbd4b98db1p-2, 0x1.ec3b24f9c698091cp+1, 0x1.1cc56ecda7cf58e4p+4,
437 0x1.92c6f7376bf8c058p+5, 0x1.4154c25aa47519b4p+6, 0x1.1b321d3b927849eap+6,
438 0x1.403a5f5a4ce7b202p+4, 1.0
439 ];
440 
441 static immutable real[8] P2 =
442 [ 0x1.8c124a850116a6d8p-21, 0x1.534abda3c2fb90bap-13,
443 0x1.29a055ec93a4718cp-7, 0x1.6468e98aad6dd474p-3, 0x1.3dab2ef4c67a601cp+0,
444 0x1.e1fb3a1e70c67464p+1, 0x1.b6cce8035ff57b02p+2, 0x1.9f4c9e749ff35f62p+1
445 ];
446 
447 static immutable real[8] Q2 =
448 [ 0x1.af03f4fc0655e006p-21, 0x1.713192048d11fb2p-13,
449 0x1.4357e5bbf5fef536p-7, 0x1.7fdac8749985d43cp-3, 0x1.4a080c813a2d8e84p+0,
450 0x1.c3a4b423cdb41bdap+1, 0x1.8160694e24b5557ap+2, 1.0
451 ];
452 
453 static immutable real[8] P3 =
454 [ -0x1.55da447ae3806168p-34, -0x1.145635641f8778a6p-24,
455 -0x1.abf46d6b48040128p-17, -0x1.7da550945da790fcp-11, -0x1.aa0b2a31157775fap-8,
456 0x1.b11d97522eed26bcp-3, 0x1.1106d22f9ae89238p+1, 0x1.029a358e1e630f64p+1
457 ];
458 
459 static immutable real[8] Q3 =
460 [ -0x1.74022dd5523e6f84p-34, -0x1.2cb60d61e29ee836p-24,
461 -0x1.d19e6ec03a85e556p-17, -0x1.9ea2a7b4422f6502p-11, -0x1.c54b1e852f107162p-8,
462 0x1.e05268dd3c07989ep-3, 0x1.239c6aff14afbf82p+1, 1.0
463 ];
The OpenD Programming Language
