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()
poisson.rvs(poisson_mu, 1)
You're shifting the distribution over by 1. Per the documentation scipy.stats.poisson.rvs(mu, loc=0, size=1, random_state=None), so your 1 is being interpreted as the loc argument based on the position of the argument. Were you trying to make sure the size was 1? If so, be sure to pass it as a keyword argument rather than positional (even though size=1 is the default).