diff options
Diffstat (limited to 'stratified_sampling.hpp')
| -rw-r--r-- | stratified_sampling.hpp | 167 |
1 files changed, 0 insertions, 167 deletions
diff --git a/stratified_sampling.hpp b/stratified_sampling.hpp deleted file mode 100644 index 81080fc..0000000 --- a/stratified_sampling.hpp +++ /dev/null @@ -1,167 +0,0 @@ -#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 sigma2=1, int seed=0) - :a(a), b(b), mean(mean), sigma2(sigma2), 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),sigma2(other.sigma2){ - gen = gsl_rng_clone(other.gen); - }; - - double operator()() { - pair<double, double> p = rtnorm(gen, a, b, mean, sigma2); - return value = p.first; - }; - - ~gaussian_truncated() { gsl_rng_free(gen); } -private: - double a, b, mean, sigma2; - int seed; - gsl_rng *gen; -}; - -template <typename L> -struct stratified_sampling { - stratified_sampling(vector<double> p, vector<L> X) - :p(p), X(X), mean(p.size(), 0), sigma2(p.size(), 0), I(p.size()){}; - void draw(int N); - vector<double> get_mean() const; - vector<double> get_var() const; - void print_mean() const; - void print_sigma() const; - pair<double,double> estimator() const; -private: - void update(int N); - vector<double> p; - vector<L> X; - vector<int> M; - vector<int> cumM; - vector<double> mean; - vector<double> sigma2; - const int I; -}; - -//actualisation du nombre de tirages à faire par strates -template <typename L> -void stratified_sampling<L>::update(int Nk) { - 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 { - //On remplit un vecter des écarts types à parti de notre vecteur de variance - std::vector<double> sigma(p.size(),0); - for (int i=0; i < I; i++) { - sigma[i]=sqrt(sigma2[i]); - } - 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-1; i++){ - M[i] += floor(current+m[i]) - floor(current); - current += m[i]; - } - M[I-1]+=Nk-I-floor(current); -} - -template <typename L> -void stratified_sampling<L>::draw(int N) { - update(N); - 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+X[i](); - s=s+X[i].current()*X[i].current(); - } - oldmean=mean[i]; - mean[i]=(mean[i]*cumM[i]+m)/(cumM[i]+M[i]); - sigma2[i]=((sigma2[i]+oldmean*oldmean)*cumM[i] + s)/(cumM[i]+M[i]) - mean[i]*mean[i]; - } -}; - -template <typename L> -vector<double> stratified_sampling<L>::get_mean() const { - return mean; -}; - -template <typename L> -vector<double> stratified_sampling<L>::get_var() const { - return sigma2; -}; - -template <typename L> -void stratified_sampling<L>::print_mean() const { - cout<<"les espérances :"<<endl; - for(int i=0;i<I;i++){ - cout<<mean[i]<<"\t"; - } - cout<<endl; -}; - -template <typename L> -void stratified_sampling<L>::print_sigma() const { - cout<<"les écarts types :"<<endl; - for(int i=0;i<I;i++){ - cout<<sqrt(sigma2[i])<<"\t"; - } - cout<<endl; -}; - -template <typename L> -pair<double,double> stratified_sampling<L>::estimator() const { - double est_mean = 0; - double est_std = 0; - for (int i=0; i<I; i++) { - est_mean += mean[i]*p[i]; - est_std += sqrt(sigma2[i])*p[i]; - } - return {est_mean, est_std}; -} |
