aboutsummaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Makefile2
-rw-r--r--src/main.cpp67
-rw-r--r--src/opti.cpp6
-rw-r--r--src/stratified_sampling.cpp1
-rw-r--r--src/stratified_sampling.hpp44
-rw-r--r--src/var_alea.cpp5
-rw-r--r--src/var_alea.hpp7
7 files changed, 82 insertions, 50 deletions
diff --git a/Makefile b/Makefile
index 37e33e7..a0dbf06 100644
--- a/Makefile
+++ b/Makefile
@@ -5,7 +5,7 @@ CXXFLAGS=-std=c++11 -g -O3 -Wall
GSL_FLAGS:=$(shell pkg-config --libs gsl)
VPATH=src
-main: main.o rtnorm.o stratified_sampling.o
+main: main.o rtnorm.o stratified_sampling.o mt19937.o var_alea.o
$(CXX) $^ -o $@ $(GSL_FLAGS)
stratified_sampling.o: stratified_sampling.hpp
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) {
generated by cgit v1.2.3-70-g09d2 (git 2.52.0) at 2025-12-21 16:20:11 +0000