Python: How to get the convolution of two continuo

问题:

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

回答1:

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

回答2:

Apart from discretizing, this does not seem to be currently possible with scipy.stats continuous distributions, since the convolution gives rise to unique distributions. If one density function is Gaussian and the other is uniform, their convolution is a 'blurred gaussian'. This is neither Gaussian nor uniform.

There are some useful particular cases, however. If you are dealing with normal distributions, for instance, the convolution of two independent distributions will be also normal.

Just one detail not stressed in your question - the convolution formula only holds if X and Y are independent.

回答3:

To make this work with discretized pdf's you need to normalize the output of fftconvolve:

conv_pdf = signal.fftconvolve(pdf1, pdf2, 'same') * delta

Note that fftconvolve cannot do it by itself since it doesn't know the actual pdf's, only the values.