численно вычисляя свертку, когда одна функция равна нулю на большей части области - PullRequest
Новая
       37

численно вычисляя свертку, когда одна функция равна нулю на большей части области

0 голосов
/ 06 июня 2019

Я пытаюсь вычислить функцию плотности вероятности гауссиана, свернутого с логнормой.Ширина гауссиана очень мала по сравнению с шириной LogNormal, таким образом, гауссиан равен ~ 0 для большей части области интегрирования (0, np.inf). 

from scipy.integrate import quad
import numpy as np
def _gaussian(x, mu, sigma):
    return 1/np.sqrt(2*np.pi*sigma**2)*np.exp(-1*(x - mu)**2/(2*sigma**2))
def _lognormal(x, mu_log, sigma_log):
    return 1/(x*sigma_log)*1/np.sqrt(2*np.pi)*np.exp(-1*(np.log(x) - mu_log)**2/(2*sigma_log**2)  )
def _gauss_log(t, x0, mu, sigma, mu_log, sigma_log):
    return  _lognormal(t, mu_log, sigma_log) * _gaussian(x0-t, mu, sigma)

def gauss_log(x,  mu, sigma, mu_log, sigma_log):
        return  [quad(_gauss_log, 0, np.inf, args=(x0, mu, sigma, mu_log, sigma_log))[0] for x0 in x]

x = [ 0.06898463,  0.12137053,  0.21353749,  0.37569469,  0.66099163,
        1.16293883,  2.04605725,  3.59980263,  6.33343911, 11.14295839,
       19.60475495]
y = [0.00000000e+00, 8.31638354e-02, 1.65440428e-01, 4.02998983e-01,
        5.42100908e-01, 4.16612321e-01, 1.72662493e-01, 5.60788435e-02,
        1.43433519e-02, 5.43498669e-03, 2.57428324e-04]
x00 = np.linspace(0.01, 20, 100)

plt.loglog(x, y, 'o')
plt.ylim([0.0001, 10])
plt.plot(x00, gauss_log(x00, 0, 0.05, 0.1, 0.9))
plt.show()

Proper Convolution Integral

Как вы можете видеть, интегральные скачки от правильно вычисленного значения до ~ 0 для некоторых частей, где гауссовский (x) ~ 0. Чтение этого: Прерывистость в результатах при использовании scipy.integrate.quad Я думал, что помогаю числовому интегрированию, отображая область интегрирования из (0, 1) вместо (0, inf), изменяя переменную t.

------ РЕДАКТИРОВАТЬ: Первоначальная проблема была вызвана ошибкой в ​​моей математике. -----

Следовательно, интеграл свертки должен бытьизменено на: PDF(x) = INT_0^inf F(x,t)dt = INT_0^1 F(x, s/(1-s)) /(1-s)^2 ds , где F(x,t) - интегрант свертки, ср. @Juan Carlos Ramirez answer.

Функция _gauss_log становится: 

def _gauss_log(t, x0, mu, sigma, mu_log, sigma_log):
    return  _lognormal(s/(s-1), mu_log, sigma_log) * _gaussian(x0 - s/(s-1), mu, sigma)/(1-s)**2

Это исправление улучшает ситуацию (см. рисунок ниже: верхнее изображение без преобразованных границ, нижнее изображение с преобразованием), однако интеграция все еще не является удовлетворительной.Как я могу решить эту проблему?

transformed integration boundaries vs. open integral

1 Ответ

1 голос
/ 06 июня 2019

Ваше утверждение: PDF(x) = INT_0^inf F(x,t)dt = INT_0^1 F(x, s) (1+s)^2 ds не совсем верно, так как вам нужно заменить t как функцию s внутри F (x, t).Обратите внимание, что

t = 1/(1-s) - 1 = s/(1-s) поэтому 

dt = 1/(1-s)^2 ds

и поэтому ваш преобразованный интеграл должен быть

INT_0^1 F(x, s/(1-s))/(1-s)^2 ds

Функция квадрата довольно мощная,но поскольку он в основном использует кучу различных методов численного интегрирования, ему может не хватать прозрачности.Что если вы используете трапециевидное правило для вычисления интеграла от -a до a логнормального (xt) гаусса (t) путем деления [-a, a] на равномерную сетку из N точек?Поскольку распределение Гаусса затухает суперэкспоненциально, вы имеете большой контроль над ошибкой, вызванной ограничением интервала [-a, a], и аналогичным образом вы можете контролировать, насколько большим вы хотите, чтобы N было, просматривая формулу ошибки (масштабируется квадратично,и вы можете связать вторую производную от гауссовской * логнормальной).https://en.wikipedia.org/wiki/Trapezoidal_rule#Uniform_grid Это будет волнообразно, но правило трапеции может иметь очень хорошее поведение на быстро затухающих функциях, см., Например, https://www.ams.org/journals/mcom/1964-18-087/S0025-5718-1964-0185804-1/S0025-5718-1964-0185804-1.pdf.

...