From ca2758864ee33fc7aebbbf65a87bc9074afc4a3b Mon Sep 17 00:00:00 2001
From: Guillaume Horel <guillaume.horel@serenitascapital.com>
Date: Mon, 1 Feb 2016 13:47:41 -0500
Subject: compute mean and variance for each class

---
 stratified_sampling.cpp | 41 ++++++++++++++++++++++++++++++++++-------
 1 file changed, 34 insertions(+), 7 deletions(-)

(limited to 'stratified_sampling.cpp')

diff --git a/stratified_sampling.cpp b/stratified_sampling.cpp
index 4d3de48..22bdbf5 100644
--- a/stratified_sampling.cpp
+++ b/stratified_sampling.cpp
@@ -3,7 +3,7 @@
 
 #include "rtnorm.hpp"
 #include <gsl/gsl_cdf.h>
-
+#include <gsl/gsl_math.h>
 
 //--génération quantiles--
 std::vector<double> quantile_norm(int n, double sigma){
@@ -23,8 +23,18 @@ double quantile_truncate_normal (int i, int n, double mu,
 	 p = rtnorm (gen, q[i], q[i+1], mu, sigma);
 	 return p.first;
  }
-  	 
 
+std::pair<double, double> mean_var( std::vector<double> r){
+    std::pair<double, double> p;
+    for(int i=0; i<r.size(); i++){
+        p.first += r[i];
+        p.second += r[i]*r[i];
+    }
+    p.first /= r.size();
+    p.second /= r.size();
+    p.second -= p.first * p.first;
+    return p;
+}
 
 int main()
 {
@@ -35,11 +45,28 @@ int main()
   std::vector<double> q;
   q = quantile_norm(10, 1);
   std::cout<<q[5]<<std::endl;
+
+  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;
-  int i;
-  for (i=0; i<20; i++){
-  a= quantile_truncate_normal(4, 10, 0, 1, gen);
-  std::cout<<a<<std::endl;}
+  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;
+  }
+  gsl_rng_free(gen);
   return 0;
 }
-
-- 
cgit v1.2.3-70-g09d2