From 22358d075a89d44dbe8ef901a0a42d2c3f8e285a Mon Sep 17 00:00:00 2001
From: Bertrand <bertrand.horel@gmail.com>
Date: Sat, 6 Feb 2016 19:32:06 +0000
Subject: encapsulage du stratified sampling

---
 stratified_sampling.cpp | 130 +++++++++++++-----------------------------------
 1 file changed, 35 insertions(+), 95 deletions(-)

(limited to 'stratified_sampling.cpp')

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