-
-
Save jeancroy/bfca4cf12a6a20eea9c0f1aa4f2d377e to your computer and use it in GitHub Desktop.
// 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 | |
} |
jeancroy
commented
Mar 5, 2019
•
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:
- Very narrow interval. It's unlikely you'll fit inside by luck.
- 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)
}