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
|
#include <iostream>
#include <gsl/gsl_rng.h>
#include <vector>
#include "rtnorm.hpp"
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_math.h>
#define LOW_DISCREPANCY_HPP
#include <cmath>
#include <climits>
#include <list>
#include <algorithm>
#include <numeric>
//--génération quantiles--
std::vector<double> quantile_norm(int n, double sigma){
std::vector<double> q(n);
for (int i=0; i<n; i++) {
q[i] = gsl_cdf_gaussian_Pinv ((double)(i+1)/n, sigma);
}
return q;
}
//--tirage de normale tronquée entre les quantile de taille 1/n i et i+1--
double quantile_truncate_normal (int i, int n, double mu,
double sigma, gsl_rng *gen) {
std::vector<double> q;
q = quantile_norm(n, sigma);
std::pair<double, double> p;
p = rtnorm (gen, q[i], q[i+1], mu, sigma);
return p.first;
}
std::pair<double, double> mean_var( std::vector<double> r){
std::pair<double, double> p;
for(auto &x: r){
p.first += x;
p.second += x*x;
}
p.first /= r.size();
p.second /= r.size();
p.second -= p.first * p.first;
return p;
}
//actualisation du nombre de tirages à faire par strates
std::vector<int> update_sampling (std::vector<double> p,
std::vector<double> sigma, int n) {
std::vector<int> r; // notre vecteur final à retourner
std::vector<double> m; //le vecteur des Mi idéals
int s = p.size();
if (sigma.empty()) {
for (int i = 0; i<s ; i++) {
m.push_back(double(n)*p[i]);
}
}
else {
double scal = std::inner_product(p.begin(), p.end(), sigma.begin(), (double) 0);
for (int i=0; i<s; i++) {
double v = p[i]*sigma[i]/scal;
m.push_back(v*n);
std::cout<<m[i]<<std::endl;
}
}
r.push_back((int)(m[0]));
double v1 = 0;
for (int i=1; i<s-1; i++){
for (int j=0; j<i; j++) {
v1=v1+m[j];
}
r.push_back((int)v1 - (int)(v1-m[i]));
}
r.push_back(n - std::accumulate( r.begin(), r.end(), 0 ));
return r;
}
int main()
{
//--- GSL random init ---
gsl_rng_env_setup(); // Read variable environnement
const gsl_rng_type* type = gsl_rng_default; // Default algorithm 'twister'
gsl_rng *gen = gsl_rng_alloc (type); // Rand generator allocation
std::vector<double> q;
q = quantile_norm(10, 1);
std::pair<double, double> p;
std::pair<double, double> mv;
//number of classes
int I = 10;
//number of samples
int N = 10000;
std::vector<double> r(N);
double a;
for (int i=0; i<I; i++){
if(i==0){
a = GSL_NEGINF;
}else{
a = q[i-1];
}
for(int j=0; j<N; j++){
p = rtnorm (gen, a, q[i], 0, 1);
r[j] = p.first;
}
mv = mean_var(r);
//std::cout<<"mean :"<<mv.first<<" var :"<<mv.second<<std::endl;
}
std::vector<int> k;
std::vector<double> z = {(double)1/3,(double)1/3,(double)1/3};
std::vector<double> sigma = {0.1, 0.4, 0.3};
k = update_sampling(z, sigma, 10000);
for (int j=0; j<k.size(); j++){
std::cout<<k[j]<<std::endl;
}
gsl_rng_free(gen);
return 0;
}
|
