Poisson distribution

poisson.py

# Author: Seiji Emery
# from https://en.wikipedia.org/wiki/Poisson_distribution
#
# Note: uses the obvious recurrence relation
#   P(0) = exp(-lambda)
#   P(k) = P(k - 1) * lambda / k, k > 0
# Which is equivalent to, but much clearer than using factorials imo.

from math import exp

def poisson_p (rate, k = 0):
    ''' Given an average event rate, calculates the discrete probability of k events occuring per unit time '''
    p = exp(-rate)
    for i in range(1, k+1):
        p *= rate / i
    return p

def poisson_k (rate, threshold):
    ''' Given an average event rate and probability threshold on [0, 1], calculates the number of events that would have occured per unit time '''
    p = exp(-rate)
    k = 0
    while threshold > p and k < 100:
        threshold -= p
        k += 1
        p *= rate / k
    return k

if __name__ == '__main__':
    K = 20  # number of counts to use for poisson_p ([0, K])
    N = 20  # number of intervals to use for poisson_k ([0, N], increment of 1 / N)
    rates = [ 0.1, 0.25, 0.5, 1.0, 1.5, 2.5, 5.0, 8.0, 10.0 ] # rates to use / display

    # Display counts
    print("\npoisson_p(rate, k): %s"%(poisson_p.__doc__))
    print("%5s: %s"%("k", ", ".join([ "%5d" % k for k in range(0, K+1) ])))
    for r in rates:
        ps = [ (k, poisson_p(r, k)) for k in range(0, K+1) ]
        print("%5.2f: %s (total %0.4f)"%(
            r, 
            ", ".join([ "%5.3f" % p for _, p in ps ]), 
            sum([ p * k for k, p in ps ])))

    print("\npoisson_k(rate, p): %s"%(poisson_k.__doc__))
    print("%5s: %s"%("p", ", ".join([ "%4.2f" % (float(i) / N) for i in range(N) ])))
    for r in rates:
        ks = [ poisson_k(r, float(i) / N) for i in range(0, N) ]
        print("%5.2f: %s (total %0.4f, approaches %0.4f)"%(
            r, 
            ", ".join([ "%4d" % k for k in ks ]), 
            float(sum(ks)) / N, 
            float(sum([ poisson_k(r, float(i) / 10000) for i in range(10000) ])) / 10000))
    print("")

poisson_output.txt

poisson_p(rate, k):  Given an average event rate, calculates the discrete probability of k events occuring per unit time 
    k:     0,     1,     2,     3,     4,     5,     6,     7,     8,     9,    10,    11,    12,    13,    14,    15,    16,    17,    18,    19,    20
 0.10: 0.905, 0.090, 0.005, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000 (total 0.1000)
 0.25: 0.779, 0.195, 0.024, 0.002, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000 (total 0.2500)
 0.50: 0.607, 0.303, 0.076, 0.013, 0.002, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000 (total 0.5000)
 1.00: 0.368, 0.368, 0.184, 0.061, 0.015, 0.003, 0.001, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000 (total 1.0000)
 1.50: 0.223, 0.335, 0.251, 0.126, 0.047, 0.014, 0.004, 0.001, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000 (total 1.5000)
 2.50: 0.082, 0.205, 0.257, 0.214, 0.134, 0.067, 0.028, 0.010, 0.003, 0.001, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000 (total 2.5000)
 5.00: 0.007, 0.034, 0.084, 0.140, 0.175, 0.175, 0.146, 0.104, 0.065, 0.036, 0.018, 0.008, 0.003, 0.001, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000 (total 5.0000)
 8.00: 0.000, 0.003, 0.011, 0.029, 0.057, 0.092, 0.122, 0.140, 0.140, 0.124, 0.099, 0.072, 0.048, 0.030, 0.017, 0.009, 0.005, 0.002, 0.001, 0.000, 0.000 (total 7.9980)
10.00: 0.000, 0.000, 0.002, 0.008, 0.019, 0.038, 0.063, 0.090, 0.113, 0.125, 0.125, 0.114, 0.095, 0.073, 0.052, 0.035, 0.022, 0.013, 0.007, 0.004, 0.002 (total 9.9655)

poisson_k(rate, p):  Given an average event rate and probability threshold on [0, 1], calculates the number of events that would have occured per unit time 
    p: 0.00, 0.05, 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.45, 0.50, 0.55, 0.60, 0.65, 0.70, 0.75, 0.80, 0.85, 0.90, 0.95
 0.10:    0,    0,    0,    0,    0,    0,    0,    0,    0,    0,    0,    0,    0,    0,    0,    0,    0,    0,    0,    1 (total 0.0500, approaches 0.0998)
 0.25:    0,    0,    0,    0,    0,    0,    0,    0,    0,    0,    0,    0,    0,    0,    0,    0,    1,    1,    1,    1 (total 0.2000, approaches 0.2497)
 0.50:    0,    0,    0,    0,    0,    0,    0,    0,    0,    0,    0,    0,    0,    1,    1,    1,    1,    1,    1,    2 (total 0.4000, approaches 0.4997)
 1.00:    0,    0,    0,    0,    0,    0,    0,    0,    1,    1,    1,    1,    1,    1,    1,    2,    2,    2,    2,    3 (total 0.9000, approaches 0.9996)
 1.50:    0,    0,    0,    0,    0,    1,    1,    1,    1,    1,    1,    1,    2,    2,    2,    2,    2,    3,    3,    4 (total 1.3500, approaches 1.4995)
 2.50:    0,    0,    1,    1,    1,    1,    2,    2,    2,    2,    2,    3,    3,    3,    3,    3,    4,    4,    5,    5 (total 2.3500, approaches 2.4995)
 5.00:    0,    2,    2,    3,    3,    3,    4,    4,    4,    5,    5,    5,    5,    6,    6,    6,    7,    7,    8,    9 (total 4.7000, approaches 4.9991)
 8.00:    0,    4,    5,    5,    6,    6,    6,    7,    7,    7,    8,    8,    9,    9,    9,   10,   10,   11,   12,   13 (total 7.6000, approaches 7.9988)
10.00:    0,    5,    6,    7,    7,    8,    8,    9,    9,    9,   10,   10,   11,   11,   12,   12,   13,   13,   14,   15 (total 9.4500, approaches 9.9986)

[Finished in 0.7s]