scipy integrate.quad возвращает неправильное значение

Question

Я использую scipy integrate.quad для расчета cdf нормального распределения:

def nor(delta, mu, x):
    return 1 / (math.sqrt(2 * math.pi) * delta) * np.exp(-np.square(x - mu) / (2 * np.square(delta)))


delta = 0.1
mu = 0
t = np.arange(4.0, 10.0, 1)
nor_int = lambda t: integrate.quad(lambda x: nor(delta, mu, x), -np.inf, t)
nor_int_vec = np.vectorize(nor_int)

s = nor_int_vec(t)
for i in zip(s[0],s[1]): 
    print i

пока он печатает следующим образом:

(1.0000000000000002, 1.2506543424265854e-08)
(1.9563704110140217e-11, 3.5403445591955275e-11)
(1.0000000000001916, 1.2616577562700088e-08)
(1.0842532749783998e-34, 1.9621183122960244e-34)
(4.234531567162006e-09, 7.753407284370446e-09)
(1.0000000000001334, 1.757986959115912e-10)

для некоторых x возвращается значение, близкое к нулю, оно должно быть возвращено 1. Может кто-нибудь сказать мне, что не так?

Answer 1

По той же причине, что и в квадрате, возвращаются оба нуля при интегрировании простого гауссова pdf с очень малой дисперсией? но, поскольку я не могу пометить его как дубликат, вот так:

Вы интегрируете функцию с узкой локализацией (в дельте масштаба) на очень большом (фактически бесконечном) интервале. Процедура интегрирования может просто пропустить часть интервала, в которой функция существенно отличается от 0, вместо этого оценивая ее как 0. Некоторое руководство требуется. Параметр points может быть использован для этого (см. связанный вопрос), но так как quad в течение бесконечного интервала это не поддерживается, интервал должен быть разделен вручную, например так:

for t in range(4, 10):
    int1 = integrate.quad(lambda x: nor(delta, mu, x), -np.inf, mu - 10*delta)[0] 
    int2 = integrate.quad(lambda x: nor(delta, mu, x), mu - 10*delta, t)[0] 
    print(int1 + int2)

Это печатает 1 или почти 1 каждый раз. я выбрал mu-10*delta в качестве точки, на которую можно разделить, вычисление большей части функции находится справа от нее, независимо от того, что такое mu и delta.

Заметки:

1. использование np.sqrt так далее; обычно нет причин ставить math функции в коде NumPy. Версии NumPy доступны и векторизованы.
2. применение np.vectorize в quad ничего не делает, кроме того, что делает код длиннее и немного сложнее для чтения. Используйте обычный цикл Python или понимание списка. Смотрите векторизацию NumPy с интеграцией