Purpose: Show how to improve the output of the math function rand to reduce sequential correlation.
The sequence of numbers returned by rand is prone to sequential correlations if large runs of numbers are taken. For example, if clusters of six sequential results from
expr { rand () }are analyzed, it will turn out that they lie on only 40 hyperplanes in six-dimensional space. This limitation means that rand() is not suitable by itself for applications like high-resolution Monte Carlo integration.
The following code implements a "better" rand that reduces this problem by shuffling the sequence that is returned. For further discussion, see Donald E. Knuth, The Art of Computer Programming, Volume 2: Seminumerical Algorithms, which has a wealth of information about random number generators and their foibles.
namespace eval BayesDurham {
namespace export rand
variable poolSize 97; # Size of the pool of random numbers
variable v; # Pool of random numbers
variable y; # Most recently generated random number
}
#----------------------------------------------------------------------
#
# BayesDurham::rand --
#
# Generate a random floating point number between 0 and 1
#
# Parameters:
# None.
#
# Results:
# Returns a randomly chosen floating point number in the
# range [0 .. 1)
#
# Side effects:
# Stores state in several namespace variables.
#
# The BayesDurham::rand procedure attempts to improve upon the
# system-supplied rand() function by breaking up sequential
# correlations. It works by shuffling a vector of numbers supplied by
# the system's random number generator, using its own output to choose
# the next slot in the vector.
#
#----------------------------------------------------------------------
proc BayesDurham::rand {} {
variable poolSize
variable v
variable y
# Test whether we're initializing the sequence
if { ![info exists v] } {
# Discard some random numbers - in case the RNG was just
# seeded.
for { set i 0 } { $i < $poolSize } { incr i } {
expr { rand() }
}
# Save a pool of random numbers
for { set i 0 } { $i < $poolSize } { incr i } {
set v($i) [expr { rand() }]
}
# Generate one more random number
set y [expr { rand() }]
}
# Use the most recently generated number to select from the pool.
set i [expr { int( double( $poolSize ) * $y ) }]
set y $v($i)
set v($i) [expr { rand() }]
return $y
}EKB Here's an alternative implementation:
namespace eval shufflerand {
variable vals
variable vallen 98
}
proc shufflerand::init {args} {
variable vals
variable vallen
set haveseed false
foreach {arg argval} $args {
switch -- $arg {
-len {
if [string is integer $argval] {
if {$argval > 0} {
set vallen $argval
} else {
error "[namespace current]::init: Length must be greater than zero"
}
} else {
error "[namespace current]::init: Length must be an integer"
}
}
-seed {
if [string is integer $argval] {
set seed $argval
set haveseed true
} else {
error "[namespace current]::init: Seed must be an integer"
}
}
default {
error "[namespace current]::init: Possible arguments are: -seed -len"
}
}
}
set vals {}
if $haveseed {
lappend vals [expr {srand($seed)}]
} else {
lappend vals [expr {rand()}]
}
for {set i 1} {$i < $vallen} {incr i} {
lappend vals [expr {rand()}]
}
}
proc shufflerand::rand {} {
variable vals
variable vallen
set retval [lindex $vals end]
set ndx [expr {int($vallen * $retval)}]
lset vals end [lindex $vals $ndx]
lset vals $ndx [expr {rand()}]
return $retval
}
# Test
console show
shufflerand::init
puts "A list of 10 random numbers from shuffled list:"
for {set i 0} {$i < 10} {incr i} {
puts [shufflerand::rand]
}
puts "\nTime to generate 1000 values using rand():"
puts [time {for {set i 0} {$i < 1000} {incr i} {expr {rand()}}}]
puts "\nTime to generate 1000 values using shufflerand:"
puts [time {for {set i 0} {$i < 1000} {incr i} {shufflerand::rand}}]Here's some sample output from the test code at the end:
A list of 10 random numbers from shuffled list: 0.00284572923688 0.897247713477 0.649551953492 0.456922465682 0.0448791366279 0.0441727326457 0.232282039818 0.645288598093 0.753969432672 0.207090688034 Time to generate 1000 values using rand(): 558 microseconds per iteration Time to generate 1000 values using shufflerand: 4052 microseconds per iteration
So, it's about 7-8 times slower than the straight linear congruential generator that's bundled with Tcl.
EKB Hm...Now using Tcl 8.5, rand() takes a little longer, shufflerand a little less long. Here are some runs for rand(), shufflerand, and Bayes-Durham:
Time to generate 1000 values using rand(): 738 microseconds per iteration Time to generate 1000 values using shufflerand: 3202 microseconds per iteration Time to generate 1000 values using Bayes-Durham rand: 3941 microseconds per iteration
So slightly longer with Bayes-Durham (about 20%). Possibly because of the test for whether the pool had been created or not. I decided not to include that in the "shufflerand" proc because typically the call to rand() is in the heart of a loop. I also cached the pool size.
EKB This is a version for Tcl 8.5 that replaces the existing rand() with a shuffled version, so code that uses rand() doesn't have to be rewritten:
namespace eval shufflerand {
variable vals
}
proc shufflerand::init {} {
variable vals
for {set i 0} {$i < 98} {incr i} {
lappend vals [expr {rand()}]
}
# lcrand -- "Linear Congruential rand"
rename ::tcl::mathfunc::rand ::tcl::mathfunc::lcrand
proc ::tcl::mathfunc::rand {} "[namespace current]::rand"
}
proc shufflerand::rand {} {
variable vals
set retval [lindex $vals end]
set ndx [expr {int(98 * $retval)}]
lset vals end [lindex $vals $ndx]
lset vals $ndx [expr {lcrand()}]
return $retval
}
# Test
console show
puts "\nTime to generate 1000 values using default rand():"
puts [time {for {set i 0} {$i < 1000} {incr i} {expr {rand()}}}]
shufflerand::init
puts "\nTime to generate 1000 values using shuffled rand():"
puts [time {for {set i 0} {$i < 1000} {incr i} {expr {rand()}}}]