From d2b133901a65244934eb642ec8e20c797efaf650 Mon Sep 17 00:00:00 2001
From: Bertrand <bertrand.horel@gmail.com>
Date: Fri, 19 Feb 2016 15:03:51 +0000
Subject: nettoyage du dépôt
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 src/stratified_sampling.hpp | 167 ++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 167 insertions(+)
 create mode 100644 src/stratified_sampling.hpp

(limited to 'src/stratified_sampling.hpp')

diff --git a/src/stratified_sampling.hpp b/src/stratified_sampling.hpp
new file mode 100644
index 0000000..81080fc
--- /dev/null
+++ b/src/stratified_sampling.hpp
@@ -0,0 +1,167 @@
+#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};
+}
-- 
cgit v1.2.3-70-g09d2