Sign in to follow this  
Followers 0
andybiochem

Ziggurat method for Normal Random Numbers

9 posts in this topic

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!


- Table UDF - create simple data tables - Line Graph UDF GDI+ - quickly create simple line graphs with x and y axes (uses GDI+ with double buffer) - Line Graph UDF - quickly create simple line graphs with x and y axes (uses AI native graphic control) - Barcode Generator Code 128 B C - Create the 1/0 code for barcodes. - WebCam as BarCode Reader - use your webcam to read barcodes - Stereograms!!! - make your own stereograms in AutoIT - Ziggurat Gaussian Distribution RNG - generate random numbers based on normal/gaussian distribution - Box-Muller Gaussian Distribution RNG - generate random numbers based on normal/gaussian distribution - Elastic Radio Buttons - faux-gravity effects in AutoIT (from javascript)- Morse Code Generator - Generate morse code by tapping your spacebar!

Share this post


Link to post
Share on other sites



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)

Share this post


Link to post
Share on other sites

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!


- Table UDF - create simple data tables - Line Graph UDF GDI+ - quickly create simple line graphs with x and y axes (uses GDI+ with double buffer) - Line Graph UDF - quickly create simple line graphs with x and y axes (uses AI native graphic control) - Barcode Generator Code 128 B C - Create the 1/0 code for barcodes. - WebCam as BarCode Reader - use your webcam to read barcodes - Stereograms!!! - make your own stereograms in AutoIT - Ziggurat Gaussian Distribution RNG - generate random numbers based on normal/gaussian distribution - Box-Muller Gaussian Distribution RNG - generate random numbers based on normal/gaussian distribution - Elastic Radio Buttons - faux-gravity effects in AutoIT (from javascript)- Morse Code Generator - Generate morse code by tapping your spacebar!

Share this post


Link to post
Share on other sites

...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.


- Table UDF - create simple data tables - Line Graph UDF GDI+ - quickly create simple line graphs with x and y axes (uses GDI+ with double buffer) - Line Graph UDF - quickly create simple line graphs with x and y axes (uses AI native graphic control) - Barcode Generator Code 128 B C - Create the 1/0 code for barcodes. - WebCam as BarCode Reader - use your webcam to read barcodes - Stereograms!!! - make your own stereograms in AutoIT - Ziggurat Gaussian Distribution RNG - generate random numbers based on normal/gaussian distribution - Box-Muller Gaussian Distribution RNG - generate random numbers based on normal/gaussian distribution - Elastic Radio Buttons - faux-gravity effects in AutoIT (from javascript)- Morse Code Generator - Generate morse code by tapping your spacebar!

Share this post


Link to post
Share on other sites

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!


The point of world view

Share this post


Link to post
Share on other sites

#6 ·  Posted (edited)

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 by andybiochem

- Table UDF - create simple data tables - Line Graph UDF GDI+ - quickly create simple line graphs with x and y axes (uses GDI+ with double buffer) - Line Graph UDF - quickly create simple line graphs with x and y axes (uses AI native graphic control) - Barcode Generator Code 128 B C - Create the 1/0 code for barcodes. - WebCam as BarCode Reader - use your webcam to read barcodes - Stereograms!!! - make your own stereograms in AutoIT - Ziggurat Gaussian Distribution RNG - generate random numbers based on normal/gaussian distribution - Box-Muller Gaussian Distribution RNG - generate random numbers based on normal/gaussian distribution - Elastic Radio Buttons - faux-gravity effects in AutoIT (from javascript)- Morse Code Generator - Generate morse code by tapping your spacebar!

Share this post


Link to post
Share on other sites

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

Share this post


Link to post
Share on other sites

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.


- Table UDF - create simple data tables - Line Graph UDF GDI+ - quickly create simple line graphs with x and y axes (uses GDI+ with double buffer) - Line Graph UDF - quickly create simple line graphs with x and y axes (uses AI native graphic control) - Barcode Generator Code 128 B C - Create the 1/0 code for barcodes. - WebCam as BarCode Reader - use your webcam to read barcodes - Stereograms!!! - make your own stereograms in AutoIT - Ziggurat Gaussian Distribution RNG - generate random numbers based on normal/gaussian distribution - Box-Muller Gaussian Distribution RNG - generate random numbers based on normal/gaussian distribution - Elastic Radio Buttons - faux-gravity effects in AutoIT (from javascript)- Morse Code Generator - Generate morse code by tapping your spacebar!

Share this post


Link to post
Share on other sites

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.


[left][sub]We're trapped in the belly of this horrible machine.[/sub][sup]And the machine is bleeding to death...[/sup][sup][/sup][/left]

Share this post


Link to post
Share on other sites

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!


Register a new account

Sign in

Already have an account? Sign in here.


Sign In Now
Sign in to follow this  
Followers 0