1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
|
#ifndef VAR_ALEA_HPP
#define VAR_ALEA_HPP
#include <cmath>
extern "C" {
#include "mt19937.h"
}
#include <ctime>
void init_alea(unsigned seed = static_cast<unsigned>(std::time(0))) {
init_genrand(seed);
};
template <typename T>
struct var_alea {
typedef T result_type;
var_alea() : value(0) {};
var_alea(T value) : value(value) {};
virtual ~var_alea() {};
virtual T operator()() = 0;
T current() const { return value; };
protected:
T value;
};
struct uniform : public var_alea<double>
{
uniform(double left = 0, double right = 1)
: left(left), size(right-left), genrand(genrand_real3) {};
double operator()() {
return value = left + size * genrand();
};
private:
double left, size;
double (*genrand)(void);
};
struct expo : public var_alea<double>
{
expo(double lambda) : inv_lambda(1./lambda), U(0,1) {};
double operator()() {
return value = - inv_lambda * log(U());
};
private:
double inv_lambda;
uniform U;
};
struct gaussian : public var_alea<double>
{
gaussian(double mean = 0, double std = 1)
: mean(mean), std(std), flag(true), unif(-1,1) {};
double operator()() {
flag = !flag;
if (!flag) {
do {
U = unif(); V = unif();
R2 = U*U + V*V;
} while (R2 > 1);
rac = sqrt(-2 * log(R2) / R2);
return value = mean + std * U * rac;
} else
return value = mean + std * V * rac;
};
private:
double mean, std, U, V, R2, rac;
uniform unif;
bool flag;
};
struct chi_deux : public var_alea<double>
{
chi_deux(int n) : n(n), G(0,1) {};
double operator()() {
value = 0;
for (int j = 0; j < n; j++) value += G()*G.current();
return value;
};
private:
int n;
gaussian G;
};
struct inverse_gaussian : public var_alea<double>
{
inverse_gaussian(double lambda, double mu)
: lambda(lambda), mu(mu), Y(1), U(0,1) {};
double operator()() {
double Z = mu + 0.5*mu*mu/lambda*Y();
double rac = sqrt(Z*Z - mu*mu);
return value = (U() < mu/(mu+Z+rac)) ? Z+rac : Z-rac;
};
private:
double lambda, mu;
chi_deux Y;
uniform U;
};
struct normal_inverse_gaussian : public var_alea<double>
{
normal_inverse_gaussian(double alpha, double beta,
double mu, double delta)
: alpha(alpha), beta(beta), mu(mu), delta(delta), G(0,1),
Y(delta/sqrt(alpha*alpha-beta*beta), delta*delta) {};
double operator()() {
double y_ = Y();
return value = mu + beta*y_ * sqrt(y_) * G();
};
private:
double alpha, beta, mu, delta;
gaussian G;
inverse_gaussian Y;
};
#endif
|
