Интегрирование гауссиана в течение очень длинного интервала - PullRequest
Новая
       20

Интегрирование гауссиана в течение очень длинного интервала

0 голосов
/ 13 мая 2019

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

from scipy.integrate import quad
from math import pi, exp, sqrt

def func(x, mean, sigma):
    return 1/(sqrt(2*pi)*sigma) * exp(-1/2*((x-mean)/sigma)**2) 

print(quad(func, 0, 1e+31, args=(1e+29, 1e+28))[0]) # case 1
print(quad(func, 0, 1e+32, args=(1e+29, 1e+28))[0]) # case 2
print(quad(func, 0, 1e+33, args=(1e+29, 1e+28))[0]) # case 3
print(quad(func, 1e+25, 1e+33, args=(1e+29, 1e+28))[0]) # case 4

, то печатаются следующие.

1.0
1.0000000000000004
0.0
0.0

Чтобы получить разумный результат, мне пришлось несколько раз попытаться изменить нижнюю / верхнюю границы интеграла и эмпирически определить его как [0, 1e + 32].Мне это кажется рискованным, так как, когда меняются среднее значение и сигма гауссовой функции, мне всегда приходится пытаться использовать разные границы.

Есть ли четкий способ интегрировать функцию от 0 до 1e + 50, не беспокоясь о границах?Если нет, как вы ожидаете от начала, какие границы дадут ненулевое значение?

Ответы [ 2 ]

1 голос
/ 13 мая 2019

Короче, вы не можете.

На этом длинном интервале область, где гауссиан отличен от нуля, является крошечной, и адаптивная процедура, которая работает под капотом integrate.quad, не видит его. Как и любая другая адаптивная рутина, если не случайно.

0 голосов
/ 14 мая 2019

Уведомление,

enter image description here

, а CDF нормальной случайной величины известен как ϕ(x), поскольку он не может быть выражен элементарная функция .Так что бери ϕ((b-m)/s) - ϕ((a-m)/s).Также обратите внимание, что ϕ(x) = 1/2(1 + erf(x/sqrt(2))), поэтому вам не нужно звонить .quad для фактического выполнения интеграции, и вам может повезти с erf из scipy.

from scipy.special import erf

def prob(mu, sigma, a, b):
    phi = lambda x: 1/2*(1 + erf((x - mu)/(sigma*np.sqrt(2))))
    return phi(b) - phi(a)

Это может дать более точные результаты (чем выше) 

>>> print(prob(0, 1e+31, 0, 1e+50))
0.5
>>> print(prob(0, 1e+32, 1e+28, 1e+29))
0.000359047985937333
>>> print(prob(0, 1e+33, 1e+28, 1e+29))
3.5904805169684195e-05
>>> print(prob(1e+25, 1e+33, 1e+28, 1e+29))
3.590480516979522e-05

и избежать возникшей интенсивной ошибки floating point.Однако области, которые вы интегрируете, настолько малы по площади, что вы все равно можете увидеть 0.

...