nettoyage du dépôt

1 files changed, 238 insertions, 0 deletions

diff --git a/src/rtnorm.cpp b/src/rtnorm.cpp
new file mode 100644
index 0000000..721d243
--- /dev/null
+++ b/src/rtnorm.cpp
@@ -0,0 +1,238 @@
+//  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 <cmath>
+#include <iostream>
+#include <gsl/gsl_rng.h>
+#include <gsl/gsl_randist.h>
+#include <gsl/gsl_sf_erf.h>
+
+#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<double, double> 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 ! ***"<<std::endl;
+    exit(1);
+  }
+
+  // Check if |a| < |b|
+  else if(fabs(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<xmin)
+  {
+    while(!stop)
+    {
+      r = gsl_ran_gaussian(gen,1);
+      stop = (r>=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<xmax))
+      {
+          
+        // Two leftmost and rightmost regions
+        sim = Rtnorm::x[k] + (Rtnorm::x[k+1]-Rtnorm::x[k]) * gsl_rng_uniform(gen);
+
+        if ((sim >= a) && (sim <= b))
+        {
+          // Accept this proposition, otherwise reject
+          simy = Rtnorm::yu[k]*gsl_rng_uniform(gen);
+          if ( (simy<yl(k)) || (sim * sim + 2*log(simy) + ALPHA) < 0 )
+          {
+            r = sim;
+            stop = true;
+          }
+        }        
+      }
+
+      else // All the other boxes
+      {    
+        u = gsl_rng_uniform(gen);
+        simy = Rtnorm::yu[k] * u;
+        d = Rtnorm::x[k+1] - Rtnorm::x[k];
+        ylk = yl(k);
+        if(simy < ylk)  // That's what happens most of the time 
+        {  
+          r = Rtnorm::x[k] + u*d*Rtnorm::yu[k]/ylk;  
+          stop = true;
+        }
+        else
+        {
+          sim = Rtnorm::x[k] + d * gsl_rng_uniform(gen);
+
+          // Otherwise, check you're below the pdf curve
+          if((sim * sim + 2*log(simy) + ALPHA) < 0)
+          {
+            r = sim;
+            stop = true;
+          }
+        }
+        
+      }
+    }  
+  }
+  
+  //-----------------------------
+
+  // Scaling
+  if(mu!=0 || sigma!=1)
+    r = r*sigma + mu;
+  
+  // Compute the probability
+  Z = sqpi *sq2 * sigma * ( gsl_sf_erf(b*sq2) - gsl_sf_erf(a*sq2) );
+  p = exp(-pow((r-mu)/sigma,2)/2) / Z; 
+
+  return std::pair<double,double>(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;
+}
+