Расчет проблем с большими биномиальными коэффициентами

Question

Вступление.

Я пытаюсь построить гипергеометрическое распределение с помощью IPython. Функция вероятности распределения содержит три биномиальных коэффициента.

Поскольку значения, которые я буду вводить в коэффициенты, очень велики, например. 1e28, я решил создать свою собственную функцию для вычисления биномиальных коэффициентов, где я использую второй порядок приближения Стирлинга.

Поскольку биномиальные коэффициенты слишком велики, чтобы помещать их в переменную и умножать напрямую, я решил вместо этого рассчитать их логарифмы и добавить их друг к другу. Чтобы получить окончательную вероятность, я просто помещаю результат в exp функция. Поскольку вероятность относительно "нормального" размера (ее максимум равен 2.7e24), проблем больше не должно быть... кроме как.

Проблема.

"Результат", на который я ссылался, - это логарифм вероятности, которая оказалась намного больше, чем должна быть, со значениями в диапазоне от -6.2e24 до 1.3e14. Для сравнения, лог математического максимума составляет примерно 56.

Другая проблема заключается в том, что кривая на графике очень сильно зазубрена при увеличении. Все выглядит хорошо при уменьшении масштаба. Кривая гладкая, ее максимум находится в среднем от распределения:

Участок по всему домену.

Но когда я увеличиваю среднее значение, где будет находиться большая часть пика функции вероятности, поскольку стандартное отклонение очень мало, я получаю это:

Увеличено в среднем.

Красная линия обозначает среднее, а черные - среднее +/- стандартное отклонение. Хотя это выглядит красиво, это не то, что мне нужно, это плавная кривая.

Вопрос.

Может кто-нибудь объяснить, почему (1) значения такие большие и (2) почему кривая зазубрена и как я могу их исправить?

Код.

import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline

#Returns the log of "n choose k" calculated with 2nd order Stirling's approximation
def l_cb(n, k):
    if (k > n) or (n < 0) or (n < 0):
        print "Invalid values for the binomial coefficient:", "n =", n, ", k =", k, "."
        return 0.0
    if (k == n) or (k == 0) or (n == 0):
        return 0.0
    A = (n + 0.5) * np.log(n) - (k + 0.5) * np.log(k) - (n - k + 0.5) * np.log(n - k)
    B = np.log(1 + 1 / (12.0 * n)) - np.log(1 + 1 / (12.0 * k)) - np.log(1 + 1 / (12.0 * (n - k)))
    return - 1/2 * np.log(2 * np.pi) + A + B

K = 2.24e28
k = 2.24e27
N = 2.7e25

#Mathematical maximum of P. np.log(MAX) is about 56. l_P is way too big.
MAX = (k + 1) * (N + 1) / (K + 2)
#Mathematical average of n
AVG = N * k / K
#Mathematical standard deviation of P.
SD = np.sqrt(N * k * (K - N) * (K - k) / K ** 2 / (K - 1))

n = np.linspace(AVG - 50e12, AVG + 50e12, 1001)
l_P = np.zeros(len(n))

#Calculating log(P).
for i in xrange(len(l_P)):
    l_P[i] = l_cb(N, n[i]) + l_cb(K - N, k - n[i]) - l_cb(K, k)

#Marking AVG, AVG - SD, AVG + SD
y = np.linspace(-4e14, 5e14, len(n))
x_AVG = np.ones(len(n))
x_SD_L = np.ones(len(n)) - SD / AVG
x_SD_R = np.ones(len(n)) + SD / AVG

plt.plot(n / AVG, l_P, x_AVG, y, 'r', x_SD_L, y, 'k', x_SD_R, y, 'k')
plt.xlabel('n / AVG')
plt.ylabel('log(P)')