diff options
| -rw-r--r-- | Makefile | 4 | ||||
| -rw-r--r-- | main.cpp | 35 | ||||
| -rw-r--r-- | stratified_sampling.cpp | 130 | ||||
| -rw-r--r-- | stratified_sampling.hpp | 51 |
4 files changed, 123 insertions, 97 deletions
@@ -4,11 +4,11 @@ RM = rm CXXFLAGS=-std=c++11 -O3 GSL_FLAGS:=$(shell pkg-config --libs gsl) -stratified_sampling: rtnorm.o stratified_sampling.o +main: main.o rtnorm.o stratified_sampling.o $(CXX) $^ -o $@ $(GSL_FLAGS) test: test.o mt19937.o $(CXX) $^ -o $@ clean: - -$(RM) -f *.o test stratified_sampling + -$(RM) -f *.o test stratified_sampling main diff --git a/main.cpp b/main.cpp new file mode 100644 index 0000000..c1eb3ec --- /dev/null +++ b/main.cpp @@ -0,0 +1,35 @@ +#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> + +//--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; +} + +int main() +{ + //--- GSL random init --- + gsl_rng_env_setup(); + std::vector<double> q; + q = quantile_norm(10, 1); + vector<double> p(10, 0.1); + vector<gaussian_truncated> gen; + gen.push_back(gaussian_truncated(GSL_NEGINF,q[0])); + for (int i=1; i<10; i++){ + gen.push_back(gaussian_truncated(q[i-1],q[i])); + } + stratified_sampling<gaussian_truncated> S(p, gen); + //S.update(1000); + return 0; +} diff --git a/stratified_sampling.cpp b/stratified_sampling.cpp index 8f408cc..c03d447 100644 --- a/stratified_sampling.cpp +++ b/stratified_sampling.cpp @@ -1,31 +1,6 @@ -#include <iostream> -#include <gsl/gsl_rng.h> -#include <vector> -#include "rtnorm.hpp" -#include <gsl/gsl_cdf.h> -#include <gsl/gsl_math.h> - -#include <cmath> +#include "stratified_sampling.hpp" #include <algorithm> - -//--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; - } +#include <iostream> std::pair<double, double> mean_var( std::vector<double> r){ std::pair<double, double> p; @@ -39,75 +14,40 @@ std::pair<double, double> mean_var( std::vector<double> r){ return p; } //actualisation du nombre de tirages à faire par strates -std::vector<int> update_sampling (std::vector<double> p, - std::vector<double> sigma, int Nk) { - int I = p.size(); - std::vector<int> M(I, 1); // notre vecteur final à retourner - 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]; - } - return M; -} - - -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; +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<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); + + 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]; + } - return 0; } + + diff --git a/stratified_sampling.hpp b/stratified_sampling.hpp new file mode 100644 index 0000000..b63cd52 --- /dev/null +++ b/stratified_sampling.hpp @@ -0,0 +1,51 @@ +#include <vector> +#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) + :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); } +private: + 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(); +private: + vector<double> p; + vector<int> M; + vector<double> sigma; + vector<Gen> gen; + +}; |