From 045fc28b3432e07f3d9702150a56f21e7650f12a Mon Sep 17 00:00:00 2001
From: Bertrand <bertrand.horel@gmail.com>
Date: Sun, 21 Feb 2016 14:26:26 +0000
Subject: modif du Makefile et de var_alea et autres

---
 src/main.cpp                | 67 ++++++++++++++++++++++++---------------------
 src/opti.cpp                |  6 ++--
 src/stratified_sampling.cpp |  1 +
 src/stratified_sampling.hpp | 44 +++++++++++++++++++++--------
 src/var_alea.cpp            |  5 ++++
 src/var_alea.hpp            |  7 ++---
 6 files changed, 81 insertions(+), 49 deletions(-)
 create mode 100644 src/var_alea.cpp

(limited to 'src')

diff --git a/src/main.cpp b/src/main.cpp
index f81f1dd..f6e998a 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -17,37 +17,42 @@ vector<double> quantile_norm(int n, double sigma){
     return q;
 }
 
+void exemple1() {
+gsl_rng_env_setup();
+vector<double> q = quantile_norm(10, 1);
+vector<double> p(10, 0.1);
+vector<gaussian_truncated> rvar;
+rvar.push_back(gaussian_truncated(GSL_NEGINF, q[0],0,1,0));
+for (int i=1; i<10; i++){
+    rvar.push_back(gaussian_truncated(q[i-1], q[i],0,1,i));
+};
+stratified_sampling<gaussian_truncated> S(p,rvar);
+S.draw(100);
+double x = 1.64*S.estimator().second;
+    cout<<"l'estimateur de la moyenne est :"<<S.estimator().first<<endl;
+    cout<<"Son intervalle de confiance à 95% est :"<<"["<<S.estimator().first-(x/10)<<" ,"<<S.estimator().first+(x/10)<<"]"<<endl;
+S.draw(1000);
+x = 1.64*S.estimator().second;
+    cout<<"l'estimateur de la moyenne est :"<<S.estimator().first<<endl;
+    cout<<"Son intervalle de confiance à 95% est :"<<"["<<S.estimator().first-(x/sqrt(1100))<<" ,"<<S.estimator().first+(x/sqrt(1100))<<"]"<<endl;
+    S.draw(10000);
+    x = 1.64*S.estimator().second;
+    cout<<"l'estimateur de la moyenne est :"<<S.estimator().first<<endl;
+    cout<<"Son intervalle de confiance à 95% est :"<<"["<<S.estimator().first-(x/sqrt(11100))<<" ,"<<S.estimator().first+(x/sqrt(11100))<<"]"<<endl;
+    
+};
+
+
+
 int main()
