CODEWITHSUNDEEP
pythonscipyconvolution
Sebastiano1991
Sebastiano1991

Reputation: 897

scipy convolve depends on x

I am trying to convolve a lognormal PDF and a gaussian PDF. I therefor defined the functions in the following way:

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

This gives me the following result:

first array domain

The red line is the log-normal PDF, the green line the gaussian. The "convolution" is the blue dashed line. When I change the x domain x = np.linspace( 0.5, 0.8, 10000 ) I get very different results:

for different domain

Clearly there is something wrong here. The result of my convolution integral F(x) = int (g(t)*f(x-t))dt should not depend on the range of "x". I then made the domain large, i.e. x = np.linspace( 0.00001, 100, 10000 ), which gives me this nonsense:

large array domain

Either there is a simple bug in my script or I missunderstand the discrete convolution.

Upvotes: 0

Views: 850

Answers (2)

Martin
Martin

Reputation: 113

The reason you need to take the same x is easy explained if you consider x to be time for example. Now if you have two signals that are described with the time t, but you do not have the same x( or here t) array from which you constructed them; well then the times do not correspond correctly. At least the time spacing needs to be the same to work in a way that makes sense. I had a similar problem, where a laser pulse needs to be convolved with a concentration depending on time. The laser pulse is very sharp, meaning also short in time, so I take a small intervall t_laser and t_concentration is the whole time interval of my experiment:

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)

Notice that for both time arrays I have the same spacing, leading to equal steps in time. If I had changed to

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

because I wanted a higher resolution for my laser pulse which should still go from t0=-20 to t1=20. But since I only convolve the laser with the signal,

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

without taking the time arrays into account; one could interpret that the laser pulse has the same temporal extent, since it has the same length of the array.

Remark: Also you need to be aware, how the signal convolve function works (scipy.signal.fftconvolve). It seems to me it places the center of the laser array on the first element of the concentration array and then calculates the sum/integral. That is basically how a convolution with arrays works (But it is different from np.convolve; that would be mode "full"). So if you choose the array for your Gaussian function not to be symmetric, the maximum of the Gaussian is not in the center of the array (regarding the index), so your result shifts in the x direction.

Upvotes: 0

Sebastiano1991
Sebastiano1991

Reputation: 897

I figured out where I had my mistake:

rather than having a proper kernel function gauss(x0,sig0,mu); I used the same x for the Gaussian.

This does what I think is correct (with the PDFs from above):

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

enter image description here

Upvotes: 1

Related Questions