Search code examples
pythonscipyprobability

Python: How to get the convolution of two continuous distributions?


Let X, Y be 2 random variables, with probability density functions pdf1 and pdf2.

Z = X + Y

Then the probability density function of Z is given by the convolution of pdf1 and pdf2. Since we can't deal with continuous distributions, we descritize the continuous distributions and deal with them.

To find the convolution of uniform distribution and normal distribution, I came up with following code.

import matplotlib.pyplot as plt
import numpy as np
import scipy.stats as stats
from scipy import signal


uniform_dist = stats.uniform(loc=2, scale=3)
std = 0.25
normal_dist = stats.norm(loc=0, scale=std)

delta = 1e-4
big_grid = np.arange(-10,10,delta)

pdf1 = uniform_dist.pdf(big_grid)
print("Integral over uniform pdf: "+str(np.trapz(pdf1, big_grid)))

pdf2 = normal_dist.pdf(big_grid)
print("Integral over normal pdf: "+str(np.trapz(pdf2, big_grid)))


conv_pdf = signal.fftconvolve(pdf1,pdf2,'same')
print("Integral over convoluted pdf: "+str(np.trapz(conv_pdf, big_grid)))

plt.plot(big_grid,pdf1, label='Tophat')
plt.plot(big_grid,pdf2, label='Gaussian error')
plt.plot(big_grid,conv_pdf, label='Sum')
plt.legend(loc='best'), plt.suptitle('PDFs')
plt.show() 

This is the output I get.

Integral over uniform pdf: 0.9999999999976696

Integral over normal pdf: 1.0

Integral over convoluted pdf: 10000.0

If the convolution was correct, I should get a value close to 1 for "Integral over convoluted pdf". So what is going wrong here? Is there a better approach to solve this problem?

Thanks


Solution

  • You should descritize your pdf into probability mass function before the convolution.

    import matplotlib.pyplot as plt
    import numpy as np
    import scipy.stats as stats
    from scipy import signal
    
    
    uniform_dist = stats.uniform(loc=2, scale=3)
    std = 0.25
    normal_dist = stats.norm(loc=0, scale=std)
    
    delta = 1e-4
    big_grid = np.arange(-10,10,delta)
    
    pmf1 = uniform_dist.pdf(big_grid)*delta
    print("Sum of uniform pmf: "+str(sum(pmf1)))
    
    pmf2 = normal_dist.pdf(big_grid)*delta
    print("Sum of normal pmf: "+str(sum(pmf2)))
    
    
    conv_pmf = signal.fftconvolve(pmf1,pmf2,'same')
    print("Sum of convoluted pmf: "+str(sum(conv_pmf)))
    
    pdf1 = pmf1/delta
    pdf2 = pmf2/delta
    conv_pdf = conv_pmf/delta
    print("Integration of convoluted pdf: " + str(np.trapz(conv_pdf, big_grid)))
    
    
    plt.plot(big_grid,pdf1, label='Uniform')
    plt.plot(big_grid,pdf2, label='Gaussian')
    plt.plot(big_grid,conv_pdf, label='Sum')
    plt.legend(loc='best'), plt.suptitle('PDFs')
    plt.show()