时间戳

首页 > 解决方案 > 当一个函数在大部分域上为零时数值计算卷积

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

问题描述

我正在尝试计算与 LogNormal 卷积的高斯概率密度函数。与 LogNormal 的宽度相比，高斯的宽度非常小，因此对于大部分积分域（0，np.inf），高斯为 ~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()

适当的卷积积分

正如您所看到的，对于 Gaussian(x) ~ 0 的某些部分，积分从正确计算的值跳到 ~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 的回答。

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

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

转换积分边界与开放积分

标签： pythonnumerical-integration

解决方案

您的陈述： 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 个点的均匀网格来计算对数正态 (xt)gauss(t) 从 -a 到 a 的积分会怎样？由于高斯分布以超指数方式衰减您可以很好地控制将区间限制为 [-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 。

推荐阅读