Provide a simple, efficient, and correct implementation to sample Gaussian random variables

[gkronber]

We currently use the Ziggurat method which is fast but not straight-forward. abeham made some tests and found significant deviations of the empirical distributions of samples generated with our implementation compared to other implementations.

In this image an offset is added to the probability values to prevent overlapping of points. The issue is that the blue line stops at approx. -3.5 and 3.5. So the current implementation does not produce correct distributions in the tails.

[gkronber]

Code from Sim#:

    const double NormalMagicConst = 4 * Math.Exp(-0.5) / Math.Sqrt(2.0);
    public double RandNormal(double mu, double sigma) {
      double z, zz, u1, u2;
      do {
        u1 = Random.NextDouble();
        u2 = 1 - Random.NextDouble();
        z = NormalMagicConst * (u1 - 0.5) / u2;
        zz = z * z / 4.0;
      } while (zz > -Math.Log(u2));
      return mu + z * sigma;
    }

[gkronber]

Performance comparison (1mio samples, release build):

FastRandom:
N(0,1)-hl: 28,7
N(0,1)-Sim#: 53,9
N(0,1)-alglib: 69,7

MersenneTwister:
N(0,1)-hl: 92,3
N(0,1)-Sim#: 174,1

[gkronber]

r17243: create branch.

[gkronber]

r17244: added class to generate normally distributed numbers using the polar method.

[gkronber]

Test script:

       
using System;
using System.Linq;
using System.Collections.Generic;

using HeuristicLab.Core;
using HeuristicLab.Common;
using HeuristicLab.Analysis;
using HeuristicLab.Analysis.Views;
using HeuristicLab.Random;

public class MyScript : HeuristicLab.Scripting.CSharpScriptBase {
  public override void Main() {
    var r1 = new NormalDistributedRandom(new FastRandom(1234), 0, 1 );
    var r2 = new NormalDistributedRandomPolar(new FastRandom(1234), 0, 1 );

    var dt = new DataTable();
    int N = 1000000;
    var s_r1 = new double[N];
    var s_r2 = new double[N];
    var s_r3 = new double[N];
    
    var stopwatch = new System.Diagnostics.Stopwatch();
    stopwatch.Start();          
    for(int i=0;i<N;i++) {
      s_r1[i] = r1.NextDouble();
    }
    stopwatch.Stop();
    Console.WriteLine("{0}", stopwatch.ElapsedMilliseconds);
    stopwatch.Restart();
    for(int i=0;i<N;i++) {
      s_r2[i] = r2.NextDouble();
    }
    stopwatch.Stop();
    Console.WriteLine("{0}", stopwatch.ElapsedMilliseconds);
    stopwatch.Restart();
    alglib.hqrndstate state;
    alglib.hqrndseed(1234, 0, out state);
    for(int i=0;i<N;i++) {
      s_r3[i] = alglib.hqrndnormal(state);
    }
    stopwatch.Stop();
    Console.WriteLine("{0}", stopwatch.ElapsedMilliseconds);
    
  }
}

[gkronber]

Runtime (1mio samples, release mode):

 -  24ms, NormalDistributedRandom / FastRandom
 -  36ms, NormalDistributedRandomPolar / FastRandom
 -  61ms, alglib
 -  70ms, NormalDistributedRandom / MersenneTwister
 - 129ms, NormalDistributedRandomPolar / MersenneTwister

[abeham]

FYI, I also changed it to use Marsaglia-Polar method in one of the recent Sim# releases.

[gkronber]

r17245: change code to match ​https://gcc.gnu.org/onlinedocs/gcc-4.6.0/libstdc++/api/a01001_source.html#l01649

[gkronber]

New result:

[gkronber]

It seems there is still a mismatch to alglib in the tails even with MersenneTwister

[abeham]

In Sim# I use a newer RNG for uniformly distributed random numbers called permuted congruential generators or ​PCG for short.

I looked at r17245, good spot, I think that's necessary in order to produce exactly mu as random number.

[abeham]

Sorry, for the long delay! Ok.

[gkronber]

• r17806: merged branch to trunk
• r17807: deleted branch

[gkronber]

r17810: add StorableType attribute

[gkronber]

r17864: merged r17806, r17810 from trunk to stable