Opened 12 days ago

Last modified 11 days ago

#3027 reviewing enhancement

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

Reported by: gkronber Owned by: abeham
Priority: medium Milestone: HeuristicLab 3.3.17
Component: Random Version: branch
Keywords: Cc:

Description (last modified by 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.

Attachments (2)

empirical distributions.png (33.7 KB) - added by gkronber 12 days ago.
empirical distributions_polar.png (28.0 KB) - added by gkronber 12 days ago.

Download all attachments as: .zip

Change History (18)

comment:1 Changed 12 days ago by 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;

comment:2 Changed 12 days ago by gkronber

Performance comparison (1mio samples, release build):

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

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

Changed 12 days ago by gkronber

comment:3 Changed 12 days ago by gkronber

  • Description modified (diff)

comment:4 Changed 12 days ago by gkronber

  • Description modified (diff)

comment:5 Changed 12 days ago by gkronber

  • Owner set to gkronber
  • Status changed from new to accepted

comment:6 Changed 12 days ago by gkronber

r17243: create branch.

comment:7 Changed 12 days ago by gkronber

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

comment:8 Changed 12 days ago by 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();
    for(int i=0;i<N;i++) {
      s_r1[i] = r1.NextDouble();
    Console.WriteLine("{0}", stopwatch.ElapsedMilliseconds);
    for(int i=0;i<N;i++) {
      s_r2[i] = r2.NextDouble();
    Console.WriteLine("{0}", stopwatch.ElapsedMilliseconds);
    alglib.hqrndstate state;
    alglib.hqrndseed(1234, 0, out state);
    for(int i=0;i<N;i++) {
      s_r3[i] = alglib.hqrndnormal(state);
    Console.WriteLine("{0}", stopwatch.ElapsedMilliseconds);

comment:9 Changed 12 days ago by gkronber

Runtime (1mio samples, release mode):

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

comment:10 Changed 12 days ago by abeham

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

comment:12 Changed 12 days ago by gkronber

New result:

Last edited 12 days ago by gkronber (previous) (diff)

comment:13 Changed 12 days ago by gkronber

  • Owner changed from gkronber to abeham
  • Status changed from accepted to reviewing
  • Version set to branch

Changed 12 days ago by gkronber

comment:14 Changed 12 days ago by gkronber

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

comment:15 Changed 12 days ago by 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.

comment:16 Changed 11 days ago by abeham

  • Milestone changed from HeuristicLab 3.3.x Backlog to HeuristicLab 3.3.17
Note: See TracTickets for help on using tickets.