Welcome, guest | Sign In | My Account | Store | Cart

This recipe produces random numbers from a normal (bell-shaped) distribution. (The rand() command in Tcl produces uniformly distributed random numbers.)

Tcl, 23 lines
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
# The following procedure generates two independent normally distributed 
# random numbers with mean 0 and vaviance stdDev^2.  If you only need 1
# random number, choose either one.

proc randNormal {stdDev} {
    global dRandNormal1
    global dRandNormal2

    set u1 rand()
    set u2 rand()

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

# The following code is only used to demonstrate the procedure and should
# be replaced by your own code.

for {set i 1} {$i<10} {incr i} {
	randNormal 1
	puts $dRandNormal1
	puts $dRandNormal2
}

This recipe is useful when you need random numbers from a normal distribution, such as simulating measurement noise. This solution is a very simple, straightforward way to convert Tcl's uniformly distributed random numbers to normally distributed random numbers. If you want to use this for Monte Carlo simulation you should pre-process rand() as described at http://mini.net/tcl/1551

7 comments

andreas kupries 19 years, 3 months ago  # | flag

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())}]
William Kappele (author) 19 years, 3 months ago  # | flag

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)

William Kappele (author) 19 years, 3 months ago  # | flag

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

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

to function correctly.

Kregg Kemper 18 years, 2 months ago  # | flag

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

set temp [expr rand()]

set u1 $temp

set u2 $temp

Matej Kristan 17 years, 8 months ago  # | flag

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.

Magnus Åhman 9 years, 10 months ago  # | flag

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]
}
Magnus Åhman 9 years, 10 months ago  # | flag

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."
Created by William Kappele on Tue, 6 Aug 2002 (MIT)
Tcl recipes (162)
William Kappele's recipes (1)

Required Modules

  • (none specified)

Other Information and Tasks