Opened 7 months ago

Last modified 5 weeks ago

#3027 readytorelease enhancement

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

Reported by: gkronber Owned by: gkronber
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 7 months ago.
empirical distributions_polar.png (28.0 KB) - added by gkronber 7 months ago.

Download all attachments as: .zip

Change History (19)

comment:1 Changed 7 months 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 7 months ago by 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

Changed 7 months ago by gkronber

comment:3 Changed 7 months ago by gkronber

  • Description modified (diff)

comment:4 Changed 7 months ago by gkronber

  • Description modified (diff)

comment:5 Changed 7 months ago by gkronber

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

comment:6 Changed 7 months ago by gkronber

r17243: create branch.

comment:7 Changed 7 months ago by gkronber

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

comment:8 Changed 7 months 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();
    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);
    
  }
}

comment:9 Changed 7 months 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 7 months ago by abeham

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

comment:12 Changed 7 months ago by gkronber

New result:

Last edited 7 months ago by gkronber (previous) (diff)

comment:13 Changed 7 months ago by gkronber

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

Changed 7 months ago by gkronber

comment:14 Changed 7 months ago by gkronber

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

comment:15 Changed 7 months 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 7 months ago by abeham

  • Milestone changed from HeuristicLab 3.3.x Backlog to HeuristicLab 3.3.17

comment:17 Changed 5 weeks ago by abeham

  • Owner changed from abeham to gkronber
  • Status changed from reviewing to readytorelease

Sorry, for the long delay! Ok.

Note: See TracTickets for help on using tickets.