Python numpy.convolve 积分限制答案

【问题标题】：Python numpy.convolve integral limitsPython numpy.convolve 积分限制
【发布时间】：2018-02-10 00:11:12
【问题描述】：

我有一个卷积积分，它与这篇文章中描述的形式相似：(link)

我需要从负无穷积分到 t，而不是从 0 积分到 t。但是，我无法使用numpy.convolve 做到这一点，因为它总是返回从负无穷到正无穷的结果。使用scipy.integrate.quad 会非常慢，因为我必须遍历每个t，并且它只适用于具有解析表达式的被积函数。

有没有办法指定numpy.convolve 的下限和上限？非常感谢。

这是代码（我很抱歉无法在这里输入 LaTeX 方程）：

import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate

def gaussian(tau, mu, sigma):

    pdf = 1/(sigma * np.sqrt(2 * np.pi)) * np.exp( - (tau - mu)**2 / (2 * sigma**2))
    return pdf

def gaussian_deriv(tau, mu, sigma):
    """
    derivative of gaussian function
    """
    pdf = -(tau - mu)/sigma**2 * gaussian(tau, mu, sigma)
    return pdf

def integral_kernel(tau):
    return np.cbrt(1/tau)

def integrand(tau, t, mu, sigma):

    return gaussian_deriv(tau, mu, sigma) * integral_kernel(t - tau + 1E-28)

tau = np.linspace(-7, 7, 1000)
dtau = tau[1] - tau[0]
lower_lim = tau[0]
g_deriv = gaussian_deriv(tau, mu=0, sigma=1)

result = np.zeros_like(tau)
for idx, t in np.ndenumerate(tau):
    result[idx], err = integrate.quad(integrand, lower_lim, t, args=(t, mu, sigma), points=[t])

result_convolve = np.convolve(g_deriv, integral_kernel(tau), mode='same') * dtau

fig, ax = plt.subplots(2,1, figsize=(10, 6))
ax[0].plot(tau, result, 'r-', label='scipy quad')
ax[1].plot(tau, result_convolve, '.', label='numpy convolve')
ax[0].legend()
ax[1].legend()
plt.show()

【问题讨论】：

标签： python numpy scipy convolution

【解决方案1】：

如果我理解正确，你想要积分

int{-inf...T} f(tau) g(t-tau) dtau

写H(t) = 1/2 (sign(t) + 1) 这和

int{-inf...inf} f(tau) g(t-tau) [1 - H(tau-T)] dtau

或

int{-inf...inf} f^(tau) g(t-tau) dtau = f^*g(t)

其中 f^(tau) = f(tau) * [1 - H(tau-T)]。

所以你需要做的就是将 T 右边的 f 归零以获得 f^，然后计算常规卷积 f^*g。

更新

如果 t 和 T 是同一个数字，那就更简单了：

int{-inf...t} f(tau)g(t-tau) dtau

在这种情况下，您可以将积分变量移动 t，因此您得到：

int{-inf...0} f(tau+t}g(-tau) dtau

接下来将 tau 替换为 -tau

int{0...inf} g(tau) f(t-tau) dtau

现在你或多或少地像上面那样进行，只有 T 现在是固定的并且等于零，而且你不是在 T 的右边而是在 T 的左边归零。

【讨论】：

当 tau > T 时如何将 f(tau) 归零？我尝试了步进函数numpy.piecewise，如果 tau > 0，f =0，但它不起作用。相反，我得到的结果非常接近实际结果，但略有偏差。
你可以这样做：def cut(tau, f, T):out=np.zeros(np.shape(tau))tau = np.asanyarray(tau)mask = tau <= Tout[mask] = f(tau[mask])return out
我有点困惑。函数numpy.convolve 将结果作为 T 的函数返回（与 T 大小相同的数组）。如果我按照您上面的建议继续，我将不得不遍历 T 的每个元素并评估积分，不是吗？如果可能的话，我想坚持使用numpy.convolve。
@misterdone 我写的假设T 固定标量，tau 积分变量，t 卷积偏移量（一个向量）一旦你计算了 f^，相当于调用@987654333 @ 一次在网格上 t 你打电话给convolve 在结果和g 上也一次，就是这样
@misterdone 如果你的意思不是 t=T，结果不同但相似，我也更新了帖子以涵盖这种情况。

【解决方案2】：

这又是一篇相当老的帖子，我仍然在这里展示我的分析。

请参阅this post 我已经回答了 OP 引用的问题。有必要通过这篇文章来了解这里使用的符号和符号。

np.convolve(A, B) 给你的是

而你正在寻找的是

这里第二项将为零，因为负系数的B 被外推为零。所以你看到你不需要做任何事情

【讨论】：