aboutsummaryrefslogtreecommitdiffstats
path: root/stratified_sampling.hpp
diff options
context:
space:
mode:
Diffstat (limited to 'stratified_sampling.hpp')
-rw-r--r--stratified_sampling.hpp167
1 files changed, 0 insertions, 167 deletions
diff --git a/stratified_sampling.hpp b/stratified_sampling.hpp
deleted file mode 100644
index 81080fc..0000000
--- a/stratified_sampling.hpp
+++ /dev/null
@@ -1,167 +0,0 @@
-#include <vector>
-#include <algorithm>
-#include <iostream>
-#include <gsl/gsl_rng.h>
-#include "rtnorm.hpp"
-
-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;
-
-struct gaussian_truncated : public var_alea_real
-{
- gaussian_truncated(double a, double b, double mean=0, double sigma2=1, int seed=0)
- :a(a), b(b), mean(mean), sigma2(sigma2), seed(seed) {
- const gsl_rng_type* type = gsl_rng_default;
- gen = gsl_rng_alloc(type);
- gsl_rng_set(gen, seed);
- };
- gaussian_truncated(gaussian_truncated const &other)
- :a(other.a),b(other.b), mean(other.mean),sigma2(other.sigma2){
- gen = gsl_rng_clone(other.gen);
- };
-
- double operator()() {
- pair<double, double> p = rtnorm(gen, a, b, mean, sigma2);
- return value = p.first;
- };
-
- ~gaussian_truncated() { gsl_rng_free(gen); }
-private:
- double a, b, mean, sigma2;
- int seed;
- gsl_rng *gen;
-};
-
-template <typename L>
-struct stratified_sampling {
- stratified_sampling(vector<double> p, vector<L> X)
- :p(p), X(X), mean(p.size(), 0), sigma2(p.size(), 0), I(p.size()){};
- void draw(int N);
- vector<double> get_mean() const;
- vector<double> get_var() const;
- void print_mean() const;
- void print_sigma() const;
- pair<double,double> estimator() const;
-private:
- void update(int N);
- vector<double> p;
- vector<L> X;
- vector<int> M;
- vector<int> cumM;
- vector<double> mean;
- vector<double> sigma2;
- const int I;
-};
-
-//actualisation du nombre de tirages à faire par strates
-template <typename L>
-void stratified_sampling<L>::update(int Nk) {
- bool first_step = M.empty();
- //reinitialistation du vecteur M du nombre de tirages par strates
- if (first_step) {
- M.resize(I,1);
- cumM.resize(I,0);
- }
- else {
- for(int i=0; i<I; i++){
- cumM[i]=cumM[i]+M[i];
- M[i]=1;
- }
- }
-
- std::vector<double> m(I, 0); //le vecteur des m_i idéals
-
- if (first_step) {
- for (int i=0; i<I; i++) {
- m[i] = (Nk-I)*p[i];
- }
- }
- else {
- //On remplit un vecter des écarts types à parti de notre vecteur de variance
- std::vector<double> sigma(p.size(),0);
- for (int i=0; i < I; i++) {
- sigma[i]=sqrt(sigma2[i]);
- }
- 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-1; i++){
- M[i] += floor(current+m[i]) - floor(current);
- current += m[i];
- }
- M[I-1]+=Nk-I-floor(current);
-}
-
-template <typename L>
-void stratified_sampling<L>::draw(int N) {
- update(N);
- double m, s, oldmean;
- for(int i=0;i<I;i++){
- m=0;
- s=0;
- for(int j=0;j<M[i];j++){
- m=m+X[i]();
- s=s+X[i].current()*X[i].current();
- }
- oldmean=mean[i];
- mean[i]=(mean[i]*cumM[i]+m)/(cumM[i]+M[i]);
- sigma2[i]=((sigma2[i]+oldmean*oldmean)*cumM[i] + s)/(cumM[i]+M[i]) - mean[i]*mean[i];
- }
-};
-
-template <typename L>
-vector<double> stratified_sampling<L>::get_mean() const {
- return mean;
-};
-
-template <typename L>
-vector<double> stratified_sampling<L>::get_var() const {
- return sigma2;
-};
-
-template <typename L>
-void stratified_sampling<L>::print_mean() const {
- cout<<"les espérances :"<<endl;
- for(int i=0;i<I;i++){
- cout<<mean[i]<<"\t";
- }
- cout<<endl;
-};
-
-template <typename L>
-void stratified_sampling<L>::print_sigma() const {
- cout<<"les écarts types :"<<endl;
- for(int i=0;i<I;i++){
- cout<<sqrt(sigma2[i])<<"\t";
- }
- cout<<endl;
-};
-
-template <typename L>
-pair<double,double> stratified_sampling<L>::estimator() const {
- double est_mean = 0;
- double est_std = 0;
- for (int i=0; i<I; i++) {
- est_mean += mean[i]*p[i];
- est_std += sqrt(sigma2[i])*p[i];
- }
- return {est_mean, est_std};
-}