From a86b1d918a4c31f74024aa6b42fa0021919585de Mon Sep 17 00:00:00 2001
From: Bertrand <bertrand.horel@gmail.com>
Date: Mon, 8 Feb 2016 14:40:20 +0000
Subject: débuggage et seeding
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 stratified_sampling.hpp | 126 ++++++++++++++++++++++++++++++++++++++----------
 1 file changed, 101 insertions(+), 25 deletions(-)

(limited to 'stratified_sampling.hpp')

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;
 };
-- 
cgit v1.2.3-70-g09d2