Normally Distributed Random Numbers (Tcl recipe) by William Kappele

ActiveState Code (http://code.activestate.com/recipes/143083/)

I am not sure why the commands

set u1 rand()

set u2 rand()

are there. Semantically they have no effect. The code in the recipe is equivalent to

set dRandNormal1 [expr {$stdDev * sqrt(-2 * log(rand())) * cos(2 * 3.14159 * rand())}]

set dRandNormal2 [expr {$stdDev * sqrt(-2 * log(rand())) * sin(2 * 3.14159 * rand())}]

Code Explanation. u1 and u2 need to be the same in both calculations, so I chose them up front. If you don't keep them the same, there is no guarantee that the normal random numbers will be independent.

Bill Kappele (bill@mathoptions.com)

Code Correction. Thank you, Andreas. The code should read

set u1 [expr rand()] set u2 [expr rand()]

to function correctly.

Need to be the same. Need a third variable to make both the same

set temp [expr rand()]

set u1 $temp

set u2 $temp

Two variables. I believe two variables are needed. I am not familiar with the current language you are using, but the variables in the formulas must be uncorrelated to provide gaussian like distribution. i.e.:

a = rand_0_1() ;

b = rand_0_1() ;

y1 = sigma * sqrt(-2.0 * a )*cos(2.0 * pi * b ) ;

y2 = sigma * sqrt(-2.0 * a )*sin(2.0 * pi * b ) ;

the random numbers y1 and y2 are now orthogonal and normally distributed. I think.

Thanks for the recipe. I needed to doublecheck my code.

Setting the two variables beforehand isn't needed to keep the resulting random numbers independent. When you think about it, andreas' numbers can't be dependent since there's no variables in his two expressions and each rand() result is indepedent of another. What it does is that it saves two rand() calls. That probably doesn't matter until you change rand() to something that takes entropy from your system, like reading from /dev/random, and you waste half of it.

If TCL's rand() returned a number in the range [0,1) (including 0 but not 1), which I've always thought, u1 or u2 would eventually become 0 and crash at the log() part. log(0) = infinity. The Box-Muller tranform expects a range of (0,1] (not including 0 but 1) and you could get that with 1 - rand(). I just read that the range is (0,1) (not including either 0 or 1) so it probably doesn't matter, but it's a safe thing to do in case it changes. I'm pretty sure rand() will never return 1.

Use acos(-1) for pi instead of typing digits. That way, when computers become 128-bit, 256-bit etc, your pi precision will follow.

set pi [expr {acos(-1)}]

proc randNormal {stdDev} {
    global pi

    set u1 [expr {1 - rand()}]
    set u2 [expr {1 - rand()}]

    set dRandNormal1 [expr $stdDev * sqrt(-2 * log($u1)) * cos(2 * $pi * $u2)]
    set dRandNormal2 [expr $stdDev * sqrt(-2 * log($u1)) * sin(2 * $pi * $u2)]
    return [list $dRandNormal1 $dRandNormal2]
}

This returns just one number, without wasting anything. rand() can be replaced with another function of your liking that returns a uniformly distibuted random number greater than or equal to 0 and less than 1.

set pi [expr {acos(-1)}]

set grands {}
proc grand {} {
  global grands pi

  if {[llength $grands] == 0} {
    set u1 [expr {1 - rand()}]
    set u2 [expr {1 - rand()}]
    lappend grands \
      [expr {sqrt(-2 * log($u1)) * cos(2 * $pi * $u2)}] \
      [expr {sqrt(-2 * log($u1)) * sin(2 * $pi * $u2)}]
  }
  set tmp [lindex $grands 0]
  set grands [lrange $grands 1 end]
  return $tmp
}

proc randNormal {mean stdDev} {
  return [expr {$mean + $stdDev * [grand]}]
}

# ----- USAGE -----

puts "John is [randNormal 176.0 9.2] cm tall."
puts "Jane is [randNormal 162.1 9.8] cm tall."