This is a homework of mine to plot the gaussian and poisson distribution for r = 20, mu = 20 and sigma = sqrt(20). But the problem is that the value when computing of the computer is capped somewhere, and the large number like 20^20 cannot be computed, so the graph looks horribly wrong after r >= 8. Here is my code:
import numpy as np
import matplotlib.pyplot as plt
import sympy as sp
from scipy.special import factorial
#c)
sigma = np.sqrt(20)
mu1 = 20
b = np.sqrt(2*np.pi)
def fGauss(r):
return 1/(b * sigma) * np.exp(-(r - mu1)**2 / (2*sigma**2))
def fPoisson2(r):
return mu1**r * np.exp(-mu1)/factorial(r)
x1 = np.arange(1, 21, 1)
Gauss = [fGauss(r) for r in x1]
Poisson2 = [fPoisson2(r) for r in x1]
plt.bar(x1, Poisson2, color = 'b', alpha = 1, label = 'Poisson-Verteilung')
plt.bar(x1, Gauss, color = 'r', alpha = 0.2, label = 'Gauss-Verteilung')
plt.xlim(-0.1, 20.1)
plt.xticks(x1)
plt.xlabel('r')
plt.ylabel('Wahrscheinlichkeit')
plt.title(r'Poisson- und Gauss-Verteilung für $(\mu = 20; r = 1,2,...,20)$')
plt.show()
I have asked a friend of mine for advice and he said to use np.random.poisson, but I still don't know how to work around it, so I post this here, hoping for help. In attachment is my file. If anyone has question regarding the task, I will provide as much as I could. Many thanks
Your intermediate working contains huge (but partially-cancelling) numbers like 20^20 and 20!
You can avoid this by constructing the successive probabilities incrementally, multiplying by mu/i for i = 1, ... r.
Change your Poisson probability function to
def fPoisson2(r):
answer = np.exp( -mu1 )
for i in range( 1, r + 1 ): answer *= mu1 / i
return answer
and you will be OK.
Note that the normal distribution will be a close approximation to Poisson for mu this large.