对直方图定义的分布进行反卷积

Question

我正在阅读一篇优秀的论文（这里），其中作者从随机对照试验的效果估计数据集开始。

理论上，这些数字是未知分布和标准正态分布之间的卷积。也就是说，每个数字都可以被认为是从未知分布加上一些白噪声得出的结果。

他们声称可以通过用标准高斯对数据进行反卷积来恢复未知分布。我想在一个玩具示例中自己做这件事，但很难获得合理的结果。在下面的代码中我：

• 从伽玛分布中绘制

  1e5

  随机数
• 在每次抽奖中，我都会添加白噪声
• 我计算这些新抽奖的直方图，并且
• 令“信号”为绑定边中点处直方图的高度（由 numby 和 matplotlib 定义）

我生成数据的代码如下

import numpy as np
from scipy.stats import norm, gamma
from scipy.signal import convolve, deconvolve
import matplotlib.pyplot as plt

# First, create the original signal
N = int(1e5)
X = gamma(a=4, scale=1/2).rvs(N)

# Corrupt with gaussian noise
E = norm().rvs(N)
Y = X + E

height, edge, ax = plt.hist(X, edgecolor='white', bins = 50);  
mid = (edge[:-1] + edge[1:])/2

现在我有了信号，我想用高斯对其进行反卷积。结果应该是伽玛分布（或接近我上面使用的伽玛密度的东西）。但是，我不确定如何在

scipy.signal.deconvolution

函数调用中设置“脉冲”。这应该是多长，我应该在什么点评估高斯密度？

Answer 1

让我知道这个答案是否有任何帮助。如果不是我会删除它。

scipy.signal.deconvolve

函数专为一维反卷积而设计，您需要提供脉冲响应。在您的情况下，脉冲响应是标准正态分布。您应该通过离散高斯分布来定义脉冲响应。

以下是执行反卷积的方法：

# Create the original signal
N = int(1e5)
X = gamma(a=4, scale=1/2).rvs(N)

# Corrupt with Gaussian noise
E = norm().rvs(N)
Y = X + E

# Define the impulse response (standard normal distribution)
impulse_response = norm().pdf(np.linspace(-5, 5, 1000))  # Discretized Gaussian

# Deconvolve the signal
deconvolved, remainder = deconvolve(Y, impulse_response)

# Plot the original signal, corrupted signal, and deconvolved signal
plt.figure(figsize=(12, 4))

plt.subplot(1, 3, 1)
plt.hist(X, edgecolor='white', bins=50)
plt.title('Original Signal')

plt.subplot(1, 3, 2)
plt.hist(Y, edgecolor='white', bins=50)
plt.title('Corrupted Signal')

plt.subplot(1, 3, 3)
plt.plot(deconvolved)
plt.title('Deconvolved Signal')

plt.show()

在此代码中，

impulse_response

是离散高斯分布。您可以根据您的喜好调整离散化的范围和分辨率。然后

deconvolve

函数执行反卷积。

请记住，反卷积可能是一个不适定问题，结果可能并不完美，尤其是在噪声水平很高的情况下。根据数据的特征，可能需要调整参数并尝试不同的方法。