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()}}}]Lars H: I claim in Random Number that shuffling such as the above doesn't do much for the period of the PRNG. The following is the reasoning that led me to that conclusion.
Not all of the data present in the numbers being shuffled are part of the practical state of the shuffler. Concretely, it is possible to view the numbers merely as indices into the array of numbers and still operate the shuffler exactly as above. To illustrate this follows here a modification of the first implementation which keeps the values and indices as separate arrays.
namespace eval BayesDurham {
namespace export rand
variable poolSize 97 ;# Size of the pool of random numbers
variable value ;# Pool of random numbers
variable index ;# Corresponding index values
variable y ;# Index of last output
}
proc BayesDurham::num2index {num} {
variable poolSize
return [expr { int( double($poolSize) * $num ) }]
}
proc BayesDurham::initFromRand {} {
variable poolSize
variable value
variable index
variable y
# 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 value($i) [expr rand()]
set index($i) [num2index $value($i)]
}
# Generate "last" index
set y [num2index [expr rand()]]
}
proc BayesDurham::initFromRandIfNecessary {} {
initFromRand
# Turn self into no-op.
proc initFromRandIfNecessary {} {}
}
proc BayesDurham::rand {} {
variable value
variable index
variable y
initFromRandIfNecessary
set i $y
set y $index($i)
set r $value($i)
set value($i) [expr rand()]
set index($i) [num2index $value($i)]
return $r
}Here, the value array has no influence on how the numbers are shuffled; everything is controlled by the index array and the y variable. Hence the period exhibited by BayesDurham::rand is no greater than that exhibited by
proc BayesDurham::step {} {
variable index
variable y
initFromRandIfNecessary
set i $y
set y $index($i)
set index($i) [num2index [expr rand()]]
return $y
}The poolSize entries in index take values in the range 0 through $poolSize-1, and y also takes values in this range, so the shuffling extends the PRNG state space relevant for the period by a factor pow($poolSize, $poolSize+1). It doesn't extend the period by that much, however.
For the period to grow by the same factor as the state space when adding shuffling, it would be necessary that the effect on y and the index array of
for {set n [periodOfRand]} {$n>0} {incr n -1} {BayesDurham::step}(running through one period of rand, so that we return that part of the PRNG to its original state) preserves all the information stored in these variables, but that is not the case; this operation is not a bijection. Indeed, it is generally not possible to run BayesDurham::step backwards even one step — to deduce from post-state of y and index, and knowing the entire sequence of numbers generated by rand — what the pre-state of these variables were. Suppose [num2index [expr rand()]]] was j in this call. Then some index($i) was set to j, and the old value of index($i) was stored in y, but the old value of y (which, incidentally, is the i) is lost, so the process cannot be reversed. It is possible to determine this old value i of y if only one entry in index is equal to j, but as soon as there is more than one such entry, if could have been either of them; two distinct pre-states can lead to the same post-state, and since we don't know which of the two it was, we've lost 1 bit of information.
[Monowidth formulae below are in TeX notation.] Generally, if there are k entries equal to j in index then we lose \log_2 k bits of information, but how much is lost on average depends on how probable the various arrangements are. The total number of length n lists of integers from 0 to n–1 inclusive is n^n. The number of such lists where exactly k entries are equal to some given j in this range is \binom{n}{k} (n-1)^{n-k} — this corresponds to an expansion of the binomial ((n-1) + 1)^n — and we're not interested in the case k=0 since j is always going to be a number that is actually in the list. Hence there are n^n - (n-1)^n possible values for index, and assuming all of them to be of equal probability (they are after all generated by filling the array with random numbers), we find that the probability of having j appearing exactly k times in index is
\frac{
\binom{n}{k} (n-1)^{n-k}
}{
n^n - (n-1)^n
} =
\frac{ (n-1)^n }{ n^n - (n-1)^n } \binom{n}{k} (n-1)^{-k} =
( (\frac{n}{n-1})^n - 1 )^{-1} \prod_{l=1}^k \frac{ n-l+1 }{ l (n-1) }The expected number of bits lost at each iteration of BayesDurham::step is thus
\frac{ (n-1)^n }{ n^n - (n-1)^n } \sum_{k=1}^n \binom{n}{k} (n-1)^{-k} \log_2 kwhich turns out to be slightly above 0.5 for all n>=42. A Horner-style unrolling of this sum gives the procedure
proc E_bitloss {n} {
set sum 0.0
for {set k $n} {$k >= 1} {incr k -1} {
set sum [expr {$sum * ($n-$k)/($k+1)/($n-1) + log($k)}]
}
expr { $sum/log(2) * $n/($n-1.0) / ( pow( 1-1.0/$n, -$n) - 1) }
}With a poolSize of 97, there are (slightly less than) 647 bits of information in the shuffler, so after 1294 steps one expects all of them to have been lost. Since rand has a period of 2147483646, the odds for the shuffler of keeping some information around an entire period of rand seems astronomically small…
Still, that's a purely probabilistic analysis. It should work fine at the beginning, but the more steps being run, the lower is the chance that it accurately describes what is happening. My guess would be that there are small sets of initial values for which the trajectories make it around an entire rand period without colliding, but finding even a pair of such values is not easy, and then comes the matter of what happens on the next rand period. Things have to work out very nicely to even make the shuffler exhibit a period twice that of the underlying PRNG; the typical arrangement would be a period exactly equal to that of the underlying PRNG.
For a sufficiently good underlying PRNG, it can actually be proved that there is a point along the period where the shuffler is reset to a particular state regardless of how it was initialised; from this follows that the period of the shuffled PRNG is exactly equal to that of the underlying PRNG. The idea is merely to construct a sequence of new indices that will set the shuffler to a given state, since there is a positive (although very small) probability that this particular sequence of numbers will appear on random, and any sufficiently good PRNG will therefore generate it somewhere along its period.
One useful sequence of indices for resetting the shuffler is a long sequence of all zeroes. (Remember, since we're talking indices here, this doesn't necessarily mean that rand outputs 0.0, only that rand*$poolSize < 1.) If y and index(0) are both zero, then a step in which the new index is also a zero won't change the state of the shuffler. Once in such a state, it is possible to set any individual entry in index to zero (while keeping all entries already zero unchanged) by introducing its index j into a sequence of mostly zeroes, as follows:
| step | new index | y | index(0) | index($j) |
|---|---|---|---|---|
| 0 | j | 0 | 0 | ? |
| 1 | 0 | 0 | j | ? |
| 2 | 0 | j | 0 | ? |
| 3 | 0 | ? | 0 | 0 |
Finally, in a state with unknown value of y and at most k nonzero entries in index, a sequence of k+1 zeroes as incoming new indices will ensure that y=0 and index(0)=0, again while keeping all entries already zero unchanged; the only way for y to stay nonzero is that it picks up a nonzero value from some entry in index, but doing so will increase the number of entries that are zero, so at the very latest after k+1 zeroes, there won't be any left to pick and y=0. (More probable is however that a zero is picked up much earlier, and in that case the number of nonzero elements isn't reduced all that much.)
Combining the above, one finds that a suitable constructed sequence (mostly zeroes) of n^2/2 + O(n) new indices will reset a shuffler with poolSize=n to an all zero state. For n=97, and ignoring the linear term (I just want to get the order of magnitude for the exponent right), this is a sequence of about 4700 random numbers, with an information content of about 31050 bits. In order for a PRNG to reproduce that particular sequence, it needs to reproduce all bit sequences of that length, so the period has to be somewhere in the vicinity of, say, 2^{31065}. The most common Mersenne Twister has a period of 2^{19937}-1, so it is not sufficient, but OTOH the 27th Mersenne prime is 2^{44497}-1, so a twister based on this would suffice with great margin.
Then again, even 2^{19937}-1 (roughly 4.3e6001) is so large that the chances of actually running any computation that many steps is very low (if protons are unstable [L1 ], the entire computer would probably have decayed into radiation long before the time at which the first period is expected to be complete). In these scales, the period is probably not very relevant anyway.