gfm.math.simplerng source code

1 /**
2  Translation of SimpleRNG.
3  Removed the builtin RNG to use std.random, but kept the distribution functions.
4  John D. Cook confirmed this code as public domain.
5 
6  Authors: John D. Cook.
7  See_also: $(WEB www.johndcook.com/cpp_random_number_generation.html)
8  */
9 module gfm.math.simplerng;
10 
11 public import std.random;
12 import std.math;
13 
14 /// Returns: Normal (Gaussian) random sample.
15 /// See_also: Box-Muller algorithm.
16 double randNormal(RNG)(ref RNG rng, double mean = 0.0, double standardDeviation = 1.0)
17 {
18     assert(standardDeviation > 0);
19     double u1;
20 
21     do
22     {
23         u1 = uniform(0.0, 1.0, rng);
24     } while (u1 == 0); // u1 must not be zero
25     double u2 = uniform(0.0, 1.0, rng);
26     double r = sqrt(-2.0 * log(u1));
27     double theta = 2.0 * PI * u2;
28     return mean + standardDeviation * r * sin(theta);
29 }
30 
31 /// Returns: Exponential random sample with specified mean.
32 double randExponential(RNG)(ref RNG rng, double mean = 1.0)
33 {
34     assert(mean > 0);
35     return -mean*log(uniform(0.0, 1.0, rng));
36 }
37 
38 /// Returns: Gamma random sample.
39 /// See_also: "A Simple Method for Generating Gamma Variables"
40 /// by George Marsaglia and Wai Wan Tsang.  ACM Transactions on Mathematical Software
41 /// Vol 26, No 3, September 2000, pages 363-372.
42 double randGamma(RNG)(ref RNG rng, double shape, double scale)
43 {
44     double d, c, x, xsquared, v, u;
45 
46     if (shape >= 1.0)
47     {
48         d = shape - 1.0/3.0;
49         c = 1.0/sqrt(9.0*d);
50         for (;;)
51         {
52             do
53             {
54                 x = randNormal(rng);
55                 v = 1.0 + c*x;
56             }
57             while (v <= 0.0);
58             v = v*v*v;
59             u = uniform(0.0, 1.0, rng);
60             xsquared = x*x;
61             if (u < 1.0 -.0331*xsquared*xsquared || log(u) < 0.5*xsquared + d*(1.0 - v + log(v)))
62                 return scale*d*v;
63         }
64     }
65     else
66     {
67         assert(shape > 0);
68         double g = randGamma(rng, shape+1.0, 1.0);
69         double w = uniform(0.0, 1.0, rng);
70         return scale*g*pow(w, 1.0/shape);
71     }
72 }
73 
74 /// Returns: Chi-square sample.
75 double randChiSquare(RNG)(ref RNG rng, double degreesOfFreedom)
76 {
77     // A chi squared distribution with n degrees of freedom
78     // is a gamma distribution with shape n/2 and scale 2.
79     return randGamma(rng, 0.5 * degreesOfFreedom, 2.0);
80 }
81 
82 /// Returns: Inverse-gamma sample.
83 double randInverseGamma(RNG)(ref RNG rng, double shape, double scale)
84 {
85     // If X is gamma(shape, scale) then
86     // 1/Y is inverse gamma(shape, 1/scale)
87     return 1.0 / randGamma(rng, shape, 1.0 / scale);
88 }
89 
90 /// Returns: Weibull sample.
91 double randWeibull(RNG)(ref RNG rng, double shape, double scale)
92 {
93     assert(shape > 0 && scale > 0);
94     return scale * pow(-log(uniform(0.0, 1.0, rng)), 1.0 / shape);
95 }
96 
97 /// Returns: Cauchy sample.
98 double randCauchy(RNG)(ref RNG rng, double median, double scale)
99 {
100     assert(scale > 0);
101     double p = uniform(0.0, 1.0, rng);
102 
103     // Apply inverse of the Cauchy distribution function to a uniform
104     return median + scale*tan(PI*(p - 0.5));
105 }
106 
107 /// Returns: Student-t sample.
108 /// See_also: Seminumerical Algorithms by Knuth.
109 double randStudentT(RNG)(ref RNG rng, double degreesOfFreedom)
110 {
111     assert(degreesOfFreedom > 0);
112 
113 
114     double y1 = getNormal(rng);
115     double y2 = getChiSquare(rng, degreesOfFreedom);
116     return y1 / sqrt(y2 / degreesOfFreedom);
117 }
118 
119 /// Returns: Laplace distribution random sample (also known as the double exponential distribution).
120 double randLaplace(RNG)(ref RNG rng, double mean, double scale)
121 {
122     double u = uniform(0.0, 1.0, rng);
123     return (u < 0.5) ? (mean + scale*log(2.0*u))
124                         : (mean - scale*log(2*(1-u)));
125 }
126 
127 /// Returns: Log-normal sample.
128 double randLogNormal(RNG)(ref RNG rng, double mu, double sigma)
129 {
130     return exp(getNormal(rng, mu, sigma));
131 }
132 
133 /// Returns: Beta sample.
134 double randBeta(RNG)(ref RNG rng, double a, double b)
135 {
136     assert(a > 0 && b > 0);
137 
138     // There are more efficient methods for generating beta samples.
139     // However such methods are a little more efficient and much more complicated.
140     // For an explanation of why the following method works, see
141     // http://www.johndcook.com/distribution_chart.html#gamma_beta
142 
143     double u = getGamma(rng, a, 1.0);
144     double v = getGamma(rng, b, 1.0);
145     return u / (u + v);
146 }
147 
148 /// Returns: Poisson sample.
149 int randPoisson(RNG)(ref RNG rng, double lambda)
150 {
151     return (lambda < 30.0) ? poissonSmall(rng, lambda) : poissonLarge(rng, lambda);
152 }
153 
154 private
155 {
156     int poissonSmall(RNG)(ref RNG rng, double lambda)
157     {
158         // Algorithm due to Donald Knuth, 1969.
159         double p = 1.0, L = exp(-lambda);
160         int k = 0;
161         do
162         {
163             k++;
164             p *= uniform(0.0, 1.0, rng);
165         }
166         while (p > L);
167         return k - 1;
168     }
169 
170     int poissonLarge(RNG)(ref RNG rng, double lambda)
171     {
172         // "Rejection method PA" from "The Computer Generation of Poisson Random Variables" by A. C. Atkinson
173         // Journal of the Royal Statistical Society Series C (Applied Statistics) Vol. 28, No. 1. (1979)
174         // The article is on pages 29-35. The algorithm given here is on page 32.
175 
176         double c = 0.767 - 3.36/lambda;
177         double beta = PI/sqrt(3.0*lambda);
178         double alpha = beta*lambda;
179         double k = log(c) - lambda - log(beta);
180 
181         for(;;)
182         {
183             double u = uniform(0.0, 1.0, rng);
184             double x = (alpha - log((1.0 - u)/u))/beta;
185             int n = cast(int)(floor(x + 0.5));
186             if (n < 0)
187                 continue;
188             double v = uniform(0.0, 1.0, rng);
189             double y = alpha - beta*x;
190             double temp = 1.0 + exp(y);
191             double lhs = y + log(v/(temp*temp));
192             double rhs = k + n*log(lambda) - logFactorial(n);
193             if (lhs <= rhs)
194                 return n;
195         }
196     }
197 
198     double logFactorial(int n) nothrow
199     {
200         assert(n >= 0);
201         if (n > 254)
202         {
203             double x = n + 1;
204             return (x - 0.5)*log(x) - x + 0.5*log(2*PI) + 1.0/(12.0*x);
205         }
206         else
207         {
208             return LOG_FACTORIAL[n];
209         }
210     }
211 }
212 
213 private static immutable double[255] LOG_FACTORIAL =
214 [
215     0.000000000000000,
216     0.000000000000000,
217     0.693147180559945,
218     1.791759469228055,
219     3.178053830347946,
220     4.787491742782046,
221     6.579251212010101,
222     8.525161361065415,
223     10.604602902745251,
224     12.801827480081469,
225     15.104412573075516,
226     17.502307845873887,
227     19.987214495661885,
228     22.552163853123421,
229     25.191221182738683,
230     27.899271383840894,
231     30.671860106080675,
232     33.505073450136891,
233     36.395445208033053,
234     39.339884187199495,
235     42.335616460753485,
236     45.380138898476908,
237     48.471181351835227,
238     51.606675567764377,
239     54.784729398112319,
240     58.003605222980518,
241     61.261701761002001,
242     64.557538627006323,
243     67.889743137181526,
244     71.257038967168000,
245     74.658236348830158,
246     78.092223553315307,
247     81.557959456115029,
248     85.054467017581516,
249     88.580827542197682,
250     92.136175603687079,
251     95.719694542143202,
252     99.330612454787428,
253     102.968198614513810,
254     106.631760260643450,
255     110.320639714757390,
256     114.034211781461690,
257     117.771881399745060,
258     121.533081515438640,
259     125.317271149356880,
260     129.123933639127240,
261     132.952575035616290,
262     136.802722637326350,
263     140.673923648234250,
264     144.565743946344900,
265     148.477766951773020,
266     152.409592584497350,
267     156.360836303078800,
268     160.331128216630930,
269     164.320112263195170,
270     168.327445448427650,
271     172.352797139162820,
272     176.395848406997370,
273     180.456291417543780,
274     184.533828861449510,
275     188.628173423671600,
276     192.739047287844900,
277     196.866181672889980,
278     201.009316399281570,
279     205.168199482641200,
280     209.342586752536820,
281     213.532241494563270,
282     217.736934113954250,
283     221.956441819130360,
284     226.190548323727570,
285     230.439043565776930,
286     234.701723442818260,
287     238.978389561834350,
288     243.268849002982730,
289     247.572914096186910,
290     251.890402209723190,
291     256.221135550009480,
292     260.564940971863220,
293     264.921649798552780,
294     269.291097651019810,
295     273.673124285693690,
296     278.067573440366120,
297     282.474292687630400,
298     286.893133295426990,
299     291.323950094270290,
300     295.766601350760600,
301     300.220948647014100,
302     304.686856765668720,
303     309.164193580146900,
304     313.652829949878990,
305     318.152639620209300,
306     322.663499126726210,
307     327.185287703775200,
308     331.717887196928470,
309     336.261181979198450,
310     340.815058870798960,
311     345.379407062266860,
312     349.954118040770250,
313     354.539085519440790,
314     359.134205369575340,
315     363.739375555563470,
316     368.354496072404690,
317     372.979468885689020,
318     377.614197873918670,
319     382.258588773060010,
320     386.912549123217560,
321     391.575988217329610,
322     396.248817051791490,
323     400.930948278915760,
324     405.622296161144900,
325     410.322776526937280,
326     415.032306728249580,
327     419.750805599544780,
328     424.478193418257090,
329     429.214391866651570,
330     433.959323995014870,
331     438.712914186121170,
332     443.475088120918940,
333     448.245772745384610,
334     453.024896238496130,
335     457.812387981278110,
336     462.608178526874890,
337     467.412199571608080,
338     472.224383926980520,
339     477.044665492585580,
340     481.872979229887900,
341     486.709261136839360,
342     491.553448223298010,
343     496.405478487217580,
344     501.265290891579240,
345     506.132825342034830,
346     511.008022665236070,
347     515.890824587822520,
348     520.781173716044240,
349     525.679013515995050,
350     530.584288294433580,
351     535.496943180169520,
352     540.416924105997740,
353     545.344177791154950,
354     550.278651724285620,
355     555.220294146894960,
356     560.169054037273100,
357     565.124881094874350,
358     570.087725725134190,
359     575.057539024710200,
360     580.034272767130800,
361     585.017879388839220,
362     590.008311975617860,
363     595.005524249382010,
364     600.009470555327430,
365     605.020105849423770,
366     610.037385686238740,
367     615.061266207084940,
368     620.091704128477430,
369     625.128656730891070,
370     630.172081847810200,
371     635.221937855059760,
372     640.278183660408100,
373     645.340778693435030,
374     650.409682895655240,
375     655.484856710889060,
376     660.566261075873510,
377     665.653857411105950,
378     670.747607611912710,
379     675.847474039736880,
380     680.953419513637530,
381     686.065407301994010,
382     691.183401114410800,
383     696.307365093814040,
384     701.437263808737160,
385     706.573062245787470,
386     711.714725802289990,
387     716.862220279103440,
388     722.015511873601330,
389     727.174567172815840,
390     732.339353146739310,
391     737.509837141777440,
392     742.685986874351220,
393     747.867770424643370,
394     753.055156230484160,
395     758.248113081374300,
396     763.446610112640200,
397     768.650616799717000,
398     773.860102952558460,
399     779.075038710167410,
400     784.295394535245690,
401     789.521141208958970,
402     794.752249825813460,
403     799.988691788643450,
404     805.230438803703120,
405     810.477462875863580,
406     815.729736303910160,
407     820.987231675937890,
408     826.249921864842800,
409     831.517780023906310,
410     836.790779582469900,
411     842.068894241700490,
412     847.352097970438420,
413     852.640365001133090,
414     857.933669825857460,
415     863.231987192405430,
416     868.535292100464630,
417     873.843559797865740,
418     879.156765776907600,
419     884.474885770751830,
420     889.797895749890240,
421     895.125771918679900,
422     900.458490711945270,
423     905.796028791646340,
424     911.138363043611210,
425     916.485470574328820,
426     921.837328707804890,
427     927.193914982476710,
428     932.555207148186240,
429     937.921183163208070,
430     943.291821191335660,
431     948.667099599019820,
432     954.046996952560450,
433     959.431492015349480,
434     964.820563745165940,
435     970.214191291518320,
436     975.612353993036210,
437     981.015031374908400,
438     986.422203146368590,
439     991.833849198223450,
440     997.249949600427840,
441     1002.670484599700300,
442     1008.095434617181700,
443     1013.524780246136200,
444     1018.958502249690200,
445     1024.396581558613400,
446     1029.838999269135500,
447     1035.285736640801600,
448     1040.736775094367400,
449     1046.192096209724900,
450     1051.651681723869200,
451     1057.115513528895000,
452     1062.583573670030100,
453     1068.055844343701400,
454     1073.532307895632800,
455     1079.012946818975000,
456     1084.497743752465600,
457     1089.986681478622400,
458     1095.479742921962700,
459     1100.976911147256000,
460     1106.478169357800900,
461     1111.983500893733000,
462     1117.492889230361000,
463     1123.006317976526100,
464     1128.523770872990800,
465     1134.045231790853000,
466     1139.570684729984800,
467     1145.100113817496100,
468     1150.633503306223700,
469     1156.170837573242400,
470 ];
471 
472 unittest
473 {
474     Xorshift32 rng;
475     rng.seed(unpredictableSeed());
476 
477     double x = randNormal!Xorshift32(rng, 0.0, 1.0);
478     x = randExponential!Xorshift32(rng);
479     x = randGamma!Xorshift32(rng, 1.2, 1.0);
480     x = randGamma!Xorshift32(rng, 0.8, 2.0);
481     x = randChiSquare!Xorshift32(rng, 2.0);
482     x = randInverseGamma!Xorshift32(rng, 1.1, 0.7);
483     x = randWeibull!Xorshift32(rng, 3.0, 0.7);
484     x = randCauchy!Xorshift32(rng, 5.0, 1.4);
485 }