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

Question

Когда сигма2 велика, интеграция дает правильный ответ. но при уменьшении сигма2 до 1е-8 или ниже, он возвращает 0. Есть идеи, как решить проблему?

Я проверил Matlab, который дает правильный ответ, независимо от того, насколько мала сигма2. Спасибо!

quad (лямбда x: 1/sqrt(2*pi*sigma2)*exp(-x**2/(2*sigma2)), -10, 10)

Answer 1

Короткий ответ

Используйте необязательный параметр points чтобы помочь интегратору найти пик Гаусса:

quad(lambda x: 1/sqrt(2*pi*sigma2)*exp(-x**2/(2*sigma2)), -10, 10, points=[-10*sqrt(sigma2), 10*sqrt(sigma2)])

Вывод является точным, даже если sigma2 = 1e-100,

Значение points указывает на изменение поведения функции: от 0 (для всех практических целей) он поворачивается и образует очень резкий пик в пределах нескольких стандартных отклонений от среднего значения.

Более простой способ - просто заменить интервал интеграции на -10*sqrt(sigma2), 10*sqrt(sigma2), но это может быть не то, что вы хотите, если интеграл включает в себя другие термины, кроме гауссова.

объяснение

Проблема заключается в том, что когда пик Гаусса очень тонкий (например, ширина 0,001 в интервале длины 20), и процедура интеграции, скорее всего, никогда его не увидит. Интеграторы, как правило, хорошо адаптируются к функциям функций, которые им предоставляются, но они не могут адаптироваться к функциям, которые они не видят.

Хорошая производительность MATLAB's quad на этом примере в основном удача: он использует адаптивный метод Симпсона, который включает выборку функции точно в середине интервала интегрирования. Если вы дадите ему несимметричный интервал, такой как (-10,11), он будет работать намного хуже, чем у Сципи. quad, Вот сравнение (используя мою копию MATLAB, которая R2013a):

• quad(...,-10,10) всегда возвращает 1
• quad(...,-10,11) возвращает 1 для sigma2=0.1 и не подходит для sigma2=0.01 и меньше.
• scipy.integrate.quad(..., -10, 10) возвращает 1 до 1e-4
• scipy.integrate.quad(..., -10, 11) возвращает 1 до 1e-3
• от Matlab integral(...) возвращает 1 до 1e-5. Будет ли интервал (-10,10) или (-10,11) не имеет значения для него.

Стоит отметить, что MATLAB quad устарела, и MathWorks рекомендует использовать integral что, как вы можете видеть, делает гораздо лучше, когда удача не участвует.

Разница в производительности между MATLAB integral и SciPy's quad не драматично; для сильно локализованных гауссиан (относительно интервала интегрирования) оба будут нуждаться в помощи в поиске пика. Помощь приходит в виде points параметр scipy.integrate.quad или же WayPoints параметр integral,