SciPy deconvolution function

Question

I would like to use SciPy's deconvolve function to find an unknown distribution given two Gaussian distribions. There is no documentation associated with this function in SciPy, so I'm just looking for an example as to how this function can be used in my situation. For example, given two normal distributions N(100, 1), N(300, 2), I would like to understand how I can find the distribution of the deconvolution N(200, 1).

>>> sample1 = np.round(scipy.around(scipy.stats.norm(100, 1).rvs(size=1000)))
>>> sample2 = np.round(scipy.stats.norm(300, 2).rvs(size=2000))
>>> signal.deconvolve(sample1, sample2)

The above code gives me negative values, which seems wrong. How can I recover the distribtion N(200, 1) from this deconvolution? In particular, I think my problem is that I do not understand how to get the divisor.

What I would really like is to see is an example of how I can recover ~ N(200, 1) from these samples using SciPy's deconvolution.

Answer 1

I think you are a little confused about your expectations... Since we all know that the convolution of two normal distributions is another normal distribution with mean the sum of the means, and variance the sum of the variances, you seem to expect that the convolution of two normal random samples will also be a normal random sample. And that just ain't so:

a = scipy.stats.norm(100, 1).rvs(size=1000)
b = scipy.stats.norm(200, 1).rvs(size=1000)
c = scipy.convolve(a, b)
plt.subplot(311)
plt.hist(a, bins=50)
plt.subplot(312)
plt.hist(a, bins=50)
plt.subplot(313)
plt.hist(a, bins=50)

You were probably thinking of something along the lines of:

a = scipy.stats.norm(100, 10).pdf(np.arange(500))
b = scipy.stats.norm(200, 20).pdf(np.arange(500))
c = scipy.convolve(a, b)
m_ = max(a.max(), b.max(), c.max())
plt.subplot(311)
plt.axis([0, 1000, 0, 1.25*m_])
plt.plot(a)
plt.subplot(312)
plt.axis([0, 1000, 0, 1.25*m_])
plt.plot(b)
plt.subplot(313)
plt.axis([0, 1000, 0, 1.25*m_])
plt.plot(c)

In any case, getting back to deconvolve... If you call it with two arrays of length m and n, it will return you a tuple with two arrays:

• the first of length m - n + 1 is the deconvolved array, i.e. the array you should convolve the second one with, to get the first
• the second of length m is the error in replacing the first array with the convolution of the second with the first returned one.