Sample Truncated Normal.cpp

// Sample from a normal distribution, truncated between -1 and 1

double MultiObjectiveCompactDifferentialEvolution::sampleValue(double mu, double sigma) {
	if (sigma < 1e-8) // Should avoid infinity error below
		return mu;
	
	try {
		double sqrt2Sigma = SQRT_2 * sigma;
		double erfMuNeg = boost::math::erf((mu - 1) / sqrt2Sigma);
		double erfMuPlus = boost::math::erf((mu + 1) / sqrt2Sigma);
		
		double u = randreal();
		double C = - boost::math::erf((mu + 1) / sqrt2Sigma) / (erfMuNeg - erfMuPlus);
		return mu - sqrt2Sigma * boost::math::erf_inv((u - C) * (erfMuNeg - erfMuPlus));
    
	} catch (exception& e) {
    
		cerr << "ERROR: Sampling method failed! mu: " << mu << ", sigma: " << sigma << ". Got: " << e.what() << endl;
    
	}
	
	return mu; // Safer guess
}

truncatedNormalSigma(mu, sigma, start, end){

// Rayleigh proposal truncated inside [a,b[
// https://www.iro.umontreal.ca/~lecuyer/myftp/papers/vt16truncnormal.pdf
// Original code produce sample in N(0,1) that need to be rescaled/translated to get N(mu,sigma^2)

scale = 1.0 / sigma;
a = (start - mu)*scale;
b = (end - mu)*scale;

c = 0.5*a*a;
d = 0.5*b*b;
q = 1.0 - exp( c - d );

while(true){
  u = rand(0,1.0);
  v = rand(0,1.0);
  x = c - ln(1.0-q*u);
  if( v*v*x <= a ) break;
}

return mu + sigma*sqrt(2*x)

}

truncatedNormal(mu, variance, start, end){


a = (start - mu);
b = (end - mu);
iv = 1.0 / variance;

aa = a*a*iv;
bb = b*b*iv;

c = 0.5*aa;
d = 0.5*bb;

q = 1.0 - exp( c - d );

while(true){
  u = rand(0,1.0);
  v = rand(0,1.0);
  x = c - ln(1.0-q*u);
  t = v*v*x;
  if( t*t <= aa ) break;
}

return mu + sqrt(2*variance*x)

}

// Simulate negative number as positive then flip sign.
//  if abs(b) < abs(b), swap ab, flip sign and re-flip sign at end 

  a = (a - mu) / sigma;
  b = (b - mu) / sigma;
  bma = b - a;


//  Marsaglia polar method N(0,1) between a and b
// The method generate two number and we use that as two opportunity to pass the range criteria.

       
        var u, v, s,t;
        
       while(true){
      
        while(true){
            u = Math.random() * 2 - 1;
            v = Math.random() * 2 - 1;
            s = u * u + v * v;
           if(s>0 && s<1) break;
        } 
        
        var scale = Math.sqrt(-2.0 * Math.log(s) / s); // log is ln
        
        t = v*scale 
        if( t >= a && t < b ) return t;

        t = u*scale 
       if( t >= a && t < b ) return t;
        }
      

//Robert Algorithm, for a<0<b 

 do{      
             u = Math.random();
             v = Math.random();
             r = a + (b-a)*v
      }
while( 2*ln(u) > -r*r  ); // u > exp(-0.5*r*r)
return r;

//  Robert Algorithm, for 0<a<b -- uniform generation, work best when b-a is small
//  To use when a and b are bellow 0, convert to a above 0 problem then flip sign at the end.

 do{      
             u = Math.random();
             v = Math.random();
             r = a + (b-a)*v
      }
while( 2*ln(u) > a*a -r*r  );
return r

How to generate normal number inside [a,b[ ?
The basic idea is to generate a Gaussian number, test interval and retry.

This will break under 2 conditions:

1. Very narrow interval. It's unlikely you'll fit inside by luck.
2. Very offset interval. The interval between 100 and 103 standard deviation is 3 standard deviation large, so quite large, but we are unlikely to shoot inside anyways.

To address (1) we need a process able to generate random number that's aware of [a,b]
To get there we flip the role, 1st we generate a uniform number inside [a,b] then reject if inconsistent with normal statistics.

To address (2) we'll use an algorithm used to simulate tail. (therefore whether [a,b[ pass by the mean will be an important question)
Some algorithm use both a and b to scale the problem, those are slower but uniform in speed & still faster than blind luck.
Some will generate from [a,inf[ then check b after the fact, where that strategy work they can be the fastest (but it'll still fail if too narrow)

input: mu,sigma, min,max

if(max==min) return max

if(max>min) flip


  a = (min - mu) / sigma;
  b = (max - mu) / sigma;
  bma = b - a;

// both are negative,
// flip sign, then reflip before return;
if( a<=0 && b<=0 ){
    tmp = a;
    a = -b;
    b = -tmp;
    sigma  = -sigma;
}

// Do the interval contains 0 ?
// Or is close to zero ?
// simulate bulk

if( a * b < 0 || a<2) {
   
    if(bma > 1.0){

       // interval cover at least one standard deviation ? 
      // simulate a gaussian first, then reject answers not in [a,b] 
       return sigma*GaussianProposal(a,b) + mu;

   }else{

      // interval is narrow ?
      // simulate the interval first, then reject answers not normal like.
       return sigma*UniformProposal (a,b) + mu;

   }

}
else{

    // We are far from 0, simulate right tail
   // (and flip for the left tail)

     return sigma*RaleighProposal(a,b) + mu;
}

GaussianProposal(a,b){

  var u, v, s,t;
        
       while(true){
      
        while(true){
            u = Math.random() * 2 - 1;
            v = Math.random() * 2 - 1;
            s = u * u + v * v;
           if(s>0 && s<1) break;
        } 
        
        var scale = Math.sqrt(-2.0 * Math.log(s) / s); // Math.log is natural logarithm
        
        t = v*scale 
        if( t >= a && t < b ) return t;

        t = u*scale 
       if( t >= a && t < b ) return t;
        }

}

UniformProposal (a,b) {

// k=0 if a<0<b 
// k = a*a if 0<a<b
// k = b*b if a<b<0
// k = a*b<0 ? 0 : Math.min(a*a, b*b)

  k = a<0?0:a*a;

 do{      
             u = Math.random();
             v = Math.random();
             r = a + (b-a)*v
      }
while( 2*ln(u) > k-r*r  ); // u > exp(0.5*k-0.5*r*r)
return r;

}

RaleighProposal (a,b) {

c = 0.5*a*a
q = 1.0 - exp( c - 0.5*b*b )

while(true){
  u = rand(0,1.0)
  v = rand(0,1.0)
  x = c - ln(1.0-q*u)
  if( v*v*x <= a ) break;
}
return sqrt(2*x)

}

One usage of truncated normal is
https://www.researchgate.net/publication/257680252_Compact_Differential_Evolution_Light_High_Performance_Despite_Limited_Memory_Requirement_and_Modest_Computational_Overhead