From 69ef15f06495b92bd619963eb543cb098bab479c Mon Sep 17 00:00:00 2001 From: Bertrand Date: Tue, 1 Mar 2016 14:32:31 +0000 Subject: on n'utilise plus rtnorm --- src/rtnorm.cpp | 238 --------------------------------------------------------- 1 file changed, 238 deletions(-) delete mode 100644 src/rtnorm.cpp (limited to 'src/rtnorm.cpp') diff --git a/src/rtnorm.cpp b/src/rtnorm.cpp deleted file mode 100644 index 721d243..0000000 --- a/src/rtnorm.cpp +++ /dev/null @@ -1,238 +0,0 @@ -// Pseudorandom numbers from a truncated Gaussian distribution. -// -// This implements an extension of Chopin's algorithm detailed in -// N. Chopin, "Fast simulation of truncated Gaussian distributions", -// Stat Comput (2011) 21:275-288 -// -// Copyright (C) 2012 Guillaume Dollé, Vincent Mazet -// (LSIIT, CNRS/Université de Strasbourg) -// Version 2012-07-04, Contact: vincent.mazet@unistra.fr -// -// 06/07/2012: -// - first launch of rtnorm.cpp -// -// Licence: GNU General Public License Version 2 -// This program is free software; you can redistribute it and/or modify it -// under the terms of the GNU General Public License as published by the -// Free Software Foundation; either version 2 of the License, or (at your -// option) any later version. This program is distributed in the hope that -// it will be useful, but WITHOUT ANY WARRANTY; without even the implied -// warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -// GNU General Public License for more details. You should have received a -// copy of the GNU General Public License along with this program; if not, -// see http://www.gnu.org/licenses/old-licenses/gpl-2.0.txt -// -// Depends: LibGSL -// OS: Unix based system - - -#include -#include -#include -#include -#include - -#include "rtnorm.hpp" - -int N = 4001; // Index of the right tail - -//------------------------------------------------------------ -// Pseudorandom numbers from a truncated Gaussian distribution -// The Gaussian has parameters mu (default 0) and sigma (default 1) -// and is truncated on the interval [a,b]. -// Returns the random variable x and its probability p(x). -std::pair rtnorm( - gsl_rng *gen, - double a, - double b, - const double mu, - const double sigma) -{ - // Design variables - double xmin = -2.00443204036; // Left bound - double xmax = 3.48672170399; // Right bound - int kmin = 5; // if kb-ka < kmin then use a rejection algorithm - double INVH = 1631.73284006; // = 1/h, h being the minimal interval range - int I0 = 3271; // = - floor(x(0)/h) - double ALPHA = 1.837877066409345; // = log(2*pi) - int xsize=sizeof(Rtnorm::x)/sizeof(double); // Length of table x - int stop = false; - double sq2 = 7.071067811865475e-1; // = 1/sqrt(2) - double sqpi = 1.772453850905516; // = sqrt(pi) - - double r, z, e, ylk, simy, lbound, u, d, sim, Z, p; - int i, ka, kb, k; - - // Scaling - if(mu!=0 || sigma!=1) - { - a=(a-mu)/sigma; - b=(b-mu)/sigma; - } - - //----------------------------- - - // Check if a < b - if(a>=b) - { - std::cerr<<"*** B must be greater than A ! ***"<fabs(b)) - r = -rtnorm(gen,-b,-a).first; // Pair (r,p) - - // If a in the right tail (a > xmax), use rejection algorithm with a truncated exponential proposal - else if(a>xmax) - r = rtexp(gen,a,b); - - // If a in the left tail (a < xmin), use rejection algorithm with a Gaussian proposal - else if(a=a) && (r<=b); - } - } - - // In other cases (xmin < a < xmax), use Chopin's algorithm - else - { - // Compute ka - i = I0 + floor(a*INVH); - ka = Rtnorm::ncell[i]; - - // Compute kb - (b>=xmax) ? - kb = N : - ( - i = I0 + floor(b*INVH), - kb = Rtnorm::ncell[i] - ); - - // If |b-a| is small, use rejection algorithm with a truncated exponential proposal - if(abs(kb-ka) < kmin) - { - r = rtexp(gen,a,b); - stop = true; - } - - while(!stop) - { - // Sample integer between ka and kb - k = floor(gsl_rng_uniform(gen) * (kb-ka+1)) + ka; - - if(k == N) - { - // Right tail - lbound = Rtnorm::x[xsize-1]; - z = -log(gsl_rng_uniform(gen)); - e = -log(gsl_rng_uniform(gen)); - z = z / lbound; - - if ((pow(z,2) <= 2*e) && (z < b-lbound)) - { - // Accept this proposition, otherwise reject - r = lbound + z; - stop = true; - } - } - - else if ((k <= ka+1) || (k>=kb-1 && b= a) && (sim <= b)) - { - // Accept this proposition, otherwise reject - simy = Rtnorm::yu[k]*gsl_rng_uniform(gen); - if ( (simy(r,p); -} - -//------------------------------------------------------------ -// Compute y_l from y_k -double yl(int k) -{ - double yl0 = 0.053513975472; // y_l of the leftmost rectangle - double ylN = 0.000914116389555; // y_l of the rightmost rectangle - - if (k == 0) - return yl0; - - else if(k == N-1) - return ylN; - - else if(k <= 1953) - return Rtnorm::yu[k-1]; - - else - return Rtnorm::yu[k+1]; -} - -//------------------------------------------------------------ -// Rejection algorithm with a truncated exponential proposal -double rtexp(gsl_rng *gen, double a, double b) -{ - int stop = false; - double twoasq = 2*pow(a,2); - double expab = exp(-a*(b-a)) - 1; - double z, e; - - while(!stop) - { - z = log(1 + gsl_rng_uniform(gen)*expab); - e = -log(gsl_rng_uniform(gen)); - stop = (twoasq*e > pow(z,2)); - } - return a - z/a; -} - -- cgit v1.2.3-70-g09d2