scipy 卷积取决于 x

Question

我正在尝试对一个对数正态 PDF 和一个高斯 PDF 进行卷积。因此，我按以下方式定义了函数：

def PDF_log(x,sig,mu): # log normal PDF
   mu = np.log(mu)
   return( (1/x)*(1/(sig*np.sqrt(2*np.pi))) * np.exp(-(np.log(x)-mu)**2/(2*sig**2)) )   

def gauss(x,sig,mu): # a noraml PDF
   return( 1/(np.sqrt(2*np.pi*sig**2)) * np.exp(-(x-mu)**2/(2.*sig**2)) )

def gauss_log(x, sig, mu, sig0, mu0):
    a = signal.convolve(PDF_log(x,sig,mu),gauss(x,sig0,mu0),mode='same')/np.sum(gauss(x,sig0,mu0)) 


def test():
    mu = 0.6
    sig = 0.2
    sig0 = 0.05
    mu0 = mu
    x = np.linspace( 0.5, 0.6, 10000 )
    plt.plot( x, gauss_log(x, sig, mu, sig0, mu0), '--', label='gauss_X_log', zorder=10 )
    plt.plot( x, gauss(x,sig0,mu0), label='gauss' )
    plt.plot( x, PDF_log(x,sig,mu), label='log' )
    plt.legend()
    plt.show()

这给了我以下结果：

红线是对数正态 PDF，绿线是高斯。 "convolution" 是蓝色虚线。当我更改 x 域时 x = np.linspace( 0.5, 0.8, 10000 ) 我得到非常不同的结果：

显然这里有问题。我的卷积积分结果F(x) = int (g(t)*f(x-t))dt应该不依赖于"x"的范围。然后我将域变大，即 x = np.linspace( 0.00001, 100, 10000 )，这让我胡说八道：

要么我的脚本中存在一个简单的错误，要么我误解了离散卷积。

Answer 1

我发现我错在哪里了：

而不是有一个合适的核函数gauss(x0,sig0,mu)；我对高斯使用了相同的 x。

我认为这是正确的（使用上面的 PDF）：

def gauss_log(x, x0, sig, mu, sig0, mu0):
    a = signal.convolve(gauss(x0,sig0,mu0),PDF_log(x,sig,mu),mode='same')/np.sum(gauss(x0,sig0,mu0))
    return( a )

def test(lightcurve,noisecurve):
    mu = 0.1 #can now put arbitrary small values of mu
    sig = 0.1
    sig0 = 0.05
    mu0 = 0
    x = np.linspace( 0.00001, 5, 1000 )
    x0 = np.linspace(-5,5,1000) #note that arrays need to be equal length!
    g_log = gauss_log(x, x0, sig, mu, sig0, mu0)
    plt.plot( x, g_log, '--', label='gauss_X_log', zorder=10 )
    plt.plot( x0, gauss(x0,sig0,mu0), label='gauss' )
    plt.plot( x, PDF_log(x,sig,mu), label='log' )
    plt.legend()
    plt.show()

    ###testing normalization!
    print(np.trapz(gauss(x0,sig0,mu0),x0))
    print(np.trapz(g_log,x))
    print(np.trapz(PDF_log(x,sig,mu),x))

1.0

0.999938903253

1.0

Answer 2

如果您认为 x 是时间，那么您需要取相同 x 的原因很容易解释。现在，如果你有两个用时间 t 描述的信号，但你没有相同的 x( 或这里的 t) 数组来构造它们；好吧，时间不正确。至少时间 spacing 需要相同 才能以有意义的方式工作。我有一个类似的问题，激光脉冲需要根据时间与浓度进行卷积。激光脉冲非常尖锐，意味着时间也很短，所以我取了一个小间隔 t_laser 而 t_concentration 是我实验的整个时间间隔：

t_laser = np.linspace(-20 ,20,160) 
t_concentration = np.linspace(-2000,20000,88000)
laser = laser_pulse(t_laser) ## Pulshape depending on time
concentration = ABC_model(t_concentration, some_parameters) ## Model for concentration(t)

请注意，对于两个时间数组，我有相同的间距，导致时间步长相等。如果我改成

t_laser = np.linspace(-20 ,20,88000) 
t_concentration = np.linspace(-2000,20000,88000)

因为我想要我的激光脉冲有更高的分辨率，它仍然应该从 t0=-20 到 t1=20。但是由于我只将激光与信号进行卷积，

signal.convolve(concentration, laser, mode='same')/sum(laser)

不考虑时间数组；可以解释为激光脉冲具有相同的时间范围，因为它具有相同的阵列长度。

备注：您还需要了解信号卷积函数的工作原理(scipy.signal.fftconvolve)。在我看来，它将激光阵列的中心放在集中阵列的第一个元素上，然后计算 sum/integral。这基本上就是数组卷积的工作原理（但它不同于 np.convolve；那将是模式“完整”）。因此，如果您选择的高斯函数数组不是对称的，则高斯函数的最大值不在数组的中心（关于索引），因此您的结果会在 x 方向上移动。

scipy 卷积取决于 x

scipy convolve depends on x

python

convolution

scipy