In [1]:
# Empirical deomonstration of Poisson distribution 

import matplotlib.pyplot as plt
intens = 1/3 # intensity(lambda) = expected number of hits in one unit of time 
num_trials = 10000 
time_steps = 50 # divide our unit of time into this many
# pieces.  In each piece prob of getting a hit 
# is intens/time_steps
# Poisson will be limit as time_steps -> infinity

results = [] # store number of hits for each trial
for trial in range(num_trials):
    num_hits = 0  # keeps track of the number of hits for this trial
    for t in range(time_steps):
        # in this step, prob of getting hit is intens/time_steps
        hit = 0 # will be 0 if not hit, 1 if hit
        rand = random() # uniform (0,1) random variable
        if rand < intens/time_steps: # happens with prob intens/time_steps
            hit = 1
        else:
            hit = 0
        num_hits += hit  # increase number of hits by the numb of hits from this step
    results.append(num_hits)
plt.hist(results, bins=range(0,10), align="left", rwidth=0.5)
# can set bins to be a list
Out[1]:
(array([7.140e+03, 2.414e+03, 4.010e+02, 4.300e+01, 2.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]),
 <a list of 9 Patch objects>)
In [2]:
# Increase intensity 
import matplotlib.pyplot as plt
intens = 1 # intensity(lambda) = expected number of hits in one unit of time 
num_trials = 10000 
time_steps = 50 # divide our unit of time into this many
# pieces.  In each piece prob of getting a hit 
# is intens/time_steps
# Poisson will be limit as time_steps -> infinity

results = [] # store number of hits for each trial
for trial in range(num_trials):
    num_hits = 0  # keeps track of the number of hits for this trial
    for t in range(time_steps):
        # in this step, prob of getting hit is intens/time_steps
        hit = 0 # will be 0 if not hit, 1 if hit
        rand = random() # uniform (0,1) random variable
        if rand < intens/time_steps: # happens with prob intens/time_steps
            hit = 1
        else:
            hit = 0
        num_hits += hit  # increase number of hits by the numb of hits from this step
    results.append(num_hits)
plt.hist(results, bins=range(0,10), align="left", rwidth=0.5)
# can set bins to be a list
Out[2]:
(array([3.631e+03, 3.725e+03, 1.879e+03, 6.090e+02, 1.270e+02, 2.300e+01,
        5.000e+00, 1.000e+00, 0.000e+00]),
 array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9]),
 <a list of 9 Patch objects>)
In [4]:
# Increase intensity 
import matplotlib.pyplot as plt
intens = 4 # intensity(lambda) = expected number of hits in one unit of time 
num_trials = 10000 
time_steps = 50 # divide our unit of time into this many
# pieces.  In each piece prob of getting a hit 
# is intens/time_steps
# Poisson will be limit as time_steps -> infinity

results = [] # store number of hits for each trial
for trial in range(num_trials):
    num_hits = 0  # keeps track of the number of hits for this trial
    for t in range(time_steps):
        # in this step, prob of getting hit is intens/time_steps
        hit = 0 # will be 0 if not hit, 1 if hit
        rand = random() # uniform (0,1) random variable
        if rand < intens/time_steps: # happens with prob intens/time_steps
            hit = 1
        else:
            hit = 0
        num_hits += hit  # increase number of hits by the numb of hits from this step
    results.append(num_hits)
plt.hist(results, bins=range(0,20), align="left", rwidth=0.5)
# can set bins to be a list
Out[4]:
(array([ 163.,  684., 1395., 2015., 2037., 1577., 1078.,  578.,  298.,
         129.,   37.,    9.,    0.,    0.,    0.,    0.,    0.,    0.,
           0.]),
 array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16,
        17, 18, 19]),
 <a list of 19 Patch objects>)
In [56]:
# plot the pmf of Poisson distribution

lam = 4 # intensity
pmf_list = [lam^k * exp(-lam)/factorial(k) for k in range(0,10)]
list_plot(pmf_list)
Out[56]:
In [63]:
# plot pmf and experimental dist on same list_plot
import matplotlib.pyplot as plt
import numpy as np
intens = 1 # expected number of hits in one unit of time 
num_trials = 10000 
time_steps = 50 # divide our unit of time into this many
#pieces.  Will get large. In each piece prob of getting a hit 
# is intens/time_steps
# Poisson is limit as time_steps -> infinity

hits_type = np.zeros(num_trials) # ith element will be the number of trials that 
# got i hits (divided by num_trials to give a frequency)
for trial in range(num_trials):
    num_hits = 0
    for t in range(time_steps):
        # in this step, prob of getting hit is intens/time_steps
        hit = 0 # will be 0 if not hit, if 1 hit
        rand = random()
        if rand < intens/time_steps:
            hit = 1
        else:
            hit = 0
        num_hits += hit
    hits_type[num_hits] +=1/num_trials

lam = 1 # intensity
# plot the pmf of Poisson distribution
pmf_list = [lam^k * exp(-lam)/factorial(k) for k in range(0,10)]
list_plot(pmf_list[0:30], color="red",marker="<")\
+ list_plot(hits_type[0:30])
# theoretical pmf plotted in red
# experimental results in blue
# they will be very close
Out[63]: