rmathd.normal source code

1 module rmathd.normal;
2 
3 public import rmathd.common;
4 public import rmathd.uniform;
5 
6 /*
7 ** dmd normal.d uniform.d common.d && ./normal
8 */
9 
10 T dnorm(T: double)(T x, T mu, T sigma, int log_p)
11 {
12     if (isNaN(x) || isNaN(mu) || isNaN(sigma))
13         return x + mu + sigma;
14 
15     mixin R_D__0!log_p;
16     if(!isFinite(sigma))
17         return R_D__0;
18     if(!isFinite(x) && mu == x) return T.nan;/* x-mu is NaN */
19     if (sigma <= 0) {
20         if (sigma < 0)
21             return T.nan;
22         /* sigma == 0 */
23         return (x == mu) ? T.infinity : R_D__0;
24     }
25     x = (x - mu) / sigma;
26 
27     if(!isFinite(x))
28         return R_D__0;
29 
30     x = fabs (x);
31     if (x >= 2 * sqrt(T.max))
32         return R_D__0;
33     if (log_p)
34         return -(M_LN_SQRT_2PI + 0.5 * x * x + log(sigma));
35     // more accurate, less fast :
36     if (x < 5)
37         return M_1_SQRT_2PI * exp(-0.5 * x * x) / sigma;
38 
39     if (x > sqrt(-2*M_LN2*(MIN_EXP!T + 1 - MANT_DIG!T)))
40         return 0.;
41 
42     T x1 = ldexp( nearbyint(ldexp(x, 16)), -16);
43     T x2 = x - x1;
44     return M_1_SQRT_2PI / sigma * (exp(-0.5 * x1 * x1) * exp( (-0.5 * x2 - x1) * x2 ) );
45 }
46 
47 
48 enum SIXTEN = 16;
49 
50 template swap_tail()
51 {
52     enum swap_tail = `if (x > 0.) {
53         temp = cum;
54         if(lower)
55             cum = ccum;
56         ccum = temp;
57     }`;
58 }
59 
60 template do_del(alias X)
61 {
62     enum do_del = `xsq = trunc(` ~ X.stringof ~ ` * SIXTEN) / SIXTEN;
63     del = (` ~ X.stringof ~ ` - xsq) * (` ~ X.stringof ~ ` + xsq);
64     if(log_p) {
65         cum = (-xsq * xsq * 0.5) + (-del * 0.5) + log(temp);
66         if((lower && x > 0.) || (upper && x <= 0.))
67           ccum = log1p(-exp(-xsq * xsq * 0.5) * exp(-del * 0.5) * temp);
68     }                               
69     else {
70         cum = exp(-xsq * xsq * 0.5) * exp(-del * 0.5) * temp;
71         ccum = 1.0 - cum;
72     }`;
73 }
74 
75 void pnorm_both(T: double)(T x, ref T cum, ref T ccum, int i_tail, in int log_p)
76 {
77     const static T[5] a = [
78     2.2352520354606839287,
79     161.02823106855587881,
80     1067.6894854603709582,
81     18154.981253343561249,
82     0.065682337918207449113
83     ];
84     const static T[4] b = [
85     47.20258190468824187,
86     976.09855173777669322,
87     10260.932208618978205,
88     45507.789335026729956
89     ];
90     const static T[9] c = [
91         0.39894151208813466764,
92         8.8831497943883759412,
93         93.506656132177855979,
94         597.27027639480026226,
95         2494.5375852903726711,
96         6848.1904505362823326,
97         11602.651437647350124,
98         9842.7148383839780218,
99         1.0765576773720192317e-8
100     ];
101     const static T[8] d = [
102         22.266688044328115691,
103         235.38790178262499861,
104         1519.377599407554805,
105         6485.558298266760755,
106         18615.571640885098091,
107         34900.952721145977266,
108         38912.003286093271411,
109         19685.429676859990727
110         ];
111     const static T[6] p = [
112         0.21589853405795699,
113         0.1274011611602473639,
114         0.022235277870649807,
115         0.001421619193227893466,
116         2.9112874951168792e-5,
117         0.02307344176494017303
118         ];
119     const static T[5] q = [
120         1.28426009614491121,
121         0.468238212480865118,
122         0.0659881378689285515,
123         0.00378239633202758244,
124         7.29751555083966205e-5
125         ];
126 
127     T xden, xnum, temp, del, eps, xsq, y;
128 //#ifdef NO_DENORMS
129     T min = MIN!T;
130     int i, lower, upper;
131 
132     if(isNaN(x)){
133         cum = ccum = x;
134         return;
135     }
136 
137     /* Consider changing these : */
138     eps = EPSILON!T * 0.5;
139 
140     /* i_tail in {0,1,2} =^= {lower, upper, both} */
141     lower = i_tail != 1;
142     upper = i_tail != 0;
143 
144     mixin R_D__0!log_p;
145     mixin R_D__1!log_p;
146     y = fabs(x);
147     if (y <= 0.67448975) { /* qnorm(3/4) = .6744.... -- earlier had 0.66291 */
148         if (y > eps) {
149             xsq = x * x;
150             xnum = a[4] * xsq;
151             xden = xsq;
152             for (i = 0; i < 3; ++i) {
153             xnum = (xnum + a[i]) * xsq;
154             xden = (xden + b[i]) * xsq;
155             }
156         } else xnum = xden = 0.0;
157         
158         temp = x * (xnum + a[3]) / (xden + b[3]);
159         if(lower)
160             cum = 0.5 + temp;
161         if(upper)
162             ccum = 0.5 - temp;
163         if(log_p) {
164             if(lower)
165                 cum = log(cum);
166             if(upper)
167                 ccum = log(ccum);
168         }
169     }
170     else if (y <= M_SQRT_32) {
171 
172     /* Evaluate pnorm for 0.674.. = qnorm(3/4) < |x| <= sqrt(32) ~= 5.657 */
173 
174     xnum = c[8] * y;
175     xden = y;
176     for (i = 0; i < 7; ++i) {
177         xnum = (xnum + c[i]) * y;
178         xden = (xden + d[i]) * y;
179     }
180     temp = (xnum + c[7]) / (xden + d[7]);
181         mixin (do_del!y);
182         mixin swap_tail;
183     }
184 
185 /* else   |x| > sqrt(32) = 5.657 :
186  * the next two case differentiations were really for lower=T, log=F
187  * Particularly  *not*  for  log_p !
188 
189  * Cody had (-37.5193 < x  &&  x < 8.2924) ; R originally had y < 50
190  *
191  * Note that we do want symmetry(0), lower/upper -> hence use y
192  */
193     else if((log_p && y < 1e170) 
194         || (lower && -37.5193 < x  &&  x < 8.2924)
195         || (upper && -8.2924  < x  &&  x < 37.5193)
196     ) {
197 
198     /* Evaluate pnorm for x in (-37.5, -5.657) union (5.657, 37.5) */
199     xsq = 1.0 / (x * x); /* (1./x)*(1./x) might be better */
200     xnum = p[5] * xsq;
201     xden = xsq;
202     for (i = 0; i < 4; ++i) {
203         xnum = (xnum + p[i]) * xsq;
204         xden = (xden + q[i]) * xsq;
205     }
206     temp = xsq * (xnum + p[4]) / (xden + q[4]);
207     temp = (M_1_SQRT_2PI - temp) / y;
208 
209     mixin (do_del!x);
210     mixin swap_tail;
211     } else { /* large x such that probs are 0 or 1 */
212         if(x > 0) {
213             cum = R_D__1;
214             ccum = R_D__0;
215         } else {
216                 cum = R_D__0;
217                 ccum = R_D__1;
218             }
219     }
220 
221 
222 //#ifdef NO_DENORMS
223     /* do not return "denormalized" -- we do in R */
224     if(log_p) {
225         if(cum > -min)
226             cum = -0.;
227         if(ccum > -min)
228             ccum = -0.;
229     }
230     else {
231         if(cum < min)
232             cum = 0.;
233         if(ccum < min)
234             ccum = 0.;
235     }
236 
237     return;
238 }
239 
240 T pnorm(T: double)(T x, T mu, T sigma, int lower_tail, int log_p)
241 {
242     T p, cp;
243 
244     /* Note: The structure of these checks has been carefully thought through.
245      * For example, if x == mu and sigma == 0, we get the correct answer 1.
246      */
247     if(isNaN(x) || isNaN(mu) || isNaN(sigma))
248         return x + mu + sigma;
249 
250     if(!x.isFinite && mu == x)
251         return T.nan;
252 
253     if (sigma <= 0) {
254         if(sigma < 0)
255             return T.nan;
256         /* sigma = 0 : */
257         return (x < mu) ? R_DT_0!T(lower_tail, log_p) : R_DT_1!T(lower_tail, log_p);
258     }
259     p = (x - mu) / sigma;
260     if(!p.isFinite)
261         return (x < mu) ? R_DT_0!T(lower_tail, log_p) : R_DT_1!T(lower_tail, log_p);
262     x = p;
263     pnorm_both(x, p, cp, (lower_tail ? 0 : 1), log_p);
264     return(lower_tail ? p : cp);
265 }
266 
267 /*template R_Q_P01_boundaries(alias p, alias LEFT, alias RIGHT)
268 {
269     enum R_Q_P01_boundaries = `if (log_p) {
270         if(p > 0)
271             return T.nan;
272         if(p == 0)
273             return lower_tail ? ` ~ RIGHT.stringof ~ ` : ` ~ LEFT.stringof ~ `;
274         if(p == -T.infinity)
275             return lower_tail ? ` ~ LEFT.stringof ~ ` : ` ~ RIGHT.stringof ~ `;
276     } else {
277         if(p < 0 || p > 1)
278             return T.nan;
279         if(p == 0)
280             return lower_tail ? ` ~ LEFT.stringof ~` : ` ~ RIGHT.stringof ~ `;
281         if(p == 1)
282             return lower_tail ? ` ~ RIGHT.stringof ~ ` : ` ~ LEFT.stringof ~ `;
283     }`;
284 }*/
285 
286 
287 T qnorm(T: double)(T p, T mu, T sigma, int lower_tail, int log_p)
288 {
289     T p_, q, r, val;
290 
291 
292     /* Convert this code snippet back to a mixin noe that you know how! */
293     if (isNaN(p) || isNaN(mu) || isNaN(sigma))
294         return p + mu + sigma;
295 
296     immutable(T) PINF = T.infinity;
297     immutable(T) NINF = -T.infinity;
298     mixin (R_Q_P01_boundaries!(p, NINF, PINF)());
299 
300     if(sigma  < 0)
301         return T.nan;
302     if(sigma == 0)
303         return mu;
304 
305     mixin R_DT_qIv!p;
306     p_ = R_DT_qIv;/* real lower_tail prob. p */
307     q = p_ - 0.5;
308 
309 /*-- use AS 241 --- */
310 /* double ppnd16_(double *p, long *ifault)*/
311 /*      ALGORITHM AS241  APPL. STATIST. (1988) VOL. 37, NO. 3
312 
313         Produces the normal deviate Z corresponding to a given lower
314         tail area of P; Z is accurate to about 1 part in 10**16.
315 
316         (original fortran code used PARAMETER(..) for the coefficients
317          and provided hash codes for checking them...)
318 */
319     if (fabs(q) <= .425) {/* 0.075 <= p <= 0.925 */
320         r = .180625 - q * q;
321     val =
322             q * (((((((r * 2509.0809287301226727 +
323                        33430.575583588128105) * r + 67265.770927008700853) * r +
324                      45921.953931549871457) * r + 13731.693765509461125) * r +
325                    1971.5909503065514427) * r + 133.14166789178437745) * r +
326                  3.387132872796366608)
327             / (((((((r * 5226.495278852854561 +
328                      28729.085735721942674) * r + 39307.89580009271061) * r +
329                    21213.794301586595867) * r + 5394.1960214247511077) * r +
330                  687.1870074920579083) * r + 42.313330701600911252) * r + 1.);
331     }
332     else { /* closer than 0.075 from {0,1} boundary */
333 
334     /* r = min(p, 1-p) < 0.075 */
335     mixin R_DT_CIv!p;
336     if (q > 0)
337         r = R_DT_CIv;/* 1-p */
338     else
339         r = p_;/* = R_DT_Iv(p) ^=  p */
340 
341     r = sqrt(- ((log_p &&
342              ((lower_tail && q <= 0) || (!lower_tail && q > 0))) ?
343             p : /* else */ log(r)));
344         /* r = sqrt(-log(r))  <==>  min(p, 1-p) = exp( - r^2 ) */
345 
346         if (r <= 5.) { /* <==> min(p,1-p) >= exp(-25) ~= 1.3888e-11 */
347             r += -1.6;
348             val = (((((((r * 7.7454501427834140764e-4 +
349                        .0227238449892691845833) * r + .24178072517745061177) *
350                      r + 1.27045825245236838258) * r +
351                     3.64784832476320460504) * r + 5.7694972214606914055) *
352                   r + 4.6303378461565452959) * r +
353                  1.42343711074968357734)
354                 / (((((((r *
355                          1.05075007164441684324e-9 + 5.475938084995344946e-4) *
356                         r + .0151986665636164571966) * r +
357                        .14810397642748007459) * r + .68976733498510000455) *
358                      r + 1.6763848301838038494) * r +
359                     2.05319162663775882187) * r + 1.);
360         }
361         else { /* very close to  0 or 1 */
362             r += -5.;
363             val = (((((((r * 2.01033439929228813265e-7 +
364                        2.71155556874348757815e-5) * r +
365                       .0012426609473880784386) * r + .026532189526576123093) *
366                     r + .29656057182850489123) * r +
367                    1.7848265399172913358) * r + 5.4637849111641143699) *
368                  r + 6.6579046435011037772)
369                 / (((((((r *
370                          2.04426310338993978564e-15 + 1.4215117583164458887e-7)*
371                         r + 1.8463183175100546818e-5) * r +
372                        7.868691311456132591e-4) * r + .0148753612908506148525)
373                      * r + .13692988092273580531) * r +
374                     .59983220655588793769) * r + 1.);
375         }
376 
377     if(q < 0.0)
378         val = -val;
379         /* return (q >= 0.)? r : -r ;*/
380     }
381     return mu + sigma * val;
382 }
383 
384 double BM_norm_keep = 0.0;
385 N01type N01_kind = INVERSION;
386 
387 alias void* function() DL_FUNC;
388 extern DL_FUNC User_norm_fun;
389 
390 enum C1 = 0.398942280401433;
391 enum C2 = 0.180025191068563;
392 
393 T norm_rand(T: double)()
394 {
395 
396     const static T[32] a =
397     [
398     0.0000000, 0.03917609, 0.07841241, 0.1177699,
399     0.1573107, 0.19709910, 0.23720210, 0.2776904,
400     0.3186394, 0.36012990, 0.40225010, 0.4450965,
401     0.4887764, 0.53340970, 0.57913220, 0.6260990,
402     0.6744898, 0.72451440, 0.77642180, 0.8305109,
403     0.8871466, 0.94678180, 1.00999000, 1.0775160,
404     1.1503490, 1.22985900, 1.31801100, 1.4177970,
405     1.5341210, 1.67594000, 1.86273200, 2.1538750
406     ];
407 
408     const static T[31] d =
409     [
410     0.0000000, 0.0000000, 0.0000000, 0.0000000,
411     0.0000000, 0.2636843, 0.2425085, 0.2255674,
412     0.2116342, 0.1999243, 0.1899108, 0.1812252,
413     0.1736014, 0.1668419, 0.1607967, 0.1553497,
414     0.1504094, 0.1459026, 0.1417700, 0.1379632,
415     0.1344418, 0.1311722, 0.1281260, 0.1252791,
416     0.1226109, 0.1201036, 0.1177417, 0.1155119,
417     0.1134023, 0.1114027, 0.1095039
418     ];
419 
420     const static T[31] t =
421     [
422     7.673828e-4, 0.002306870, 0.003860618, 0.005438454,
423     0.007050699, 0.008708396, 0.010423570, 0.012209530,
424     0.014081250, 0.016055790, 0.018152900, 0.020395730,
425     0.022811770, 0.025434070, 0.028302960, 0.031468220,
426     0.034992330, 0.038954830, 0.043458780, 0.048640350,
427     0.054683340, 0.061842220, 0.070479830, 0.081131950,
428     0.094624440, 0.112300100, 0.136498000, 0.171688600,
429     0.227624100, 0.330498000, 0.584703100
430     ];
431 
432     const static T[31] h =
433     [
434     0.03920617, 0.03932705, 0.03950999, 0.03975703,
435     0.04007093, 0.04045533, 0.04091481, 0.04145507,
436     0.04208311, 0.04280748, 0.04363863, 0.04458932,
437     0.04567523, 0.04691571, 0.04833487, 0.04996298,
438     0.05183859, 0.05401138, 0.05654656, 0.05953130,
439     0.06308489, 0.06737503, 0.07264544, 0.07926471,
440     0.08781922, 0.09930398, 0.11555990, 0.14043440,
441     0.18361420, 0.27900160, 0.70104740
442     ];
443 
444     /*----------- Constants and definitions for  Kinderman - Ramage --- */
445     /*
446      *  REFERENCE
447      *
448      *    Kinderman A. J. and Ramage J. G. (1976).
449      *    Computer generation of normal random variables.
450      *    JASA 71, 893-896.
451      */
452 
453     const static T A =  2.216035867166471;
454 
455     T s, u1, w, y, u2, u3, aa, tt, theta, R;
456     int i;
457 
458 
459     template g(T)
460     {
461         auto g(T x)
462         {
463             return C1 * exp(-x*x/2.0) - C2*(A - x);
464         }
465     }
466 
467     switch(N01_kind) {
468 
469     case  AHRENS_DIETER: /* see Reference above */
470 
471     u1 = unif_rand!T();
472     s = 0.0;
473     if (u1 > 0.5)
474         s = 1.0;
475     u1 = u1 + u1 - s;
476     u1 *= 32.0;
477     i = cast(int) u1;
478     if (i == 32)
479         i = 31;
480     if (i != 0) {
481         u2 = u1 - i;
482         aa = a[i - 1];
483         while (u2 <= t[i - 1]) {
484             u1 = unif_rand!T();
485             w = u1 * (a[i] - aa);
486             tt = (w * 0.5 + aa) * w;
487             for(;;) {
488                 if (u2 > tt)
489                     goto deliver;
490                 u1 = unif_rand!T();
491                 if (u2 < u1)
492                     break;
493                 tt = u1;
494                 u2 = unif_rand!T();
495             }
496             u2 = unif_rand!T();
497         }
498         w = (u2 - t[i - 1]) * h[i - 1];
499     } else {
500         i = 6;
501         aa = a[31];
502         for(;;) {
503         u1 = u1 + u1;
504         if (u1 >= 1.0)
505             break;
506         aa = aa + d[i - 1];
507         i = i + 1;
508         }
509         u1 = u1 - 1.0;
510         for(;;) {
511         w = u1 * d[i - 1];
512         tt = (w * 0.5 + aa) * w;
513         for(;;) {
514             u2 = unif_rand!T();
515             if (u2 > tt)
516                 goto jump;
517             u1 = unif_rand!T();
518             if (u2 < u1)
519                 break;
520             tt = u1;
521         }
522         u1 = unif_rand!T();
523         }
524       jump:;
525     }
526 
527       deliver:
528     y = aa + w;
529     return (s == 1.0) ? -y : y;
530 
531     /*-----------------------------------------------------------*/
532 
533     case BUGGY_KINDERMAN_RAMAGE: /* see Reference above */
534     /* note: this has problems, but is retained for
535      * reproducibility of older codes, with the same
536      * numeric code */
537     u1 = unif_rand!T();
538     if(u1 < 0.884070402298758) {
539         u2 = unif_rand!T();
540         return A*(1.13113163544180 * u1 + u2 - 1);
541     }
542 
543     if(u1 >= 0.973310954173898) { /* tail: */
544         for(;;) {
545         u2 = unif_rand!T();
546         u3 = unif_rand!T();
547         tt = (A*A - 2*log(u3));
548         if( (u2*u2) < (A*A)/tt )
549             return (u1 < 0.986655477086949) ? sqrt(tt) : -sqrt(tt);
550         }
551     }
552 
553     if(u1 >= 0.958720824790463) { /* region3: */
554         for(;;) {
555         u2 = unif_rand!T();
556         u3 = unif_rand!T();
557         tt = A - 0.630834801921960 * fmin2(u2, u3);
558         if(fmax2(u2, u3) <= 0.755591531667601)
559             return (u2 < u3) ? tt : -tt;
560         if(0.034240503750111 * fabs(u2 - u3) <= g(tt))
561             return (u2 < u3) ? tt : -tt;
562         }
563     }
564 
565     if(u1 >= 0.911312780288703) { /* region2: */
566         for(;;) {
567         u2 = unif_rand!T();
568         u3 = unif_rand!T();
569         tt = 0.479727404222441 + 1.105473661022070 * fmin2(u2, u3);
570         if( fmax2(u2, u3) <= 0.872834976671790 )
571             return (u2 < u3) ? tt : -tt;
572         if( 0.049264496373128 * fabs(u2 - u3) <= g(tt) )
573             return (u2 < u3) ? tt : -tt;
574         }
575     }
576 
577     /* ELSE  region1: */
578     for(;;) {
579         u2 = unif_rand!T();
580         u3 = unif_rand!T();
581         tt = 0.479727404222441 - 0.595507138015940 * fmin2(u2, u3);
582         if(fmax2(u2, u3) <= 0.805577924423817)
583             return (u2 < u3) ? tt : -tt;
584     }
585     case BOX_MULLER:
586     if(BM_norm_keep != 0.0) { /* An exact test is intentional */
587         s = BM_norm_keep;
588         BM_norm_keep = 0.0;
589         return s;
590     } else {
591         theta = 2 * M_PI * unif_rand!T();
592         R = sqrt(-2 * log(unif_rand!T())) + 10 * MIN!T; /* ensure non-zero */
593         BM_norm_keep = R * sin(theta);
594         return R * cos(theta);
595     }
596 
597     //case USER_NORM:
598     //return cast(T)User_norm_fun();
599 
600     case INVERSION:
601     enum BIG = 134217728; /* 2^27 */
602     /* unif_rand() alone is not of high enough precision */
603     u1 = unif_rand!T();
604     u1 = cast(int)(BIG*u1) + unif_rand!T();
605     return qnorm(u1/BIG, 0.0, 1.0, 1, 0);
606     case KINDERMAN_RAMAGE: /* see Reference above */
607     /* corrected version from Josef Leydold
608      * */
609     u1 = unif_rand!T();
610     if(u1 < 0.884070402298758) {
611         u2 = unif_rand!T();
612         return A*(1.131131635444180*u1 + u2 - 1);
613     }
614 
615     if(u1 >= 0.973310954173898) { /* tail: */
616         for(;;) {
617             u2 = unif_rand!T();
618             u3 = unif_rand!T();
619             tt = (A*A - 2*log(u3));
620             if( u2*u2<(A*A)/tt )
621                 return (u1 < 0.986655477086949) ? sqrt(tt) : -sqrt(tt);
622         }
623     }
624 
625     if(u1 >= 0.958720824790463) { /* region3: */
626         for(;;) {
627         u2 = unif_rand!T();
628         u3 = unif_rand!T();
629         tt = A - 0.630834801921960*fmin2(u2, u3);
630         if(fmax2(u2, u3) <= 0.755591531667601)
631             return (u2 < u3) ? tt : -tt;
632         if(0.034240503750111*fabs(u2 - u3) <= g(tt))
633             return (u2 < u3) ? tt : -tt;
634         }
635     }
636 
637     if(u1 >= 0.911312780288703) { /* region2: */
638         for(;;) {
639             u2 = unif_rand!T();
640             u3 = unif_rand!T();
641             tt = 0.479727404222441+1.105473661022070*fmin2(u2, u3);
642         if( fmax2(u2,u3)<=0.872834976671790 )
643             return (u2 < u3) ? tt : -tt;
644         if(0.049264496373128*fabs(u2 - u3)<=g(tt))
645             return (u2 < u3) ? tt : -tt;
646         }
647     }
648 
649     /* ELSE  region1: */
650     for(;;) {
651         u2 = unif_rand!T();
652         u3 = unif_rand!T();
653         tt = 0.479727404222441-0.595507138015940*fmin2(u2, u3);
654         if (tt < 0.) continue;
655         if(fmax2(u2, u3) <= 0.805577924423817)
656         return (u2 < u3) ? tt : -tt;
657             if(0.053377549506886*fabs(u2 - u3) <= g(tt))
658         return (u2 < u3) ? tt : -tt;
659     }
660     default:
661         //assert(0, "norm_rand(): invalid N01_kind: " ~ N01_kind.stringof);
662         return 0.0;/*- -Wall */
663     }/*switch*/
664 }
665 
666 
667 T rnorm(T: double)(T mu, T sigma)
668 {
669     if (isNaN(mu) || !isFinite(sigma) || sigma < 0.)
670         return T.nan;
671     if (sigma == 0. || !isFinite(mu))
672         return mu;
673     else
674         return mu + sigma * norm_rand!T();
675 }
676 
677 
678 
679 void test_normal()
680 {
681     import std.stdio: write, writeln;
682     writeln("dnorm: ", "0: ", dnorm(0., 0., 1., 0), ", 1.96: ", dnorm(1.96, 0., 1., 0), 
683             ", 1.644854: ", dnorm(1.644854, 0., 1., 0));
684     writeln("pnorm: ", "z = 1.96: ", pnorm(1.96, 0, 1, 0, 0),
685         ", z = 1.644854: ", pnorm(1.644854, 0, 1, 0, 0));
686     writeln("qnorm: ", "p = 0.975: ", qnorm(0.975, 0, 1, 0, 0),
687         ", p = 0.95: ", qnorm(0.95, 0, 1, 0, 0));
688     writeln("rnorm(0., 1.): ", rnorm(0., 1.), "; rnorm(0., 1.): ", rnorm(0., 1.),
689         "; rnorm(0., 1.): ", rnorm(0., 1.), "; rnorm(0., 1.): ", rnorm(0., 1.));
690 }
691 
692 
693 
694 
695