summaryrefslogtreecommitdiffstats
path: root/tcllib/modules/simulation/random.tcl
diff options
context:
space:
mode:
Diffstat (limited to 'tcllib/modules/simulation/random.tcl')
-rwxr-xr-xtcllib/modules/simulation/random.tcl577
1 files changed, 577 insertions, 0 deletions
diff --git a/tcllib/modules/simulation/random.tcl b/tcllib/modules/simulation/random.tcl
new file mode 100755
index 0000000..4541761
--- /dev/null
+++ b/tcllib/modules/simulation/random.tcl
@@ -0,0 +1,577 @@
+# random.tcl --
+# Create procedures that return various types of pseudo-random
+# number generators (PRNGs)
+#
+# Copyright (c) 2007 by Arjen Markus <arjenmarkus@users.sourceforge.net>
+#
+# See the file "license.terms" for information on usage and redistribution
+# of this file, and for a DISCLAIMER OF ALL WARRANTIES.
+#
+# TODO:
+# - Beta
+# - Weighted discrete
+# - Poisson
+# - Cauchy
+# - Binomial
+#
+# Note:
+# Several formulae and algorithms come from "Monte Carlo Simulation"
+# by C. Mooney (Sage Publications, 1997)
+#
+# RCS: @(#) $Id: random.tcl,v 1.5 2012/08/15 04:38:48 arjenmarkus Exp $
+#------------------------------------------------------------------------------
+
+package require Tcl 8.4
+
+# ::simulation::random --
+# Create the namespace
+#
+namespace eval ::simulation::random {
+ variable count 0
+ variable pi [expr {4.0*atan(1.0)}]
+}
+
+
+# prng_Bernoulli --
+# Create a PRNG with a Bernoulli distribution
+#
+# Arguments:
+# p Probability that the outcome will be 1
+#
+# Result:
+# Name of a procedure that returns a Bernoulli-distributed random number
+#
+proc ::simulation::random::prng_Bernoulli {p} {
+ variable count
+
+ incr count
+
+ set name ::simulation::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 ::simulation::random::prng_Uniform {min max} {
+ variable count
+
+ incr count
+
+ set name ::simulation::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:
+# min Minimum value
+# mean Mean value
+#
+# Result:
+# Name of a procedure that returns an exponentially distributed
+# random number
+#
+proc ::simulation::random::prng_Exponential {min mean} {
+ variable count
+
+ incr count
+
+ set name ::simulation::random::PRNG_$count
+ set b [expr {$mean-$min}]
+ proc $name {} [string map [list MIN $min B $b] {return [expr {MIN-B*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 ::simulation::random::prng_Discrete {n} {
+ variable count
+
+ incr count
+
+ set name ::simulation::random::PRNG_$count
+ proc $name {} [string map [list N $n] {return [expr {int(N*rand())}]}]
+
+ return $name
+}
+
+
+# prng_Poisson --
+# Create a PRNG with a Poisson distribution
+#
+# Arguments:
+# lambda The one parameter of the Poisson distribution
+#
+# Result:
+# Name of a procedure that returns such a random number
+#
+proc ::simulation::random::prng_Poisson {lambda} {
+ variable count
+
+ incr count
+
+ set name ::simulation::random::PRNG_$count
+
+ set explambda [expr {exp(-$lambda)}]
+
+ proc $name {} [string map [list INIT $explambda LAMBDA $lambda] {
+ set r [expr {rand()}]
+ set number 0
+ set sum INIT
+ set rfact INIT
+
+ while { $r > $sum } {
+ set rfact [expr {$rfact * LAMBDA /($number+1.0)}]
+ set sum [expr {$sum + $rfact}]
+ incr number
+ }
+ return $number
+ }]
+
+ 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 ::simulation::random::prng_Normal {mean stdev} {
+ variable count
+
+ incr count
+
+ set name ::simulation::random::PRNG_$count
+ proc $name {} [string map [list MEAN $mean STDEV $stdev] \
+ {
+ variable pi
+ set rad [expr {sqrt(-2.0*log(rand()))}]
+ set phi [expr {2.0*$pi*rand()}]
+ set r [expr {$rad*cos($phi)}]
+ return [expr {MEAN + STDEV*$r}]
+ }]
+
+ return $name
+}
+
+
+# prng_Pareto --
+# Create a PRNG with a Pareto distribution
+#
+# Arguments:
+# min Minimum value for the distribution
+# steep Steepness of the descent
+#
+# Result:
+# Name of a procedure that returns a Pareto-distributed number
+#
+proc ::simulation::random::prng_Pareto {min steep} {
+ variable count
+
+ incr count
+
+ set name ::simulation::random::PRNG_$count
+
+ set rsteep [expr {1.0/$steep}]
+ proc $name {} [string map [list MIN $min RSTEEP $rsteep] \
+ {
+ return [expr {MIN * pow(1.0-rand(),RSTEEP)}]
+ }]
+
+ return $name
+}
+
+
+# prng_Gumbel --
+# Create a PRNG with a Gumbel distribution
+#
+# Arguments:
+# min Minimum value for the distribution
+# f Factor to scale the value
+#
+# Result:
+# Name of a procedure that returns a Gumbel-distributed number
+#
+# Note:
+# The chance P(v) = exp( -exp( f*(v-min) ) )
+#
+proc ::simulation::random::prng_Gumbel {min f} {
+ variable count
+
+ incr count
+
+ set name ::simulation::random::PRNG_$count
+
+ proc $name {} [string map [list MIN $min F $f] \
+ {
+ return [expr {MIN + log( -log(1.0-rand()) ) / F}]
+ }]
+
+ return $name
+}
+
+
+# prng_chiSquared --
+# Create a PRNG with a chi-squared distribution
+#
+# Arguments:
+# df Degrees of freedom
+#
+# Result:
+# Name of a procedure that returns a chi-squared distributed number
+# with mean 0 and standard deviation 1
+#
+proc ::simulation::random::prng_chiSquared {df} {
+ variable count
+
+ incr count
+
+ set name ::simulation::random::PRNG_$count
+ proc $name {} [string map [list DF $df] \
+ {
+ variable pi
+ set y 0.0
+ for { set i 0 } { $i < DF } { incr i } {
+ set rad [expr {sqrt(-log(rand()))}]
+ set phi [expr {2.0*$pi*rand()}]
+ set r [expr {$rad*cos($phi)}]
+ set y [expr {$y+$r*$r}]
+ }
+ return [expr {($y-DF)/sqrt(2.0*DF)}]
+ }]
+
+ 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 ::simulation::random::prng_Disk {rad} {
+ variable count
+
+ incr count
+
+ set name ::simulation::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
+}
+
+
+# prng_Ball --
+# Create a PRNG with a uniform distribution of points within a ball
+#
+# Arguments:
+# rad Radius of the ball
+#
+# Result:
+# Name of a procedure that returns the x-, y- and z-coordinates of
+# such a random point
+#
+proc ::simulation::random::prng_Ball {rad} {
+ variable count
+
+ incr count
+
+ set name ::simulation::random::PRNG_$count
+ proc $name {} [string map [list RAD $rad] \
+ {
+ variable pi
+ set rad [expr {RAD*pow(rand(),0.333333333333)}]
+ set phi [expr {2.0*$pi*rand()}]
+ set theta [expr {acos(2.0*rand()-1.0)}]
+ set x [expr {$rad*cos($phi)*cos($theta)}]
+ set y [expr {$rad*sin($phi)*cos($theta)}]
+ set z [expr {$rad*sin($theta)}]
+ return [list $x $y $z]
+ }]
+
+ return $name
+}
+
+
+# prng_Sphere --
+# Create a PRNG with a uniform distribution of points on the surface
+# of a sphere
+#
+# Arguments:
+# rad Radius of the sphere
+#
+# Result:
+# Name of a procedure that returns the x-, y- and z-coordinates of
+# such a random point
+#
+proc ::simulation::random::prng_Sphere {rad} {
+ variable count
+
+ incr count
+
+ set name ::simulation::random::PRNG_$count
+ proc $name {} [string map [list RAD $rad] \
+ {
+ variable pi
+ set phi [expr {2.0*$pi*rand()}]
+ set theta [expr {acos(2.0*rand()-1.0)}]
+ set x [expr {RAD*cos($phi)*cos($theta)}]
+ set y [expr {RAD*sin($phi)*cos($theta)}]
+ set z [expr {RAD*sin($theta)}]
+ return [list $x $y $z]
+ }]
+
+ return $name
+}
+
+
+# prng_Rectangle --
+# Create a PRNG with a uniform distribution of points in a rectangle
+#
+# Arguments:
+# length Length of the rectangle (x-direction)
+# width Width of the rectangle (y-direction)
+#
+# Result:
+# Name of a procedure that returns the x- and y-coordinates of
+# such a random point
+#
+proc ::simulation::random::prng_Rectangle {length width} {
+ variable count
+
+ incr count
+
+ set name ::simulation::random::PRNG_$count
+
+ proc $name {} [string map [list LENGTH $length WIDTH $width] \
+ {
+ set x [expr {LENGTH*rand()}]
+ set y [expr {WIDTH*rand()}]
+ return [list $x $y]
+ }]
+
+ return $name
+}
+
+
+# prng_Block --
+# Create a PRNG with a uniform distribution of points in a block
+#
+# Arguments:
+# length Length of the block (x-direction)
+# width Width of the block (y-direction)
+# depth Depth of the block (y-direction)
+#
+# Result:
+# Name of a procedure that returns the x-, y- and z-coordinates of
+# such a random point
+#
+proc ::simulation::random::prng_Block {length width depth} {
+ variable count
+
+ incr count
+
+ set name ::simulation::random::PRNG_$count
+
+ proc $name {} [string map [list LENGTH $length WIDTH $width DEPTH $depth] \
+ {
+ set x [expr {LENGTH*rand()}]
+ set y [expr {WIDTH*rand()}]
+ set z [expr {DEPTH*rand()}]
+ return [list $x $y $z]
+ }]
+
+ return $name
+}
+
+# Announce the package
+#
+package provide simulation::random 0.3.1
+
+
+# main --
+# Test code
+#
+if { 0 } {
+set bin [::simulation::random::prng_Bernoulli 0.2]
+
+set ones 0
+set zeros 0
+for { set i 0} {$i < 100000} {incr i} {
+ if { [$bin] } {
+ incr ones
+ } else {
+ incr zeros
+ }
+}
+
+puts "Bernoulli: $ones - $zeros"
+
+set discrete [::simulation::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 [::simulation::random::prng_Rectangle 10 3]
+
+puts "Rectangle:"
+for { set i 0} {$i < 10} {incr i} {
+ puts [$rect]
+}
+
+set normal [::simulation::random::prng_Normal 0 1]
+
+puts "Normal:"
+for { set i 0} {$i < 10} {incr i} {
+ puts [$normal]
+}
+
+#
+# Timing: how fast is the normal random number generator?
+#
+# Surprising speed: 15 million numbers per minute!
+#
+puts "Normal random number generator:"
+puts "[time {set value [$normal]} 30000]"
+set result [lindex [time {set value [$normal]} 30000] 0]
+puts "[expr {60.0e6/$result}] numbers per minute"
+puts "Creating a long list: [time {lappend value [$normal]} 30000]"
+puts "[lrange $value 0 20] - [llength $value] numbers in total"
+set value {}
+set result [lindex [time {lappend value [$normal]} 30000] 0]
+puts "[expr {60.0e6/$result}] numbers per minute"
+
+puts "Points in a rectangle:"
+puts "[time {set value [$rect]} 30000]"
+set result [lindex [time {set value [$rect]} 30000] 0]
+puts "[expr {60.0e6/$result}] numbers per minute"
+
+#
+# A more formal test
+#
+package require math
+unset count
+set samples 100000
+
+set lambda 10.0
+set poisson [::simulation::random::prng_Poisson $lambda]
+for { set i 0 } { $i < $samples } { incr i } {
+ set number [$poisson]
+ if { [info exists count($number)] } {
+ incr count($number)
+ } else {
+ set count($number) 1
+ }
+}
+parray count
+for { set i 0 } { $i < 30 } { incr i } {
+ set expected [expr {int($samples * pow($lambda,$i) * exp(-$lambda) / [::math::factorial $i])}]
+ set exp_error [expr {sqrt($expected)}]
+ if { [info exists count($i)] } {
+ if { $expected-$exp_error < $count($i) &&
+ $expected+$exp_error > $count($i) } {
+ set okay "okay"
+ } else {
+ set okay "difference too large"
+ }
+ puts "$i $expected $count($i) - [expr {$expected/double($count($i))}] - $okay"
+ } else {
+ puts "$i $expected none"
+ }
+}
+}
+
+#
+# Test hypothesis concerning rectangle
+#
+if { 0 } {
+set r2 [::simulation::random::prng_Rectangle2 10 1]
+
+set count_down 0
+set count_up 0
+set count_left 0
+set count_right 0
+for { set i 0 } { $i < 1000000 } { incr i } {
+
+ foreach {x y} [$r2] {
+ if { $x < 2.0 } { incr count_left }
+ if { $x > 8.0 } { incr count_right }
+ if { $y < 0.2 } { incr count_down }
+ if { $y > 0.8 } { incr count_up }
+ }
+}
+puts "Left-right:\t$count_left\t$count_right"
+puts "Up-down: \t$count_up\t$count_down"
+}
+
+#
+# Check normal distribution
+#
+if { 0 } {
+ package require math::statistics
+ set normal [::simulation::random::prng_Normal 0 1]
+
+ for { set i 0} {$i < 1000} {incr i} {
+ lappend numbers [$normal]
+ }
+ puts "Mean: [::math::statistics::mean $numbers]"
+ puts "Stdev: [::math::statistics::stdev $numbers]"
+}
generated by cgit v0.12 at 2024-12-27 06:32:43 (GMT)