1 files changed, 33 insertions, 11 deletions
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)