当一个函数在大部分域上为零时，数值计算卷积

Question

我正在尝试计算与 LogNormal 卷积的高斯概率密度函数。与 LogNormal 的宽度相比，Gaussian 的宽度非常小，因此对于大部分积分域 (0, np.inf).

，Gaussian 约为 0

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()

如您所见，对于高斯 (x) ~ 0 的某些部分，积分从正确计算的值跳到 ~0。阅读此文：Discontinuity in results when using scipy.integrate.quad 我的想法是我通过以下方式帮助数值积分通过更改变量 t 从 (0, 1) 而不是 (0, inf) 映射集成域。

------ 编辑：最初的问题是我的数学错误引起的。 -----

因此卷积积分必须改为： PDF(x) = INT_0^inf F(x,t)dt = INT_0^1 F(x, s/(1-s)) /(1-s)^2 ds 其中 F(x,t) 是卷积被积函数，c.f @Juan Carlos Ramirez 回答。

函数_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

此修复改善了这种情况（见下图：顶部图像未转换边界，底部图像有转换），但是，集成仍然不令人满意。我该如何解决这个问题？

Answer 1

您的陈述： PDF(x) = INT_0^inf F(x,t)dt = INT_0^1 F(x, s) (1+s)^2 ds 不太正确，因为您需要将 t 替换为 F(x,t) 中 s 的函数。请注意

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

quad 函数非常强大，但由于它基本上尝试了一堆不同的数值积分技术，因此可能缺乏透明度。如果您使用梯形规则通过将 [-a,a] 划分为 N 个点的均匀网格来计算从 -a 到 lognormal(x-t)gauss(t) 的 a 的积分会怎么样？由于高斯分布以超指数方式衰减，您可以通过将区间限制为 [-a,a] 来控制误差，同样地，您可以通过查看误差公式（按二次比例缩放，你可以绑定高斯*对数正态的二阶导数）。 https://en.wikipedia.org/wiki/Trapezoidal_rule#Uniform_grid 这将是手波式的，但梯形法则在快速衰减函数上可以有非常好的行为，例如参见 [=15=]。

当一个函数在大部分域上为零时，数值计算卷积

numerically computing convolution when one function is zero over large part of domain

python

numerical-integration