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
|
#include <vector>
#include <algorithm>
#include <iostream>
#include <gsl/gsl_rng.h>
#include "rtnorm.hpp"
using namespace std;
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;
};
typedef var_alea<double> var_alea_real;
struct gaussian_truncated : public var_alea_real
{
gaussian_truncated(double a, double b, double mean=0, double sigma=1, int seed=0)
:a(a), b(b), mean(mean), sigma(sigma), seed(seed) {
const gsl_rng_type* type = gsl_rng_default;
gen = gsl_rng_alloc(type);
gsl_rng_set(gen, seed);
};
gaussian_truncated(gaussian_truncated const &other)
:a(other.a),b(other.b), mean(other.mean),sigma(other.sigma){
gen = gsl_rng_clone(other.gen);
};
double operator()() {
pair<double, double> p = rtnorm(gen, a, b, mean, sigma);
return value = p.first;
};
~gaussian_truncated() { gsl_rng_free(gen); }
private:
int seed;
double mean, sigma, a, b;
gsl_rng *gen;
};
template <typename Gen>
struct stratified_sampling {
stratified_sampling(vector<double> p, vector<Gen> gen)
:p(p), gen(gen), mean(p.size(), 0), sigma(p.size(), 0){};
void update(int N);
void draw();
vector<double> get_mean();
//double estimator();
private:
vector<double> p;
vector<int> M;
vector<int> cumM;
vector<double> mean;
vector<double> sigma;
vector<Gen> gen;
};
//actualisation du nombre de tirages à faire par strates
template <typename Gen>
void stratified_sampling<Gen>::update(int Nk) {
int I = p.size();
bool first_step = M.empty();
//reinitialistation du vecteur M du nombre de tirages par strates
if (first_step) {
M.resize(I,1);
cumM.resize(I,0);
}
else {
for(int i=0; i<I; i++){
cumM[i]=cumM[i]+M[i];
M[i]=1;
}
}
std::vector<double> m(I, 0); //le vecteur des m_i idéals
if (first_step) {
for (int i=0; i<I; i++) {
m[i] = (Nk-I)*p[i];
}
}
else {
double scal = std::inner_product(p.begin(), p.end(), sigma.begin(), (double) 0);
for (int i=0; i < I; i++) {
m[i] = (Nk-I)*p[i]*sigma[i]/scal;
std::cout<<m[i]<<std::endl;
}
}
M[0]+=floor(m[0]);
double current = m[0];
for (int i=1; i<I; i++){
M[i] += floor(current+m[i]) - floor(current);
current += m[i];
cout<<M[i]<<"\t";
}
cout<<endl;
}
template <typename Gen>
void stratified_sampling<Gen>::draw() {
int I = p.size();
double m, s, oldmean;
for(int i=0;i<I;i++){
m=0;
s=0;
for(int j=0;j<M[i];j++){
m=m+gen[i]();
s=s+gen[i].current()*gen[i].current();
}
oldmean=mean[i];
mean[i]=(mean[i]*cumM[i]+m)/(cumM[i]+M[i]);
sigma[i]=((sigma[i]+oldmean*oldmean)*cumM[i] + s)/(cumM[i]+M[i]) - mean[i]*mean[i];
}
};
template <typename Gen>
vector<double> stratified_sampling<Gen>::get_mean() {
return mean;
};
|
