Statistics
| Branch: | Revision:

root / src / sim / distrib.cc @ e1750c09

History | View | Annotate | Download (9.58 KB)

1 01873262 Georg Kunz
//==========================================================================
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) * (b-a);
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.0-v+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, b-a+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