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
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
|
#include "low_discrepancy.hpp"
#include <vector>
#include <gsl/gsl_cdf.h>
#include "var_alea.hpp"
double frac_part(double x);
//fonctions de compositions de monte_carlo.hpp
template <class _Result>
struct generator
{
typedef _Result result_type;
};
template <typename Fct, typename VA>
struct compose_t : public generator< typename Fct::result_type >
{
compose_t(Fct f, VA X) : f(f), X(X) {};
typename Fct::result_type operator()() {
return f(X());
};
private:
Fct f; VA X;
};
template <typename Fct, typename VA>
inline compose_t<Fct, VA>
compose(Fct f, VA X) {
return compose_t<Fct, VA>(f, X);
};
template <typename VA>
struct mean_t: public generator<double>
{
mean_t(int n, VA X): n(n), X(X) {};
double operator()(){
double sum = 0;
for(int i=0; i<n; i++){
sum+=X();
}
return sum/n;
}
private :
int n; VA X;
};
template <typename VA>
inline mean_t<VA>
mean(int n, VA X){
return mean_t<VA>(n, X);
};
//Les classes de monte-Carlo
template <typename L>
std::vector<double> monte_carlo(int n, L X, double p=0.05)
{
std::vector<double> result(3,0);
double x;
for (int j = 0; j < n; j++) {
x = X();
result[0] += x;
result[1] += x*x;
}
result[0] /= (double) n;
result[1] = (result[1] - n*result[0]*result[0])/(double)(n-1);
result[2] = gsl_cdf_gaussian_Pinv(1-p/2,1)*sqrt(result[1]/(double) n);
return result;
};
template <typename L, typename Fct>
std::vector<double> monte_carlo(int n, Fct f, L X, double p=0.05)
{
return monte_carlo(n, compose(f, X), p);
};
//Les classes de quasi-mean
template <typename LDS>
struct quasi_gaussian : public var_alea<std::vector<double> >
{
quasi_gaussian (int d, double mean = 0, double std =1)
: d(d), mean(mean), std(std), s(d), U(0,1), seed(d) {
for(int i=0; i<d; i++) {seed[i]=U();}
};
std::vector<double> operator ()(){
std::vector<double> result = s();
for(int i=0;i<d; i++){
result[i] = mean + gsl_cdf_gaussian_Pinv(frac_part(result[i]+seed[i]), std);
}
return value = result;
}
private:
int d;
double mean, std;
LDS s;
uniform U;
std::vector<double> seed;
};
struct gaussian_d : public var_alea<std::vector<double> >
{
gaussian_d(int d, double mean = 0, double std =1)
: d(d), mean(mean), std(std), G(mean, std) {};
std::vector<double> operator ()(){
std::vector<double> result(d);
for(int i=0;i<d; i++){
result[i] = G();
}
return value = result;
}
private:
int d;
double mean, std;
gaussian G;
};
template <typename Fct, typename LDS>
struct quasi_mean : public generator<typename Fct::result_type>
{
quasi_mean(int n, int d, Fct f) : n(n), d(d), f(f), G(d) {};
typename Fct::result_type operator()() {
double sum =0;
for(int i=0; i<n; i++){
sum += f(G());
}
return sum/n;
};
private:
int n, d;
Fct f;
quasi_gaussian<LDS> G;
};
|
