From 0bf1bbf3975a752675677a03b73b8b485640e3c5 Mon Sep 17 00:00:00 2001 From: Alexis David Jacq Date: Tue, 11 Jul 2017 16:12:31 +0100 Subject: [PATCH] Using Box-Muller for truncated normal sampling We may want to sample a normal truncated in a disk of ray (2*standard deviation) in order to initialize weights. The implementation is immediate with box-muller method. --- THRandom.c | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/THRandom.c b/THRandom.c index 55ee943..26fe56c 100644 --- a/THRandom.c +++ b/THRandom.c @@ -241,6 +241,31 @@ double THRandom_normal(THGenerator *_generator, double mean, double stdv) return _generator->normal_rho*sin(2.*M_PI*_generator->normal_x)*stdv+mean; } +double THRandom_normal_truncated(THGenerator *_generator, double mean, double stdv) +{ + THArgCheck(stdv > 0, 2, "standard deviation must be strictly positive"); + + /* Using Box-Muller for truncated normal sampling + between -2*std and 2*std. + (based on Martino, L.; Luengo, D.; Míguez, J. + "Efficient sampling from truncated bivariate Gaussians + via Box-Muller transformation")*/ + if(!_generator->normal_is_valid) + { + _generator->normal_x = __uniform__(_generator); + _generator->normal_y = __uniform__(_generator, 0, 1-exp(-2)); + _generator->normal_rho = sqrt(-2. * log(1.0-_generator->normal_y)); + _generator->normal_is_valid = 1; + } + else + _generator->normal_is_valid = 0; + + if(_generator->normal_is_valid) + return _generator->normal_rho*cos(2.*M_PI*_generator->normal_x)*stdv+mean; + else + return _generator->normal_rho*sin(2.*M_PI*_generator->normal_x)*stdv+mean; +} + double THRandom_exponential(THGenerator *_generator, double lambda) { return(-1. / lambda * log(1-__uniform__(_generator)));