From a858140e7775d52586c11fff665318ac3837cf54 Mon Sep 17 00:00:00 2001 From: Bertrand Date: Mon, 8 Feb 2016 21:16:09 +0000 Subject: vecteur des variances et des écarts types mis à jour MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- stratified_sampling.hpp | 21 +++++++++++++-------- 1 file 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 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 p = rtnorm(gen, a, b, mean, sigma); + pair 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 struct stratified_sampling { stratified_sampling(vector p, vector 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 get_mean(); @@ -58,7 +58,7 @@ private: vector M; vector cumM; vector mean; - vector sigma; + vector sigma2; vector gen; }; @@ -88,6 +88,11 @@ void stratified_sampling::update(int Nk) { } } else { + //On remplit un vecter des écarts types à parti de notre vecteur de variance + std::vector 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::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]; } }; -- cgit v1.2.3-70-g09d2