diff options
Diffstat (limited to 'stratified_sampling.hpp')
| -rw-r--r-- | stratified_sampling.hpp | 21 |
1 files changed, 13 insertions, 8 deletions
diff --git a/stratified_sampling.hpp b/stratified_sampling.hpp index 0b0b4ef..e25ac63 100644 --- a/stratified_sampling.hpp +++ b/stratified_sampling.hpp @@ -22,33 +22,33 @@ typedef var_alea<double> var_alea_real; struct gaussian_truncated : public var_alea_real { - gaussian_truncated(double a, double b, double mean=0, double sigma=1, int seed=0) - :a(a), b(b), mean(mean), sigma(sigma), seed(seed) { + 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),sigma(other.sigma){ + :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, sigma); + pair<double, double> p = rtnorm(gen, a, b, mean, sigma2); return value = p.first; }; ~gaussian_truncated() { gsl_rng_free(gen); } private: int seed; - double mean, sigma, a, b; + double mean, sigma2, a, b; gsl_rng *gen; }; template <typename Gen> struct stratified_sampling { stratified_sampling(vector<double> p, vector<Gen> gen) - :p(p), gen(gen), mean(p.size(), 0), sigma(p.size(), 0){}; + :p(p), gen(gen), mean(p.size(), 0), sigma2(p.size(), 0){}; void update(int N); void draw(); vector<double> get_mean(); @@ -58,7 +58,7 @@ private: vector<int> M; vector<int> cumM; vector<double> mean; - vector<double> sigma; + vector<double> sigma2; vector<Gen> gen; }; @@ -88,6 +88,11 @@ void stratified_sampling<Gen>::update(int Nk) { } } 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; @@ -117,7 +122,7 @@ void stratified_sampling<Gen>::draw() { } oldmean=mean[i]; mean[i]=(mean[i]*cumM[i]+m)/(cumM[i]+M[i]); - sigma[i]=((sigma[i]+oldmean*oldmean)*cumM[i] + s)/(cumM[i]+M[i]) - mean[i]*mean[i]; + sigma2[i]=((sigma2[i]+oldmean*oldmean)*cumM[i] + s)/(cumM[i]+M[i]) - mean[i]*mean[i]; } }; |
