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]:

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]:

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]:

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]: