Version 5 of Poisson distribution

Updated 2007-12-23 11:51:27 by EKB

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). What is the probability, that the Wiki will have 20 visitors a day? The answer is calculated using the Poisson distribution:

         20    -12
       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

See also: Statistical Distributions


Category Mathematics | Category Statistics