Python 3.4 SciPy integrate.quad

Question

оставите машину Я пытаюсь вычислить интеграл гауссовой в питона так:

from math import exp 
from scipy import stats, integrate 
import scipy.interpolate as interpolate 
from numpy import cumsum, random, histogram, linspace, zeros, inf, pi,sqrt 
import matplotlib.pyplot as plt 

A = 1 
mu = 0 
sigma = 1 
p = lambda x: A * exp(-(((x-mu)**2))/(2*(sigma**2))) 
F = lambda x: integrate.quad(p, -inf, x)[0] 
Ns = 1000; 
x = linspace(-50,50,Ns); 
y = zeros(Ns) 
yy = zeros(Ns) 
for i in range(Ns): 
    y[i] = F(x[i]) 
    yy[i]= p(x[i]) 

plt.plot(x,y) 
plt.plot(x,yy) 
plt.show() 

, но если посмотреть на участке, есть падение до нуля в диапазоне между 21,0 до 22, а после 38+. Кто-нибудь знает, почему это так? Возможно, ошибки округления?

спасибо !!

Answer 1

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

Гауссово быстро переходит в ноль при отклонении от среднего значения, поэтому в основном на интервале между (-50, 50) большинство значений функции равны нулю. Если метод интеграции не может отбирать точки из вашей небольшой области, где функция отлична от нуля, она увидит, что вся функция полностью плоская и, таким образом, дает вам интеграл 0.

Итак, что вы можете сделать?

Вместо того, чтобы выбирать фиксированный интервал (-50,50), выберите интервал, основанный на меньших значениях сигма, чтобы избежать интегрирования по слишком большому интервалу нулей.

Если вы идете только на 5, 10 или 20 стандартных отклонений влево и вправо, вы не увидите эту проблему, и у вас все еще есть очень точный результат интеграции.

Это результат, если вы интегрируете 10 стандартных отклонений влево и вправо.