## Version 16 of Poisson distribution

Updated 2007-12-23 16:24:30 by dkf

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 [L1 ] 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```