diff options
| -rw-r--r-- | main.cpp | 28 | ||||
| -rw-r--r-- | stratified_sampling.cpp | 38 | ||||
| -rw-r--r-- | stratified_sampling.hpp | 126 |
3 files changed, 119 insertions, 73 deletions
@@ -1,16 +1,16 @@ #include <iostream> #include <gsl/gsl_rng.h> #include <vector> -//#include "rtnorm.hpp" #include <gsl/gsl_cdf.h> #include <gsl/gsl_math.h> #include "stratified_sampling.hpp" #include <cmath> #include <algorithm> +using namespace std; //--génération quantiles-- -std::vector<double> quantile_norm(int n, double sigma){ - std::vector<double> q(n); +vector<double> quantile_norm(int n, double sigma){ + vector<double> q(n); for (int i=0; i<n; i++) { q[i] = gsl_cdf_gaussian_Pinv ((double)(i+1)/n, sigma); } @@ -21,15 +21,23 @@ int main() { //--- GSL random init --- gsl_rng_env_setup(); - std::vector<double> q; - q = quantile_norm(10, 1); + vector<double> q = quantile_norm(10, 1); vector<double> p(10, 0.1); - vector<gaussian_truncated> gen; - gen.push_back(gaussian_truncated(GSL_NEGINF,q[0])); + vector<gaussian_truncated> rvar; + rvar.push_back(gaussian_truncated(GSL_NEGINF, q[0],0,1,0)); for (int i=1; i<10; i++){ - gen.push_back(gaussian_truncated(q[i-1],q[i])); + rvar.push_back(gaussian_truncated(q[i-1], q[i],0,1,i)); + } + stratified_sampling<gaussian_truncated> S(p,rvar); + S.update(100); + S.draw(); + for(int i=0;i<10;i++){ + cout<<S.get_mean()[i]<<endl; + } + S.update(500); + S.draw(); + for(int i=0;i<10;i++){ + cout<<S.get_mean()[i]<<endl; } - stratified_sampling<gaussian_truncated> S(p, gen); - //S.update(1000); return 0; } diff --git a/stratified_sampling.cpp b/stratified_sampling.cpp index b959776..c1c8dc1 100644 --- a/stratified_sampling.cpp +++ b/stratified_sampling.cpp @@ -1,6 +1,4 @@ #include "stratified_sampling.hpp" -#include <algorithm> -#include <iostream> std::pair<double, double> mean_var( std::vector<double> r){ std::pair<double, double> p; @@ -13,39 +11,3 @@ std::pair<double, double> mean_var( std::vector<double> r){ p.second -= p.first * p.first; return p; } - -//actualisation du nombre de tirages à faire par strates -template <typename Gen> -void stratified_sampling<Gen>::update(int Nk) { - int I = p.size(); - //reinitialistation du vecteur M du nombre de tirages par strates - if (M.empty()) { - M.resize(I,1); - } - else { - for(int i=0; i<I; i++){ - M[i]=1; - } - } - - std::vector<double> m(I, 0); //le vecteur des m_i idéals - - if (sigma.empty()) { - 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]; - } -} diff --git a/stratified_sampling.hpp b/stratified_sampling.hpp index b63cd52..0b0b4ef 100644 --- a/stratified_sampling.hpp +++ b/stratified_sampling.hpp @@ -1,4 +1,6 @@ #include <vector> +#include <algorithm> +#include <iostream> #include <gsl/gsl_rng.h> #include "rtnorm.hpp" @@ -9,43 +11,117 @@ struct var_alea { typedef T result_type; var_alea() : value(0) {}; var_alea(T value) : value(value) {}; - virtual ~var_alea() {}; + virtual ~var_alea() {}; virtual T operator()() = 0; T current() const { return value; }; - protected: - T 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) - :a(a), b(b), mean(mean), sigma(sigma) { - const gsl_rng_type* type = gsl_rng_default; - gen = gsl_rng_alloc(type); - }; - double operator()() { - pair<double, double> p = rtnorm(gen, a, b, mean, sigma); - return value = p.first; - }; - ~gaussian_truncated() { gsl_rng_free(gen); } + 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: - double mean, sigma, a, b; - gsl_rng *gen; + 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) {}; - void update(int N); - //vector<double> get_mean(); - //double estimator(); + 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<double> sigma; - vector<Gen> gen; - + 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; }; |