8 replies to this topic

[andybiochem]

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   2.92K   211 downloads

AutoIt        expandcollapse  popup

;#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]

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        expandcollapse  popup

#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]

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]

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

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]

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.

[Mat]

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]

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]

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.