The distribution of the sum of two normal independent random variables can be computed with aconvolution. I.e. Z=X+Y and f_z = convolution(f_x, f_y).
I know that the convolution of two normal distributions is also a normal distribution with mean mu_1 + mu_2 and variance var_1 + var_2. Now, I have a normal distribution Z = X + Y with N(0, std=5) and I also have the distribution X with N(20, std=4). I apply deconvolution to find the distribution Y, which should be N(-20, std=3). However when I apply scipy.signal.deconvolve(Z, X) it doesn't work as expected.
import numpy as np
import scipy.signal
from matplotlib.pyplot import draw, figure, show
dx = 1e-1
grid = np.arange(-20, 30, dx)
dist_X = scipy.stats.norm(loc=20, scale=4)
dist_Z = scipy.stats.norm(loc=0, scale=5)
# Dist Y should be Gaussian : N(loc=-20, scale=np.sqrt(25-16))
Z_signal = dist_Z.pdf(grid)
X_signal = dist_X.pdf(grid)
# Reconstructed distribution
Y_signal = scipy.signal.deconvolve(Z_signal, X_signal)
plt.figure(figsize=(8, 6))
plt.plot(grid, Y_signal[1]/Y_signal[0], label=r'$Y$')
plt.plot(grid, X_signal, label=r'$X$')
plt.plot(grid, Z_signal, label=r'$Z$')
plt.legend(fontsize='16')
plt.grid(axis='x', color='0.95')
plt.grid(axis='y', color='0.95')
plt.ylabel(r'PDF', fontsize='20')
plt.xlabel(r'Distribution', fontsize='20')
draw()
show()
This is the result
Should I use another method? or am I missing something about this deconvolve tool? This case is pretty easy, but I intend to apply this method to other probability distributions which are not Gaussian.
