Opened 2 years ago

Closed 7 months ago

#3027 closed enhancement (done)

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: trunk
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 2 years ago.
empirical distributions_polar.png (28.0 KB) - added by gkronber 2 years ago.

Download all attachments as: .zip

Change History (24)

comment:1 Changed 2 years 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 2 years 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 2 years ago by gkronber

comment:3 Changed 2 years ago by gkronber

  • Description modified (diff)

comment:4 Changed 2 years ago by gkronber

  • Description modified (diff)

comment:5 Changed 2 years ago by gkronber

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

comment:6 Changed 2 years ago by gkronber

r17243: create branch.

comment:7 Changed 2 years ago by gkronber

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

comment:8 Changed 2 years 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 2 years 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 2 years ago by abeham

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

comment:12 Changed 2 years ago by gkronber

New result:

Last edited 2 years ago by gkronber (previous) (diff)

comment:13 Changed 2 years ago by gkronber

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

Changed 2 years ago by gkronber

comment:14 Changed 2 years ago by gkronber

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

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

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

comment:17 Changed 20 months ago by abeham

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

Sorry, for the long delay! Ok.

comment:18 Changed 10 months ago by gkronber

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

comment:19 Changed 10 months ago by gkronber

  • Version changed from branch to trunk

comment:20 Changed 10 months ago by gkronber

r17810: add StorableType attribute

comment:21 Changed 8 months ago by gkronber

r17864: merged r17806, r17810 from trunk to stable

comment:22 Changed 7 months ago by gkronber

  • Resolution set to done
  • Status changed from readytorelease to closed
Note: See TracTickets for help on using tickets.