root / src / sim / distrib.cc @ e1750c09
History  View  Annotate  Download (9.58 KB)
1 
//==========================================================================


2 
// DISTRIB.CC

3 
//

4 
// OMNeT++/OMNEST

5 
// Discrete System Simulation in C++

6 
//

7 
// Random variate generation

8 
//

9 
// @author Werner Sandmann (ws), Kay Michael Masslow (kmm), Kyeong Soo (Joseph) Kim (jk)

10 
//

11 
// @date 11/26/2002 "doxygenification" (kmm)

12 
// @date 11/20/2002 some final comments (ws)

13 
// @date 10/22/2002 implemented various discrete distributions (kmm)

14 
//

15 
//==========================================================================

16  
17 
#include <float.h> 
18 
#include <math.h> 
19 
#include "distrib.h" 
20 
#include "globals.h" 
21 
#include "cmathfunction.h" 
22 
#include "cexception.h" 
23  
24 
//

25 
// Linux  <math.h> supplies a PI, MS does not...

26 
//

27 
#ifndef M_PI

28 
#define M_PI 3.1415926535897932384626433832795 
29 
#endif

30  
31 
#ifndef M_E

32 
#define M_E 2.7182818284590452353602874713527 
33 
#endif

34  
35 
NAMESPACE_BEGIN 
36  
37  
38 
//

39 
//

40 
// C O N T I N U O U S

41 
//

42 
//

43  
44 
double uniform(double a, double b, int rng) 
45 
{ 
46 
return a + genk_dblrand(rng) * (ba);

47 
} 
48  
49 
double exponential(double p, int rng) 
50 
{ 
51 
return p * log(1.0  genk_dblrand(rng)); 
52 
} 
53  
54 
double unit_normal(int rng) 
55 
{ 
56 
double U = 1.0  genk_dblrand(rng); 
57 
double V = 1.0  genk_dblrand(rng); 
58 
return sqrt(2.0*log(U)) * cos(PI*2*V); 
59 
} 
60  
61 
double normal(double m, double d, int rng) 
62 
{ 
63 
double U = 1.0  genk_dblrand(rng); 
64 
double V = 1.0  genk_dblrand(rng); 
65 
return m + d * sqrt(2.0*log(U)) * cos(PI*2*V); 
66 
} 
67  
68 
double truncnormal(double m, double d, int rng) 
69 
{ 
70 
double res;

71 
do {

72 
res = normal(m,d,rng); 
73 
} while(res<0); 
74  
75 
return res;

76 
} 
77  
78 
/*

79 
* internal, for alpha<1. THIS IMPLEMENTATION SEEMS TO BE BOGUS, we use

80 
* gamma_MarsagliaTransf() instead.

81 
*/

82 
/*

83 
static double gamma_AhrensDieter74(double alpha, int rng)

84 
{

85 
ASSERT(alpha<1);

86 

87 
double b = (M_E + alpha) / M_E;

88 

89 
double Y;

90 
for (;;)

91 
{

92 
double U1, U2, P;

93 

94 
// step 1

95 
U1 = normal(0, 1, rng);

96 
P = b * U1;

97 
if (P > 1)

98 
{

99 
// step 3

100 
Y = log((b  P) / alpha);

101 
U2 = normal(0, 1, rng);

102 
if (U2 <= pow(Y, alpha  1.0))

103 
break; // accept Y

104 
}

105 
else

106 
{

107 
// step 2

108 
Y = pow(P, (1 / alpha));

109 
U2 = normal(0, 1, rng);

110 
if (U2 <= exp(Y))

111 
break; // accept Y

112 
}

113 
}

114 
return Y;

115 
}

116 
*/

117  
118 
/*

119 
* internal, for alpha>1. Currently unused, we use gamma_Marsaglia2000() instead.

120 
*/

121 
/*

122 
static double gamma_ChengFeast79(double alpha, int rng)

123 
{

124 
ASSERT(alpha>1);

125 

126 
double a = 1.0 / sqrt(2.0 * alpha  1);

127 
double b = alpha  log(4.0);

128 
double q = alpha + 1.0 / a;

129 
double theta = 4.5;

130 
double d = 1 + log(theta);

131 

132 
double Y;

133 
for (;;)

134 
{

135 
double U1, U2, V, Z, W;

136 

137 
// step 1

138 
U1 = genk_dblrand(rng);

139 
U2 = genk_dblrand(rng);

140 

141 
// step 2

142 
V = a * log(U1 / (1.0  U1));

143 
Y = alpha * exp(V);

144 
Z = U1 * U1 * U2;

145 
W = b + q * V  Y;

146 

147 
// step 3

148 
if (W + d  theta * Z >= .0)

149 
break; // accept Y

150 

151 
// step 4

152 
if (W >= log(Z))

153 
break; // accept Y

154 
}

155 
return Y;

156 
}

157 
*/

