py23. Poisson Distribution in Python

Given a rate of lambda for an event in some interval, the Poisson distribution finds k such events, in the same interval. We know that the expected value of k has to be lambda.

For a particular k, we have the formula lambda^k/factorial(k) multiplied by a normalization constant. We can rewrite lambda^k/factorial(k) as:

(lambda/k)*(lambda/(k-1))*(lambda/(k-2))*...*(lambda/2)*(lambda) for k not 0, and 1 for k equal to 0

The function get_num finds this value multiplied by the normalization constant, which is exp(-lambda). With the normalization, the total probability for all k is 1.

Besides calculating, we can use the scipy.stats.poisson.pmf function. Since lambda is a keyword in Python, it was written as lamb_da, in the example.

# ex23.py
from __future__ import print_function, division
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import poisson
N = 20
k = np.arange(0,N+1)

def get_num(lamb_da,k):
    coeff = np.exp(-lamb_da)
    if k==0: return coeff
    cv = lamb_da/np.arange(1,k+1)
    return coeff*np.prod(cv)
    

def pois(lamb_da, calc):
    dis_calc = np.array([get_num(lamb_da,i) for i in range(N+1)])
    dis_R = poisson.pmf(k,lamb_da)
    if calc: return dis_calc
    else: return dis_R

c0 = 0.8
pr = 0
st = ['o','v','^','s','d']
plt.plot(k,pois(1, False), st[pr] + '--', color = (0,1-c0,c0))
for i in 3,5,7,9:
    c0 -= .1
    pr = (i-1)//2
    plt.plot(k,pois(i, True), st[pr] + '--', color = (0,1-c0,c0))
    plt.xlabel('k')
    plt.ylabel('p')
    plt.title('Poisson Distribution (lambda = 1,3,5,7,9)')

err = pois(9,True) - pois(9,False)
print('err for lambda of 9 =', err.std())
# err for lambda of 9 = 9.4458398513e-17

Output: