# Ziggurat method for Normal Random Numbers

8 replies to this topic

### #1 andybiochem

andybiochem

Universalist

• Active Members
• 308 posts

Posted 23 December 2009 - 06:18 PM

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:

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

; #VARIABLES# ===============================================================================
Global \$s_adZigX[128 + 1]
_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

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

### #2 Malkey

Malkey

• Active Members
• 1,297 posts

Posted 24 December 2009 - 06:37 AM

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.

AutoIt
```#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)```

### #3 andybiochem

andybiochem

Universalist

• Active Members
• 308 posts

Posted 24 December 2009 - 01:42 PM

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!

### #4 andybiochem

andybiochem

Universalist

• Active Members
• 308 posts

Posted 24 December 2009 - 02:02 PM

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

### #5 ValeryVal

ValeryVal

Dig the way

• Active Members
• 422 posts

Posted 24 December 2009 - 02:05 PM

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

### #6 andybiochem

andybiochem

Universalist

• Active Members
• 308 posts

Posted 24 December 2009 - 02:43 PM

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.

Edited by andybiochem, 24 December 2009 - 02:53 PM.

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

### #7 Mat

Mat

43 38 48 31 30 4E 34 4F 32

• MVPs
• 4,042 posts

Posted 24 December 2009 - 03:15 PM

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

I don't know where I'm going, but I'm on my way.

### #8 andybiochem

andybiochem

Universalist

• Active Members
• 308 posts

Posted 24 December 2009 - 05:25 PM

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!

### #9 Skrip

Skrip

Psychonaut

• Active Members
• 2,340 posts

Posted 25 December 2009 - 05:53 AM

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.

We're trapped in the belly of this horrible machine.And the machine is bleeding to death...

#### 0 user(s) are reading this topic

0 members, 0 guests, 0 anonymous users