158  
159 
/*

160 
* internal, for alpha>1.

161 
*

162 
* From: "A Simple Method for Generating Gamma Variables", George Marsaglia and

163 
* Wai Wan Tsang, ACM Transactions on Mathematical Software, Vol. 26, No. 3,

164 
* September 2000. Available online.

165 
*/

166 
static double gamma_Marsaglia2000(double a, int rng) 
167 
{ 
168 
ASSERT(a>1);

169  
170 
double d,c,x,v,u;

171 
d = a  1.0/3.0; 
172 
c = 1.0/sqrt(9.0*d); 
173 
for(;;)

174 
{ 
175 
do {x = unit_normal(rng); v = 1.0 + c*x;} while (v<=0); 
176 
v = v*v*v; u = genk_dblrand(rng); 
177 
if (u < 1.0  0.0331*(x*x)*(x*x)) 
178 
return d*v;

179 
if (log(u)<0.5*x*x + d*(1.0v+log(v))) 
180 
return d*v;

181 
} 
182 
} 
183  
184 
/*

185 
* internal, for alpha<1.

186 
*

187 
* We can derive the alpha<1 case from the alpha>1 case. See "Note" at the

188 
* end of Section 6 of the Marsaglia2000 paper (see above).

189 
*/

190 
static double gamma_MarsagliaTransf(double alpha, int rng) 
191 
{ 
192 
ASSERT(alpha<1);

193  
194 
//return gamma_ChengFeast79(1+alpha,rng) * pow(genk_dblrand(rng), 1/alpha);

195 
return gamma_Marsaglia2000(1+alpha,rng) * pow(genk_dblrand(rng), 1/alpha); // faster 
196 
} 
197  
198 
double gamma_d(double alpha, double theta, int rng) 
199 
{ 
200 
if (alpha<=0  theta<=0) 
201 
throw cRuntimeError("gamma(): alpha and theta params must be positive " 
202 
"(alpha=%g, theta=%g)", alpha, theta);

203  
204 
if (fabs(alpha  1.0) <= DBL_EPSILON) 
205 
{ 
206 
return exponential(theta, rng);

207 
} 
208 
else if (alpha < 1.0) 
209 
{ 
210 
//return theta * gamma_AhrensDieter74(alpha, rng); // implementation is bogus, see above

211 
return theta * gamma_MarsagliaTransf(alpha, rng);

212 
} 
213 
else // if (alpha > 1.0) 
214 
{ 
215 
//return theta * gamma_ChengFeast79(alpha, rng);

216 
return theta * gamma_Marsaglia2000(alpha, rng); // faster 
217 
} 
218 
} 
219  
220  
221 
double beta(double alpha1, double alpha2, int rng) 
222 
{ 
223 
if (alpha1<=0  alpha2<=0) 
224 
throw cRuntimeError("beta(): alpha1 and alpha2 parameters must be positive " 
225 
"(alpha1=%g, alpha2=%g)", alpha1, alpha2);

226  
227 
double Y1 = gamma_d(alpha1, 1.0, rng); 
228 
double Y2 = gamma_d(alpha2, 1.0, rng); 
229  
230 
return Y1 / (Y1 + Y2);

231 
} 
232  
233  
234 
double erlang_k(unsigned int k, double m, int rng) 
235 
{ 
236 
double U = 1.0; 
237 
for (unsigned int i = 0; i < k; i++) 
238 
U *= (1.0  genk_dblrand(rng)); 
239  
240 
return (m / (double) k) * log(U); 
241 
} 
242  
243  
244 
double chi_square(unsigned int k, int rng) 
245 
{ 
246 
if (!(k % 2)) 
247 
return erlang_k(k >> 1, k, rng); 
248 
else

249 
return gamma_d((double) k / 2.0, 2.0, rng); 
250 
} 
251  
252  
253 
double student_t(unsigned int i, int rng) 
254 
{ 
255 
double Z = normal(0, 1, rng); 
256 
double W = sqrt(chi_square(i, rng) / (double) i); 
257 
return Z / W;

258 
} 
259  
260  
261 
double cauchy(double a, double b, int rng) 
262 
{ 
263 
if (b<=0) 
264 
throw cRuntimeError("cauchy(): parameters must be b>0 (a=%g, b=%g)", a, b); 
265  
266 
return a + b * tan(M_PI * genk_dblrand(rng));

267 
} 
268  
269  
270 
double triang(double a, double b, double c, int rng) 
271 
{ 
272 
if (b<a  c<b  a==c)

273 
throw cRuntimeError("triang(): parameters must be a<=b<=c, a<c (a=%g, b=%g, c=%g)", a, b, c); 
274  
275 
double U, beta, T;

276  
277 
U = genk_dblrand(rng); 
278 
beta = (b  a) / (c  a); 
279  
280 
if (U < beta)

281 
T = sqrt(beta * U); 
282 
else

283 
T = 1.0  sqrt((1.0  beta) * (1.0  U)); 
284  
285 
return a + (c  a) * T;

286 
} 
287  
288  
289 
// lognormal() is inline

