[Python] Comparing numpy versus scipy when sampling random values from a poisson distribution

I am generating random samples from both numpy and scipy for a Poisson distribution, but it appears that scipy never produces values of x=0. Scipy seems to be off by x=1, i.e., it is shifted to the right by 1. Can someone check this? Thanks.

# Check Poisson samples
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import gamma, poisson

N = 100000
plot_low=0
plot_high=40
x = np.arange(plot_low, plot_high)
mission_time = 40
lambda_val = 0.2
poisson_mu = lambda_val*mission_time

Poisson_np = []
Poisson_scipy = []
for i in range(N):
x_simulated_scipy = poisson.rvs(poisson_mu, 1) # NOTE this function is off by x=1
x_simulated_np = np.random.poisson(lam=poisson_mu)
# Now randomly sample from the Poisson model
Poisson_np.append(x_simulated_np)
Poisson_scipy.append(x_simulated_scipy)

Poisson_np_pmf = []
Poisson_scipy_pmf = []
for i in range(plot_high):
# Count how many times we saw X=i then divided by N to get the probability X=i
np_pmf = Poisson_np.count(i)/N
Poisson_np_pmf.append(np_pmf)
scipy_pmf = Poisson_scipy.count(i)/N
Poisson_scipy_pmf.append(scipy_pmf)
if i == 0:
print("NUMPY zero count =", Poisson_np.count(0))
print("SCIPY zero count =", Poisson_scipy.count(0))

plt.plot(x, Poisson_np_pmf, '.', ms=5, mec='b', alpha=0.7)
plt.vlines(x, 0, Poisson_np_pmf, colors='b', lw=3, label='Poisson via NUMPY', alpha=0.5)
plt.plot(x, Poisson_scipy_pmf, '.', ms=5, mec='r', alpha=0.5)
plt.vlines(x, 0, Poisson_scipy_pmf, colors='r', lw=3, label='Poisson via SCIPY', alpha=0.5)
plt.xlabel('x')
plt.ylabel('f(x)')
plt.title('Poisson Distribution (40 year duration)')
plt.grid(True)
plt.legend()
plt.show()


Output Plot from Test Script

Continue reading...