SciPy 去卷积函数

Question

我想用SciPy的deconvolve函数来找出一个未知的分布，前提是我有两个高斯分布。不过这个函数在SciPy里没有相关的文档，所以我只是想找个例子，看看这个函数在我的情况下怎么用。比如说，给定两个正态分布N(100, 1)和N(300, 2)，我想知道怎么能找到去卷积后的分布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)

上面的代码给了我一些负值，这看起来不太对。我该怎么从这个去卷积中恢复出分布N(200, 1)呢？特别是，我觉得我的问题在于我不太明白怎么得到除数。

我真正想要的是看到一个例子，展示我如何用SciPy的去卷积从这些样本中恢复出大约N(200, 1)。

Answer 1

我觉得你对自己的期望有点困惑……因为我们都知道，两个正态分布的卷积结果也是一个正态分布，它的均值是两个均值的和，方差是两个方差的和。你似乎期待两个正态随机样本的卷积结果也会是一个正态随机样本，但事实并不是这样：

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)

你可能在想的是这样的情况：

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)

无论如何，回到deconvolve这个话题……如果你用两个长度分别为m和n的数组来调用它，它会返回一个包含两个数组的元组：

• 第一个数组的长度是m - n + 1，这是解卷积后的数组，也就是你需要用第二个数组去卷积，才能得到第一个数组。
• 第二个数组的长度是m，它表示用第二个数组和第一个返回的数组进行卷积时，替换第一个数组所产生的误差。