AM (21 june 2007) I came across a relatively simple PRNG recommended by George Marsaglia, arguably an authority on the subject:
# marsaglia.tcl --
# Implementation of a PRNG according to George Marsaglia
#
namespace eval ::PRNG {
variable mod [expr {wide(256)*wide(256)*wide(256)*wide(256)-5}]
variable fac [expr {wide(256)*wide(32)}]
variable x1 [expr {wide($mod*rand())}]
variable x2 [expr {wide($mod*rand())}]
variable x3 [expr {wide($mod*rand())}]
puts $mod
}
proc ::PRNG::marsaglia {} {
variable mod
variable fac
variable x1
variable x2
variable x3
set xn [expr {($fac*($x3+$x2+$x1))%$mod}]
set x1 $x2
set x2 $x3
set x3 $xn
return [expr {$xn/double($mod)}]
}
# main --
# A small test
#
package require Tk
package require Plotchart
pack [canvas .c]
set p [::Plotchart::createXYPlot .c {0 1 0.2} {0 1 0.2}]
set x [::PRNG::marsaglia]
$p dataconfig series -type symbol -symbol dot
for {set i 0} {$i < 1000} {incr i} {
set y [::PRNG::marsaglia]
$p plot series $x $y
set x $y
}The implementation has a few disadvantages: using the builtin random number generator for seeding it means you only get 2**31 different series, but the run-length is about 2**92 for each series.
AM (22 june 2007) I am thinking about a small package/module to make Monte Carlo simulations easy to do.
Here is my first attempt at such simulations:
A simple Monte Carlo simulationAnd here is a first attempt at a package that could be useful:
# random.tcl --
# Create procedures that return various types of pseudo-random
# number generators (PRNGs)
#
# math::random --
# Create the namespace
#
namespace eval ::math::random {
variable count 0
variable pi [expr {4.0*atan(1.0)}]
}
# prng_Binomial --
# Create a PRNG with a binomial distribution
#
# Arguments:
# p Probability that the outcome will be 1
#
# Result:
# Name of a procedure that returns a binomially distributed
# random number
#
proc ::math::random::prng_Binomial {p} {
variable count
incr count
set name ::math::random::PRNG_$count
proc $name {} [string map [list P $p] {return [expr {rand()<P? 1 : 0}]}]
return $name
}
# prng_Uniform --
# Create a PRNG with a uniform distribution in a given range
#
# Arguments:
# min Minimum value
# max Maximum value
#
# Result:
# Name of a procedure that returns a uniformly distributed
# random number
#
proc ::math::random::prng_Uniform {min max} {
variable count
incr count
set name ::math::random::PRNG_$count
set range [expr {$max-$min}]
proc $name {} [string map [list MIN $min RANGE $range] {return [expr {MIN+RANGE*rand()}]}]
return $name
}
# prng_Exponential --
# Create a PRNG with an exponential distribution with given mean
#
# Arguments:
# mean Mean value
#
# Result:
# Name of a procedure that returns a uniformly distributed
# random number
#
proc ::math::random::prng_Uniform {min max} {
variable count
incr count
set name ::math::random::PRNG_$count
proc $name {} [string map [list MEAN $mean] {return [expr {-MEAN*log(rand())}]}]
return $name
}
# prng_Discrete --
# Create a PRNG with a uniform but discrete distribution
#
# Arguments:
# n Outcome is an integer between 0 and n-1
#
# Result:
# Name of a procedure that returns such a random number
#
proc ::math::random::prng_Discrete {n} {
variable count
incr count
set name ::math::random::PRNG_$count
proc $name {} [string map [list N $n] {return [expr {int(N*rand())}]}]
return $name
}
# prng_Normal --
# Create a PRNG with a normal distribution
#
# Arguments:
# mean Mean of the distribution
# stdev Standard deviation of the distribution
#
# Result:
# Name of a procedure that returns such a random number
#
# Note:
# Use the Box-Mueller method to generate a normal random number
#
proc ::math::random::prng_Normal {mean stdev} {
variable count
incr count
set name ::math::random::PRNG_$count
proc $name {} [string map [list MEAN $mean STDEV $stdev] \
{
variable pi
set rad [expr {sqrt(-log(rand()))}]
set phi [expr {2.0*$pi*rand()}]
set r [expr {$rad*cos($phi)}]
return [expr {MEAN + STDEV*$r}]
}]
return $name
}
# prng_Disk --
# Create a PRNG with a uniform distribution of points on a disk
#
# Arguments:
# rad Radius of the disk
#
# Result:
# Name of a procedure that returns the x- and y-coordinates of
# such a random point
#
proc ::math::random::prng_Disk {rad} {
variable count
incr count
set name ::math::random::PRNG_$count
proc $name {} [string map [list RAD $rad] \
{
variable pi
set rad [expr {RAD*sqrt((rand())}]
set phi [expr {2.0*$pi*rand()}]
set x [expr {$rad*cos($phi)}]
set y [expr {$rad*sin($phi)}]
return [list $x $y]
}]
return $name
}
# main --
# Test code
#
set bin [::math::random::prng_Binomial 0.2]
set ones 0
set zeros 0
for { set i 0} {$i < 100000} {incr i} {
if { [$bin] } {
incr ones
} else {
incr zeros
}
}
puts "Binomial: $ones - $zeros"
set discrete [::math::random::prng_Discrete 10]
for { set i 0} {$i < 100000} {incr i} {
set v [$discrete]
if { [info exists count($v)] } {
incr count($v)
} else {
set count($v) 1
}
}
puts "Discrete:"
parray count
set rect [::math::random::prng_Rectangle 10 3]
puts "Rectangle:"
for { set i 0} {$i < 10} {incr i} {
puts [$rect]
}
set normal [::math::random::prng_Normal 0 1]
puts "Normal:"
for { set i 0} {$i < 10} {incr i} {
puts [$normal]
}The idea is that each random variable will be associated with its own command. That way there can be no mixup of distributions or parameters, as there is a command specific for that variable.
