where random numbers flop

Testing random numbers, as in simulating throwing a dice:


proc roll {} {
        set seed [expr {[clock clicks] + [pid] + [clock seconds]}]
        set seed [expr {$seed - ($seed * 2)}]
        #set seed [expr sqrt($seed)]
        expr (srand ($seed)) 
        set seed 0
        set number [expr {rand() * 6}]
        return [set number [expr {round($number)}]]
}

set ::a 0
set ::1 0
set ::2 0
set ::3 0
set ::4 0
set ::5 0
set ::6 0 

proc statistics {} {
        set a 0
        while {$a != 100000} {
                set roll [roll]

                while {$roll == 0} {
                        set roll [roll]
                }

                foreach value {1 2 3 4 5 6} {
                        if {$value == $roll} {
                                incr ::$value
                        }
                }
                incr a
        } 
        set a 0

        puts ">1= $::1" 
        puts ">2= $::2" 
        puts ">3= $::3" 
        puts ">4= $::4" 
        puts ">5= $::5" 
        puts ">6= $::6"

        set ::1 0
        set ::2 0
        set ::3 0
        set ::4 0
        set ::5 0
        set ::6 0
}

This snippet of code should produce random numbers in the range of 1-6, but if you apply a large enough sample, you notice that the output becomes rather predictable. The generator tends to favor 3 and 4 in this instance, and favors 1 rather than 6. This might have something to do with the way round() handles its math, but given the source code above, theres a rather large degree of predictability going on. If I were a gambler Id have my money on 3 everytime.

 100000
 >1= 17471
 >2= 6785
 >3= 25872
 >4= 24633
 >5= 19663
 >6= 5576

 100000
 >1= 17226
 >2= 7011
 >3= 25862
 >4= 25001
 >5= 19008
 >6= 5892

 100000
 >1= 17361
 >2= 7183
 >3= 24947
 >4= 24568
 >5= 19996
 >6= 5945

originally created by Davou


AM 2007-02-15: Why do you seed and reseed rand() every time? That is not the proper way to use rand(). I think much of the non-randomness you are experiencing could be attributed to that (the seeds tend to be correlated!)

MS Notes that there is one more problem with that code: you should not use round(), use round(ceil()) instead. The reason is that

   expr {6*rand()}

will produce numbers uniformly distributed in (0,6). Rounding them, you get

   (0,0.5)   -> 0
   (0.5,1.5) -> 1
   (1.5,2.5) -> 2
   ...
   (4.5,5.5) -> 5
   (5.5,6)   -> 6

In other words: you would get a whole lot of zeros (which you are not counting), and half as many sixes as the other numbers.

So: do not reseed at each call (as noted by AM), and use instead

proc roll {} {
    expr {round(ceil(rand()*6))}
}

and see the stats behave rather well.


CJL couldn't resist reducing the code, while maintaining the original structural theme so that the original author can appreciate the differences:

proc roll {} {
    expr {1+int(6*rand())}
}

proc statistics {} {
    set 1 0
    set 2 0
    set 3 0
    set 4 0
    set 5 0
    set 6 0
  
    set a 0
  
    while {$a != 100000} {
        incr [roll]
        incr a
    }
  
    puts ">1= $1"
    puts ">2= $2"
    puts ">3= $3"
    puts ">4= $4"
    puts ">5= $5"
    puts ">6= $6"
}

A few inefficiencies that have been removed:

  • As srand($seed) returns the first number in the series for $seed, that value could have been used rather than the additional call to rand(). But as specifying our own seed for every call to roll is more likely to reduce randomness than improve it, I've left it out altogether.
  • Generating numbers in the range 0..6 and throwing away the zeros
  • Looping over all possible values to increment just one of them
  • Setting a local variable back to zero just before exiting a procedure

slebetman: For the sake of completeness here's a sampling of the results of the new code above:

  % statistics
  >1= 16553
  >2= 16593
  >3= 16951
  >4= 16518
  >5= 16544
  >6= 16841
  % statistics
  >1= 16783
  >2= 16536
  >3= 16731
  >4= 16691
  >5= 16597
  >6= 16662
  % statistics
  >1= 16496
  >2= 16813
  >3= 16697
  >4= 16700
  >5= 16601
  >6= 16693

Looks good to me.


EB: you get better distribution with rand() when multiplying by a large value and then using modulo:

set n [expr {round(rand()*65535)%6+1}]