In [1]:
# Goal: approximate true Poisson as limit of Binomial distributions
In [2]:
# simulate binomial distribution
def binom(n,p):
    """Sample from Binomial(n,p) distribution. Number of heads
    when flip n coins (each with prob p of landing heads)."""
    run_sum = 0 # running total of number of heads so far
    for i in range(n):
        if random()<p: # random() generates uniform r.v. in (0,1)
            # get here with probability p 
            run_sum += 1 # increase count of heads so far
    return run_sum
In [15]:
import matplotlib.pyplot as plt
plt.hist([binom(1,1/2) for i in range(1000)])
Out[15]:
(array([513.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0., 487.]),
 array([0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1. ]),
 <a list of 10 Patch objects>)
In [18]:
import matplotlib.pyplot as plt
plt.hist([binom(3,1/4) for i in range(1000)])
Out[18]:
(array([423.,   0.,   0., 425.,   0.,   0., 143.,   0.,   0.,   9.]),
 array([0. , 0.3, 0.6, 0.9, 1.2, 1.5, 1.8, 2.1, 2.4, 2.7, 3. ]),
 <a list of 10 Patch objects>)
In [21]:
import matplotlib.pyplot as plt
plt.hist([binom(100,1/4) for i in range(10000)], bins=range(10,45), align="left", rwidth=0.5)
Out[21]:
(array([  2.,   1.,   6.,  11.,  28.,  52., 108., 163., 261., 376., 478.,
        593., 768., 846., 944., 899., 838., 808., 713., 587., 441., 362.,
        249., 167., 107.,  84.,  44.,  29.,  15.,   9.,   7.,   1.,   3.,
          0.]),
 array([10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26,
        27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43,
        44]),
 <a list of 34 Patch objects>)
In [22]:
# if we take binom(n,p), for some fixed p, as n gets large converges to
# normal distribution (Central Limit Theorem)
In [23]:
# Get Poisson(lam) distribution by taking limit of binom(n,lam/n), 
# as n-> infinity
In [26]:
lam=3.0 # lambda, intensity of Poisson distribution
n=100
plt.hist([binom(n,lam/n) for i in range(10000)], bins=range(0,15), align="left", rwidth=0.5)
Out[26]:
(array([4.640e+02, 1.512e+03, 2.190e+03, 2.258e+03, 1.733e+03, 1.053e+03,
        4.860e+02, 1.940e+02, 7.300e+01, 2.500e+01, 1.000e+01, 2.000e+00,
        0.000e+00, 0.000e+00]),
 array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14]),
 <a list of 14 Patch objects>)
In [27]:
lam=3.0 # lambda, intensity of Poisson distribution
n=1000
plt.hist([binom(n,lam/n) for i in range(10000)], bins=range(0,15), align="left", rwidth=0.5)
Out[27]:
(array([ 494., 1491., 2219., 2254., 1677.,  976.,  570.,  199.,   82.,
          32.,    6.,    0.,    0.,    0.]),
 array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14]),
 <a list of 14 Patch objects>)
In [28]:
# Observation: the above two distributions (n=100,1000) are very close 
# Suggest that as n-> infinity, the disributions binom(n,lam/n)
# converge to some limiting distribution.  
# This limit is the Poisson distribution.
In [29]:
# Try out different values of lam
In [31]:
lam=1.0 # lambda, intensity of Poisson distribution
n=1000
plt.hist([binom(n,lam/n) for i in range(10000)], bins=range(0,15), align="left", rwidth=0.5)
Out[31]:
(array([3.692e+03, 3.727e+03, 1.790e+03, 6.170e+02, 1.520e+02, 1.900e+01,
        3.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00,
        0.000e+00, 0.000e+00]),
 array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14]),
 <a list of 14 Patch objects>)
In [33]:
lam=8.0 # lambda, intensity of Poisson distribution
n=1000
plt.hist([binom(n,lam/n) for i in range(10000)], bins=range(0,30), align="left", rwidth=0.5)
Out[33]:
(array([3.000e+00, 4.000e+01, 9.400e+01, 2.530e+02, 5.790e+02, 9.330e+02,
        1.160e+03, 1.455e+03, 1.402e+03, 1.267e+03, 1.011e+03, 7.180e+02,
        4.740e+02, 2.900e+02, 1.500e+02, 8.700e+01, 4.800e+01, 2.100e+01,
        7.000e+00, 5.000e+00, 1.000e+00, 1.000e+00, 1.000e+00, 0.000e+00,
        0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00]),
 array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16,
        17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29]),
 <a list of 29 Patch objects>)
In [ ]:
# Poisson(lam) tends to a normal distribution (with proper scaling)
# as lam->infinity
In [34]:
# Poisson can also be characterized by its Probability Mass Function (pmf)
# (analogue of Probability Density Functions, but for discrete distribution)
# PMF is a function f, such that f(k) gives prob that distribution takes 
# the value k. 
# For instance, {0,1}-discrete r.v. (i.e. randint(0,1))
# PMF is f(0)=1/2, f(1)=1/2, f(k) =0 for any other k.
In [35]:
#PMF for Poisson(lam):
# f(k) = lam^k * exp(-lam)/k! 
# e.g. for Poisson(1), f(k)= exp(-1)/k! 
# (By convention 0!=1. Note f is not even defined for k not 
# non-negative integer)
f(k)=lam^k * exp(-lam)/factorial(k)
In [37]:
lam = 1
f(k)=lam^k * exp(-lam)/factorial(k)
list_plot([f(k) for k in range(0,10)])
Out[37]:
In [39]:
N(f(0)) # prob that Poisson(1) takes value 0
Out[39]:
0.367879441171442
In [40]:
N(f(1))
Out[40]:
0.367879441171442
In [41]:
N(f(6))
Out[41]:
0.000510943668293670
In [42]:
N(f(7))
Out[42]:
0.0000729919526133814
In [45]:
N(f(30)) # non-zero, but very small
Out[45]:
1.38690094211205e-33
In [46]:
lam = 3
f(k)=lam^k * exp(-lam)/factorial(k)
list_plot([f(k) for k in range(0,10)])
Out[46]:
In [ ]:
# degrees of vertices in random graphs 
# (for n vertices, n edges, with n large)
# will be close to Poisson distributed