-{
-    //--- GSL random init ---
-    gsl_rng_env_setup();
-    vector<double> q = quantile_norm(10, 1);
-    vector<double> p(10, 0.1);
-    vector<gaussian_truncated> rvar;
-    rvar.push_back(gaussian_truncated(GSL_NEGINF, q[0],0,1,0));
-    for (int i=1; i<10; i++){
-        rvar.push_back(gaussian_truncated(q[i-1], q[i],0,1,i));
+{   
+    std::vector<double> u {sqrt(0.2), sqrt(0.2), sqrt(0.2), sqrt(0.2), sqrt(0.2)};
+    multi_gaussian_truncated G(0, 2, u);
+    std::vector<double> r(5);
+    r = G();
+    for (int i=0; i<5; i++){
+        std::cout<<r[i]<<std::endl;
     }
-    stratified_sampling<gaussian_truncated> S(p,rvar);
-    //S.update(100);
-    S.draw(100);
-    double x = 1.64*S.estimator().second;
-		cout<<"l'estimateur de la moyenne est :"<<S.estimator().first<<endl;
-		cout<<"Son intervalle de confiance à 95% est :"<<"["<<S.estimator().first-(x/10)<<" ,"<<S.estimator().first+(x/10)<<"]"<<endl;
-    //S.update(1000);
-    S.draw(1000);
-		x = 1.64*S.estimator().second;
-		cout<<"l'estimateur de la moyenne est :"<<S.estimator().first<<endl;
-		cout<<"Son intervalle de confiance à 95% est :"<<"["<<S.estimator().first-(x/sqrt(1100))<<" ,"<<S.estimator().first+(x/sqrt(1100))<<"]"<<endl;
-		//S.update(10000);
-		S.draw(10000);
-		x = 1.64*S.estimator().second;
-		cout<<"l'estimateur de la moyenne est :"<<S.estimator().first<<endl;
-		cout<<"Son intervalle de confiance à 95% est :"<<"["<<S.estimator().first-(x/sqrt(11100))<<" ,"<<S.estimator().first+(x/sqrt(11100))<<"]"<<endl;
-
-    //~ S.draw();
-    //~ for(int i=0;i<10;i++){
-        //~ cout<<S.get_mean()[i]<<endl;
-    //~ }
-    return 0;
+        
+    return 0; 
 }
diff --git a/src/opti.cpp b/src/opti.cpp
index 1b5df1d..702e8c2 100644
--- a/src/opti.cpp
+++ b/src/opti.cpp
@@ -41,17 +41,17 @@ int main() {
     opt.set_xtol_rel(1e-4);
 
     std::vector<double> x(16,0);
-    std::vector<double> g();
+    std::vector<double> g(0);
 
     
-    std::cout<<"valeur au début"<<f(x, g, &params)<<std::endl;
+    std::cout<<"valeur au début : "<<f(x, g, &params)<<std::endl;
 
     double maxf;
     nlopt::result result = opt.optimize(x, maxf);
     for(int i=0; i<16; i++){
         std::cout<<x[i]<<std::endl;
     }
-    std::cout<<"valeur à la fin"<<maxf<<std::endl;
+    std::cout<<"valeur à la fin : "<<maxf<<std::endl;
 
     return 0;
 }
diff --git a/src/stratified_sampling.cpp b/src/stratified_sampling.cpp
index c1c8dc1..156751f 100644
--- a/src/stratified_sampling.cpp
+++ b/src/stratified_sampling.cpp
@@ -11,3 +11,4 @@ std::pair<double, double> mean_var( std::vector<double> r){
     p.second -= p.first * p.first;
     return p;
 }
+
diff --git a/src/stratified_sampling.hpp b/src/stratified_sampling.hpp
index 81080fc..cebee8d 100644
--- a/src/stratified_sampling.hpp
+++ b/src/stratified_sampling.hpp
@@ -3,20 +3,11 @@
 #include <iostream>
 #include <gsl/gsl_rng.h>
 #include "rtnorm.hpp"
+#include "var_alea.hpp"
+#include <gsl/gsl_cdf.h>
 
 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;
 
@@ -45,6 +36,37 @@ private:
     gsl_rng *gen;
 };
 
+struct multi_gaussian_truncated : public var_alea<std::vector<double> >
+{
+    multi_gaussian_truncated(double a, double b, const std::vector<double> u)
+        :a(a), b(b), V(gsl_cdf_ugaussian_P(a), gsl_cdf_ugaussian_P(b)), G(0,1), u(u), d(u.size()) {};
+    
+    std::vector<double> operator()() {
+        double v = V();
+        double Z = gsl_cdf_gaussian_Pinv(v,1);
+        std::vector<double> Y(d);
+        double scal = 0;
+        for(int i=0; i<d; i++){
+            Y[i] = G();
+        }
+        std::vector<double> X(d);
+        for(int j=0; j<d; j++){
+            scal += Y[j]*u[j];
+        }
+        for(int i=0; i<d; i++){
+            X[i] = u[i]*Z + Y[i] - u[i]*scal;
+        }
+        return X;
+    }
+    
+    private:
+        double a, b;        
+        uniform V;
+        gaussian G;
+        std::vector<double> u;
+        int d;
+    };
+
 template <typename L>
 struct stratified_sampling {
     stratified_sampling(vector<double> p, vector<L> X)
diff --git a/src/var_alea.cpp b/src/var_alea.cpp
new file mode 100644
index 0000000..958d8fc
--- /dev/null
+++ b/src/var_alea.cpp
@@ -0,0 +1,5 @@
+#include "var_alea.hpp"
+
+void init_alea(unsigned seed) {
+	init_genrand(seed);
+};
diff --git a/src/var_alea.hpp b/src/var_alea.hpp
index 2d1e8a6..cba365a 100644
--- a/src/var_alea.hpp
+++ b/src/var_alea.hpp
@@ -7,9 +7,8 @@ extern "C" {
 #include <ctime>
 
 
-void init_alea(unsigned seed = static_cast<unsigned>(std::time(0))) {
-	init_genrand(seed);
-};
+void init_alea(unsigned seed = static_cast<unsigned>(std::time(0)));
+
 
 template <typename T>
 struct var_alea {
@@ -49,7 +48,7 @@ struct expo : public var_alea<double>
 struct gaussian : public var_alea<double>
 {
     gaussian(double mean = 0, double std = 1) 
-        : mean(mean), std(std), flag(true), unif(-1,1) {};
+        : mean(mean), std(std), unif(-1,1), flag(true) {};
     double operator()() {
         flag = !flag;
         if (!flag) {
-- 
cgit v1.2.3-70-g09d2