This is a discrete distribution of a random variable. An example may tell us what it can be used for:
Suppose this Wiki is visited by 12 people daily (as an average ''- [DKF]: very far on the low side, btw''). What is the probability, that the Wiki will have 20 visitors a day? The answer is calculated using the Poisson distribution:
!!!!!!
12²⁰ * e⁻¹² <
> p = ------------ <
> 20!
!!!!!!
The answer is: p = 0.0097 which is 0.97%. Another example is: [Simulating a server system].
Here is an implementation that draws integer numbers being Poisson distributed. It is based on an ecology modelling book I am currently reading and this thread on c.l.t [http://groups.google.com/group/comp.lang.tcl/browse_thread/thread/bdfb44464fe80b5c] helped finding a good Tcl implementation. That's a discussion on how to draw random numbers from a given probability distribution (not just the onw shown here).
proc RandomPoisson {lambda count} {
#
# generate random numbers that are poisson distributed
#
# lambda -> expected value = "mean value" (which happens to also be the variance)
# count -> number of numbers to be generated
#
# (adapted from: "Parameter Estimation in Ecology" by O. Richter & D. Söndgerath)
# factorial f (here: the factorial of 0):
set f 1
# poisson probability of the value 0:
set su [expr {exp(-$lambda)}]
set p(0) $su
# probabilities of the integers up to the math limits of Tcl:
# (computed as the discrete density function)
set i 0
while {1} {
incr i
set f [expr {$f*$i}]
set su [expr {$su + exp(-$lambda) * pow($lambda,$i)/double($f)}]
if {$su > 1} {
# we cannot calculate more precisely here,
# so we assume 1 is ok for the final density:
set p($i) 1
break
}
set p($i) $su
}
# calculate random values according to the
# given discrete probability density function:
# (this does work for all dicrete distributions,
# not only for poisson)
for {set c 0} {$c < $count} {incr c} {
# random number in the interval [0,1]:
set x [expr {rand()}]
# transform this number to the correct target interval:
for {set j 0} {$j <= $i} {incr j} {
if {$p($j) > $x} {
lappend result $j
break
}
}
}
return $result
}
We run into a problem here when using Tcl's normal [expr]. The numbers in the formula quickly put [expr] to its limits and therefor I have added the case where the probability is just set to 1, when this point arrives. You could use [mpexpr] instead if you need a better approximation of the very improbable part of the Poisson distribution.
Here is a test:
# produce 10,000 random integer numbers according to the Poisson distribution:
set data [RandomPoisson 3.5 10000]
# count their frequencies:
foreach el $data {
if {[info exists count($el)]} {
incr count($el) 1
} else {
set count($el) 1
}
}
# and display their frequencies:
foreach el [lsort -integer [array names count]] {
puts "$el => $count($el)"
}
----
[EKB] I just found this page after doing an implementation of my own, based on a book I'm reading on simulation & Monte Carlo techniques. I don't think my code is preferable to this, but just for the sake of variety, I'm pasting it here. It also generates a little test code that shows the pdf overlaying a sample of random deviates, to see how they compare.
Note that this does ''not'' work well for large values of mu. Faster algorithms are available for larger mu.
# Generate a poisson-distributed random deviate
# Use algorithm in section 4.9 of Dagpunar, J.S,
# "Simulation and Monte Carlo: With Applications
# in Finance and MCMC", pub. 2007 by Wiley
proc random-poisson {mu number} {
set W0 [expr {exp(-$mu)}]
set retval {}
for {set i 0} {$i < $number} {incr i} {
set W $W0
set R [expr {rand()}]
set X 0
while {$R > $W} {
set R [expr {$R - $W}]
incr X
set W [expr {$W * $mu/double($X)}]
}
lappend retval $X
}
return $retval
}
proc pdf-poisson {mu x} {
# To avoid dividing large values, build up this way:
set W [expr {exp(-$mu)}]
for {set z 1} {$z <= int($x)} {incr z} {
set W [expr {$W * $mu/double($z)}]
}
return $W
}
proc cdf-poission {mu x} {
# To avoid dividing large values, build up:
set W [expr {exp(-$mu)}]
set retval $W
for {set z 1} {$z <= int($x)} {incr z} {
set W [expr {$W * $mu/double($z)}]
set retval [expr {$retval + $W}]
}
return $retval
}
##########################################################################################
##
## TESTING
##
## Can test pdf & cdf by running in a console. For random numbers, generate histograms:
##
##########################################################################################
package require math::statistics
canvas .c
pack .c -side top
frame .f
pack .f -side bottom
label .f.mul -text "mu"
entry .f.mue -textvariable mu
pack .f.mul -side left
pack .f.mue -side left
button .f.run -text "Run" -command runtest
pack .f.run -side left
proc runtest {} {
set vals [random-poisson $::mu 5000]
set remainder 5000
for {set x 0.0} {$x < 20.0} {set x [expr {$x + 1}]} {
lappend bins $x
set distval [pdf-poisson $::mu $x]
set distval [expr {int(5000 * $distval)}]
lappend distcounts $distval
}
# Assume none are left
lappend distcounts 0.0
set bincounts [::math::statistics::histogram $bins $vals]
.c delete hist
.c delete dist
math::statistics::plot-scale .c 0 20 0 [math::statistics::max $bincounts]
math::statistics::plot-histogram .c $bincounts $bins hist
math::statistics::plot-histogram .c $distcounts $bins dist
.c itemconfigure dist -fill {} -outline green
}
set mu 5
runtest
----
See also: [Statistical Distributions]
----
[Category Mathematics] | [Category Statistics]