1 | /// <summary> |
2 | /// SimpleRNG is a simple random number generator based on |
3 | /// George Marsaglia's MWC (multiply with carry) generator. |
4 | /// Although it is very simple, it passes Marsaglia's DIEHARD |
5 | /// series of random number generator tests. |
6 | /// |
7 | /// Written by John D. Cook |
8 | /// http://www.johndcook.com |
9 | /// Converted to C++ by Pietro Incardona |
10 | /// </summary> |
11 | |
12 | #ifndef SIMPLE_RNG_HPP |
13 | #define SIMPLE_RNG_HPP |
14 | |
15 | #include <string> |
16 | #include <cmath> |
17 | |
18 | class SimpleRNG |
19 | { |
20 | private: |
21 | unsigned int m_w; |
22 | unsigned int m_z; |
23 | |
24 | // This is the heart of the generator. |
25 | // It uses George Marsaglia's MWC algorithm to produce an unsigned integer. |
26 | // See http://www.bobwheeler.com/statistics/Password/MarsagliaPost.txt |
27 | unsigned int GetUint() |
28 | { |
29 | m_z = 36969 * (m_z & 65535) + (m_z >> 16); |
30 | m_w = 18000 * (m_w & 65535) + (m_w >> 16); |
31 | return (m_z << 16) + m_w; |
32 | } |
33 | |
34 | public: |
35 | |
36 | SimpleRNG() |
37 | { |
38 | // These values are not magical, just the default values Marsaglia used. |
39 | // Any pair of unsigned integers should be fine. |
40 | m_w = 521288629; |
41 | m_z = 362436069; |
42 | } |
43 | |
44 | // The random generator seed can be set three ways: |
45 | // 1) specifying two non-zero unsigned integers |
46 | // 2) specifying one non-zero unsigned integer and taking a default value for the second |
47 | // 3) setting the seed from the system time |
48 | |
49 | void SetSeed(unsigned int u, unsigned int v) |
50 | { |
51 | if (u != 0) m_w = u; |
52 | if (v != 0) m_z = v; |
53 | } |
54 | |
55 | void SetSeed(unsigned int u) |
56 | { |
57 | m_w = u; |
58 | } |
59 | |
60 | void SetSeedFromSystemTime() |
61 | { |
62 | long x = clock(); |
63 | SetSeed((unsigned int)(x >> 16), (unsigned int)(x % 4294967296)); |
64 | } |
65 | |
66 | // Produce a uniform random sample from the open interval (0, 1). |
67 | // The method will not return either end point. |
68 | double GetUniform() |
69 | { |
70 | // 0 <= u < 2^32 |
71 | unsigned int u = GetUint(); |
72 | // The magic number below is 1/(2^32 + 2). |
73 | // The result is strictly between 0 and 1. |
74 | return (u + 1.0) * 2.328306435454494e-10; |
75 | } |
76 | |
77 | // Get normal (Gaussian) random sample with mean 0 and standard deviation 1 |
78 | double GetNormal() |
79 | { |
80 | // Use Box-Muller algorithm |
81 | double u1 = GetUniform(); |
82 | double u2 = GetUniform(); |
83 | double r = std::sqrt( -2.0*std::log(u1) ); |
84 | double theta = 2.0*3.14159265359*u2; |
85 | return r*std::sin(theta); |
86 | } |
87 | |
88 | // Get normal (Gaussian) random sample with specified mean and standard deviation |
89 | double GetNormal(double mean, double standardDeviation) |
90 | { |
91 | if (standardDeviation <= 0.0) |
92 | { |
93 | std::cerr << __FILE__ << ":" << __LINE__ << " Shape must be positive. Received." << std::to_string(standardDeviation); |
94 | throw 100001; |
95 | } |
96 | return mean + standardDeviation*GetNormal(); |
97 | } |
98 | |
99 | // Get exponential random sample with mean 1 |
100 | double GetExponential() |
101 | { |
102 | return -std::log( GetUniform() ); |
103 | } |
104 | |
105 | // Get exponential random sample with specified mean |
106 | double GetExponential(double mean) |
107 | { |
108 | if (mean <= 0.0) |
109 | { |
110 | std::cerr << __FILE__ << ":" << __LINE__ << "Mean must be positive. Received " << mean; |
111 | throw 100002; |
112 | } |
113 | return mean*GetExponential(); |
114 | } |
115 | |
116 | double GetGamma(double shape, double scale) |
117 | { |
118 | // Implementation based on "A Simple Method for Generating Gamma Variables" |
119 | // by George Marsaglia and Wai Wan Tsang. ACM Transactions on Mathematical Software |
120 | // Vol 26, No 3, September 2000, pages 363-372. |
121 | |
122 | double d, c, x, xsquared, v, u; |
123 | |
124 | if (shape >= 1.0) |
125 | { |
126 | d = shape - 1.0/3.0; |
127 | c = 1.0/std::sqrt(9.0*d); |
128 | for (;;) |
129 | { |
130 | do |
131 | { |
132 | x = GetNormal(); |
133 | v = 1.0 + c*x; |
134 | } |
135 | while (v <= 0.0); |
136 | v = v*v*v; |
137 | u = GetUniform(); |
138 | xsquared = x*x; |
139 | if (u < 1.0 -.0331*xsquared*xsquared || std::log(u) < 0.5*xsquared + d*(1.0 - v + std::log(v))) |
140 | return scale*d*v; |
141 | } |
142 | } |
143 | else if (shape <= 0.0) |
144 | { |
145 | std::cerr << __FILE__ << ":" << __LINE__ << " Shape must be positive. Received" << shape << "\n" ; |
146 | throw 100003; |
147 | } |
148 | else |
149 | { |
150 | double g = GetGamma(shape+1.0, 1.0); |
151 | double w = GetUniform(); |
152 | return scale*g*std::pow(w, 1.0/shape); |
153 | } |
154 | } |
155 | |
156 | double GetChiSquare(double degreesOfFreedom) |
157 | { |
158 | // A chi squared distribution with n degrees of freedom |
159 | // is a gamma distribution with shape n/2 and scale 2. |
160 | return GetGamma(0.5 * degreesOfFreedom, 2.0); |
161 | } |
162 | |
163 | double GetInverseGamma(double shape, double scale) |
164 | { |
165 | // If X is gamma(shape, scale) then |
166 | // 1/Y is inverse gamma(shape, 1/scale) |
167 | return 1.0 / GetGamma(shape, 1.0 / scale); |
168 | } |
169 | |
170 | double GetWeibull(double shape, double scale) |
171 | { |
172 | if (shape <= 0.0 || scale <= 0.0) |
173 | { |
174 | std::cerr << __FILE__ << ":" << __LINE__ << " Shape and scale parameters must be positive. Recieved " << shape << " and " << scale; |
175 | throw 100004; |
176 | } |
177 | return scale * std::pow(-std::log(GetUniform()), 1.0 / shape); |
178 | } |
179 | |
180 | double GetCauchy(double median, double scale) |
181 | { |
182 | if (scale <= 0) |
183 | { |
184 | std::cerr << __FILE__ << ":" << __LINE__ << "Scale must be positive. Received " << scale << "\n" ; |
185 | throw 100005; |
186 | } |
187 | |
188 | double p = GetUniform(); |
189 | |
190 | // Apply inverse of the Cauchy distribution function to a uniform |
191 | return median + scale*std::tan(3.14159265359*(p - 0.5)); |
192 | } |
193 | |
194 | double GetStudentT(double degreesOfFreedom) |
195 | { |
196 | if (degreesOfFreedom <= 0) |
197 | { |
198 | std::cerr << __FILE__ << ":" << __LINE__ << " Degrees of freedom must be positive. Received " << degreesOfFreedom; |
199 | throw 100006; |
200 | } |
201 | |
202 | // See Seminumerical Algorithms by Knuth |
203 | double y1 = GetNormal(); |
204 | double y2 = GetChiSquare(degreesOfFreedom); |
205 | return y1 / std::sqrt(y2 / degreesOfFreedom); |
206 | } |
207 | |
208 | // The Laplace distribution is also known as the double exponential distribution. |
209 | double GetLaplace(double mean, double scale) |
210 | { |
211 | double u = GetUniform(); |
212 | return (u < 0.5) ? |
213 | mean + scale*std::log(2.0*u) : |
214 | mean - scale*std::log(2*(1-u)); |
215 | } |
216 | |
217 | double GetLogNormal(double mu, double sigma) |
218 | { |
219 | return std::exp(GetNormal(mu, sigma)); |
220 | } |
221 | |
222 | double GetBeta(double a, double b) |
223 | { |
224 | if (a <= 0.0 || b <= 0.0) |
225 | { |
226 | std::cerr << __FILE__ << ":" << __LINE__ << "Beta parameters must be positive. Received " << a << " and " << b << "\n" ; |
227 | throw 100007; |
228 | } |
229 | |
230 | // There are more efficient methods for generating beta samples. |
231 | // However such methods are a little more efficient and much more complicated. |
232 | // For an explanation of why the following method works, see |
233 | // http://www.johndcook.com/distribution_chart.html#gamma_beta |
234 | |
235 | double u = GetGamma(a, 1.0); |
236 | double v = GetGamma(b, 1.0); |
237 | return u / (u + v); |
238 | } |
239 | }; |
240 | |
241 | #endif |
242 | |