n = 450000
p = 0.5
k = np.array([17.0, 225000, 226000])
b = scipy.stats.binom(n, p)
print(b.pmf(k))
n = scipy.stats.norm(n * p, np.sqrt(n * p * (1 - p)))
print(n.pdf(k))
print(b.pmf(k) - n.pdf(k))
Вычисление биномиальной вероятности для огромных чисел