Ziggurat method for Normal Random Numbers

andybiochem · December 23, 2009

Hi!

This is a small UDF to allow the creation of normally distributed random numbers using a user-specified mean and standard deviation.

I have previously released a gaussian number generator that uses the Box-Muller Polar Transform method here.

The method posted here is the 'Ziggurat' method, as described by JURGEN A DOORNIK, Marsaglia and Tsang here.

The Ziggurat method has an improvement of 25-30% in speed compared to the Box-Muller method.

Example:

- Generation of 1000 random numbers with mean = 0 and standard deviation = 1

#include <Array.au3>
#include "Ziggurat.au3"

Global $Array[1001]

For $i = 1 To 1000
    $Array[$i] = _Random_Gaussian_Zig(0, 1)
Next

_ArrayDisplay($Array,"1000 Gaussian Distribution numbers")

UDF:

Ziggurat.au3

;#AutoIt3Wrapper_Au3Check_Parameters=-d -w 1 -w 2 -w 3 -w 4 -w 5 -w 6

; #VARIABLES# ===============================================================================
Global $s_adZigX[128 + 1]
Global $s_adZigR[128]
_Random_Gaussian_Zig_Init() ;initiate tables
; ===========================================================================================



; #FUNCTION# ;===============================================================================
; Name...........: _Random_Gaussian_Zig_Init
; Description ...: Creates the look-up tables/arrays necessary for ziggurat rejection
; Syntax.........: _Random_Gaussian_Zig_Init()
; Return values .: Success - Returns two look-up arrays
;                  Failure - (under construction)
; Author ........: JURGEN A DOORNIK, Marsaglia and Tsang, AndyBiochem (AutoIt conversion)
; ===========================================================================================
Func _Random_Gaussian_Zig_Init()

    Local $iC = 128
    Local $dR = 3.442619855899
    Local $dV = 9.91256303526217E-3
    Local $i, $f

    $f = Exp(-0.5 * $dR * $dR)
    $s_adZigX[0] = $dV / $f
    $s_adZigX[1] = $dR
    $s_adZigX[$iC] = 0

    For $i = 2 To $iC - 1
        $s_adZigX[$i] = Sqrt(-2 * Log($dV / $s_adZigX[$i - 1] + $f));
        $f = Exp(-0.5 * $s_adZigX[$i] * $s_adZigX[$i]);
    Next

    For $i = 0 To $iC - 1
        $s_adZigR[$i] = $s_adZigX[$i + 1] / $s_adZigX[$i];
    Next

EndFunc   ;==>_Random_Gaussian_Zig_Init

; #FUNCTION# ;===============================================================================
; Name...........: _Random_Gaussian_Zig
; Description ...: Generates a gaussian distributed normal variable.
; Syntax.........: _Random_Gaussian_Zig($iMean,$iSD)
; Parameters ....: $iMean - The mean of the gaussian distribution
;                  $iSD - The standard deviation of the gaussian distribution
; Return values .: Success - Returns a single gaussian random number
;                  Failure - (under construction)
; Author ........: JURGEN A DOORNIK, Marsaglia and Tsang, AndyBiochem (AutoIt conversion)
; ===========================================================================================
Func _Random_Gaussian_Zig($iMean, $iSD)

    Local $i, $x, $u, $f0, $f1, $y, $xOut

    While 1
        $u = 2 * Random() - 1
        $i = Random(0, 128)

        If Abs($u) < $s_adZigR[$i] Then
            $xOut = $u * $s_adZigX[$i]
            ExitLoop
        EndIf

        If $i = 0 Then
            While (-2 * $y) < ($x * $x)
                $x = Log(Random()) / 3.442619855899
                $y = Log(Random())
            WEnd
            If $u < 0 Then
                $xOut = $x - 3.442619855899
                ExitLoop
            EndIf
            $xOut = 3.442619855899 - $x
            ExitLoop
        EndIf

        $x = $u * $s_adZigX[$i]
        $f0 = Exp(-0.5 * ($s_adZigX[$i] * $s_adZigX[$i] - $x * $x))
        $f1 = Exp(-0.5 * ($s_adZigX[$i + 1] * $s_adZigX[$i + 1] - $x * $x))
        If ($f1 + Random() * ($f0 - $f1) < 1.0) Then
            $xOut = $x
            ExitLoop
        EndIf
    WEnd

    $xOut *= $iSD

    Return $iMean + $xOut

EndFunc   ;==>_Random_Gaussian_Zig

Any comments or critisism welcome!

Malkey · December 24, 2009

Thought I would test your claims of _Random_Gaussian_Zig(0, 1) having a means of 0 and standard deviation of 1. And check its speed compared to say the built-in AutoIt Random function.

Results

Gaussian 0.8307secs Mean = 0.00970560903677185; s = 0.995471473602324
Random() 0.1025secs Mean = 0.00539179714763396; s = 0.999977443170705

Your function appear to be working satisfactorily.

Through trial and error, I noticed numbers generated with Random(-1.73, 1.73) have a standard deviation of 1. And Random(-1.73 multiplied by 2, 1.73 multiplied by 2) produce a standard deviation of 2, approximately.

#include <Array.au3>
#include "Ziggurat.au3"