290  
291  
292 
double weibull(double a, double b, int rng) 
293 
{ 
294 
if (a<=0  b<=0) 
295 
throw cRuntimeError("weibull(): a,b parameters must be positive (a=%g, b=%g)", a, b); 
296  
297 
return a * pow(log(1.0  genk_dblrand(rng)), 1.0 / b); 
298 
} 
299  
300  
301 
double pareto_shifted(double a, double b, double c, int rng) 
302 
{ 
303 
if (a==0) 
304 
throw cRuntimeError("pareto_shifted(): parameter a cannot be zero)"); 
305  
306 
double u_pow = pow(1.0  genk_dblrand(rng), 1.0 / a); 
307 
return b / u_pow  c;

308 
} 
309  
310  
311 
//

312 
//

313 
// D I S C R E T E

314 
//

315 
//

316  
317 
/*

318 
// helper function, needed for hypergeometric distribution

319 
static double _factorial(int n)

320 
{

321 
if (n<0)

322 
throw cRuntimeError("internal error: _factorial() called with n=%d",n);

323 

324 
double fact = 1.0;

325 
while (n>1)

326 
fact *= n;

327 
return fact;

328 
}

329 
*/

330  
331  
332 
int intuniform(int a, int b, int rng) 
333 
{ 
334 
return a + genk_intrand(rng, ba+1); 
335 
} 
336  
337  
338 
// bernoulli() is inline

339  
340  
341 
int binomial(int n, double p, int rng) 
342 
{ 
343 
int X = 0; 
344 
// sum up n bernoulli trials

345 
for (int i = 0; i < n; i++) 
346 
{ 
347 
double U = genk_dblrand(rng);

348 
if (p > U)

349 
X++; 
350 
} 
351 
return X;

352 
} 
353  
354  
355 
int geometric(double p, int rng) 
356 
{ 
357 
if (p<0  p>=1) 
358 
throw cRuntimeError("geometric(): parameter p=%g out of range [0,1)", p); 
359  
360 
double a = 1.0 / (log(1.0  p)); 
361 
return (int)floor(a * log(1.0  genk_dblrand(rng))); 
362 
} 
363  
364  
365 
int negbinomial(int n, double p, int rng) 
366 
{ 
367 
int X = 0; 
368 
for (int i = 0; i < n; i++) 
369 
{ 
370 
X += geometric(p, rng); 
371 
} 
372 
return X;

373 
} 
374  
375  
376 
/*

377 
* TBD: hypergeometric() doesn't work yet

378 
*

379 
int hypergeometric(int a, int b, int n, int rng)

380 
{

381 
if (a<0  b<0  n>a+b)

382 
throw cRuntimeError("hypergeometric(): params must be a>=0, b>=0, n=<a+b "

383 
"(a=%d, b=%d, n=%d)", a,b,n);

384 

385 
double U = genk_dblrand(rng);

386 
double alpha = _factorial(b) / (double) _factorial(a + b  n);

387 
double beta = _factorial(b  n) / (double) _factorial(a + b);

388 
double A = alpha / beta;

389 
double B = A;

390 
int X = 0;

391 

392 
while (U > A)

393 
{

394 
X++;

395 
B *= double(a  X) * double(n  X) / ((X + 1.0) * (b  n + X + 1.0));

396 
A += B;

397 
}

398 

399 
return X;

400 
}

401 
*/

402  
403 
int poisson(double lambda, int rng) 
404 
{ 
405 
int X;

406 
if (lambda > 30.0) 
407 
{ 
408 
double a = M_PI * sqrt(lambda / 3.0); 
409 
double b = a / lambda;

410 
double c = 0.767  3.36 / lambda; 
411 
double d = log(c)  log(b)  lambda;

412 
double U, V, Y;

413  
414 
do

415 
{ 
416 
do

417 
{ 
418 
U = genk_dblrand(rng); 
419 
Y = (a  log((1.0  U) / U)) / b; 
420 
} 
421 
while (Y <= 0.5); 
422  
423 
X = (int)floor(Y + 0.5); 
424 
V = genk_dblrand(rng); 
425 
} 
426 
while (a  b * Y + log(V / 1.0 + pow(exp(a  b * Y), 2.0)) > d + X * log(lambda)  log(double(X))); 
427 
} 
428 
else

429 
{ 
430 
double a = exp(lambda);

431 
double p = 1.0; 
432 
X = 1;

433  
434 
while (p > a)

435 
{ 
436 
p *= genk_dblrand(rng); 
437 
X++; 
438 
} 
439 
} 
440 
return X;

441 
} 
442  
443  
444 
NAMESPACE_END 