Global $Array[10000], $iSum, $iSum2, $Array2[10000], $iSumDev, $iMean, $iSumDev2, $iMean2

;========== Gaussian ======================
$begin = TimerInit()
For $i = 0 To 9999
    $Array[$i] = _Random_Gaussian_Zig(0, 1)
    $iSum += $Array[$i]
Next
;_ArrayDisplay($Array,"1000 Gaussian Distribution numbers")
$iMean = $iSum / 10000
ConsoleWrite("Gaussian " & Round(TimerDiff($begin) / 1000, 4) & "secs   Mean = " & $iMean & "; s = ")
For $i = 0 To 9999
    $iSumDev += ($Array[$i] - $iMean) ^ 2
Next
ConsoleWrite(($iSumDev / 9999) ^ 0.5 & @CRLF)

;======== Random() =========================
$begin2 = TimerInit()
For $i = 0 To 9999
    $Array2[$i] = Random(-1.73, 1.73)
    $iSum2 += $Array2[$i]
Next
$iMean2 = $iSum2 / 10000
ConsoleWrite("Random() " & Round(TimerDiff($begin2) / 1000, 4) & "secs  Mean = " & $iMean2 & "; s = ")
For $i = 0 To 9999
    $iSumDev2 += ($Array2[$i] - $iMean2) ^ 2
Next
ConsoleWrite(($iSumDev2 / 9999) ^ 0.5 & @CRLF)

andybiochem · December 24, 2009

Thought I would test your claims of _Random_Gaussian_Zig(0, 1) having a means of 0 and standard deviation of 1. And check its speed compared to say the built-in AutoIt Random function.
....etc...

Perhaps I'm not understanding your intentions correctly, but I don't see the point in comparing _Random_Gaussian_Zig() to Random().

The whole point of _Random_Gaussian_Zig() is to produce a data set with gaussian distribution.

Random() produces a data set with uniform distribution.

i.e. in the _Random_Gaussian_Zig() data set 68% of the returned values will be within 1SD of the mean, 95% within 2SD etc etc. Whilst Random(-1.73, 1.73) might have a SD of 1, this is co-incidental as the distribution of data points is still uniform...

_Random_Gaussian_Zig() would pass a test for normality, whereas Random() would not.

Thanks for testing though!

andybiochem · December 24, 2009

...Having thought about this briefly, perhaps the problem here is with the topic title "Ziggurat method for Normal Random Numbers"

In this case "Normal" does not mean 'common, usual, run-of-the-mill' etc.

"Normal" here is intended in the mathematical sense, i.e. that the data fits to a bell-curve, or Gaussian Function.

If anyone thinks the title is misleading, I'll change it.

ValeryVal · December 24, 2009

The normal number generator is very simple in AutoIt.

Global Const $pi = 4*ATan(1)

$r1 = Random() ; first random value

$r2 = Random() ; second random value

$n1 = Cos(2*$pi*$r1) * Sqrt(-2*Log($r2)) ; first normal random value

$n2 = Sin(2*$pi*$r1) * Sqrt(-2*Log($r2)) ; second normal random value

These $n1 and $n2 were compared with Statistica values.

Very good result was obtained!

andybiochem · December 24, 2009

The normal number generator is very simple in AutoIt.

Global Const $pi = 4*ATan(1)

$r1 = Random() ; first random value
$r2 = Random() ; second random value

$n1 = Cos(2*$pi*$r1) * Sqrt(-2*Log($r2)) ; first normal random value
$n2 = Sin(2*$pi*$r1) * Sqrt(-2*Log($r2)) ; second normal random value

These $n1 and $n2 were compared with Statistica values.
Very good result was obtained!

wtf! Why have I never seen this method before?

I thought the Cos/Sin/Sqrt/Log would be quite expensive, but initial tests are showing it MUCH faster than Ziggurat!

Also, SD proportions look spot on! Thanks for that!

[EDIT]...

A quick google/wiki shows it to be the NON polar version of the Box-Muller method (see the polar version link in my sig). I think I must have assumed that Marsaglia's optimisation of Box-Muller was better, and not bothered with the original Box-Muller. :evil:

Edited December 24, 2009 by andybiochem

Mat · December 24, 2009

The worst thing ever when someone comes along and shows you the real way to solve the problem...

I learnt from it though, so it's got to be pretty good

Mat

andybiochem · December 24, 2009

The worst thing ever when someone comes along and shows you the real way to solve the problem...

No,that's the BEST thing ever! Being a scientist, I'd absolutely love it if someone would pick all my theories/algorithms to bits and show a better or improved way of doing things.

Ego is nothing, perfection is everything.

Skrip · December 25, 2009

No,that's the BEST thing ever! Being a scientist, I'd absolutely love it if someone would pick all my theories/algorithms to bits and show a better or improved way of doing things.
Ego is nothing, perfection is everything.

Great response. That's the way to go right there.

Sign In

Ziggurat method for Normal Random Numbers

Recommended Posts

andybiochem

Malkey

andybiochem

andybiochem

ValeryVal

andybiochem

Mat

andybiochem

Skrip

Create an account or sign in to comment

Create an account

Sign in

Recently Browsing 0 members

Browse

AutoIt Resources

Release

Beta