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

Question

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

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

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

• 从 Gamma 分布中提取 1e5 个随机数
• 在每次抽奖中，我都会添加白噪声
• 我计算这些新抽奖的直方图，并且
• 设“信号”为绑定(bind)边中点处直方图的高度(由 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

现在我有了信号，我想用高斯对其进行反卷积。结果应该是 Gamma 分布(或接近我上面使用的 Gamma 密度的东西)。但是，我不确定如何在 scipy.signal.devolve 函数调用中设置“脉冲”。这应该是多长，我应该在什么点评估高斯密度？

Answer 1

是的，当独立随机变量X和E相加时，总和(X + E)的pdf就是X的pdf和E的pdf的卷积。因此，如果我们知道 E 的分布并有一个估计 (X + E) pdf 的直方图，则可以通过反卷积获得 X 分布的非参数估计。

需要注意的一个重要的实际限制是反卷积通常是 ill posed ，这意味着输入的微小变化可能会导致反卷积结果发生较大变化。这是这里的一个问题:直方图是 (X + E) 的精确 pdf 的不完美估计，而像 scipy.signal.devolving 这样的精确反卷积会不合理地放大这些缺陷。为了解决这个问题，需要使用 regularized反卷积方法。有很多方法可以做到这一点。一个简单的经典方法是 Wiener deconvolution .

以下是显示维纳反卷积结果的图(代码列表如下):

上图显示了我们对直方图 (X + E) 的初始估计。以下三幅图显示了不同正则化强度的维纳反卷积结果:

• (第二张图，正则化 = 1e-5)如果正则化太弱，输出中会出现严重的振荡。
• (第三张图，正则化 = 0.001)对于这个问题来说，0.001 左右似乎是合适的。
• (第 4 幅图，正则化 = 0.1)如果正则化太强，则反卷积作用不大。

一个缺陷是维纳反卷积会产生负值，但当然，pdf 应该始终是非负的。更好(但更复杂)的方法是最小化均方误差目标 in the Wiener deconvolution derivation受非负约束。除此之外，还有其他可能的技术来巧妙地进行正则化以获得更好的结果。

代码:

# Copyright 2024 Google LLC.
# SPDX-License-Identifier: Apache-2.0
import numpy as np
from scipy.stats import norm, gamma

desired_dist = gamma(a=4, scale=1/2)  # Underlying desired dist.
noise_dist = norm()                   # Distribution of additive noise.

# Generate the observed samples.
N = int(1e5)
samples = desired_dist.rvs(N) + noise_dist.rvs(N)

# Histogram to estimate (desired + noise) distribution. This is
# the "signal" to be unblurred.
bins = 64
hist, bin_edges = np.histogram(samples, bins=bins, range=[-8, 16])
dx = bin_edges[1] - bin_edges[0]
bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2
signal = hist / (N * dx)

# Sample noise pdf to get the "blur kernel."
kernel = noise_dist.pdf((np.arange(bins) - bins//2) * dx)
kernel /= kernel.sum()

def wiener_deconv(signal, kernel, regularization=0.001):
  signal_f = np.fft.rfft(signal)
  kernel_f = np.fft.rfft(np.fft.fftshift(kernel))
  deconv_f = (signal_f * np.conj(kernel_f)) / (
    np.abs(kernel_f)**2 + regularization)
  return np.fft.irfft(deconv_f, len(signal))

deconv = wiener_deconv(signal, kernel)

However, I'm not sure how to set up the "impulse" in the scipy.signal.deconvolution function call. What length should this be, and at what points should I evaluate the gaussian density?

关键点是您想要在点序列上对 E 的 pdf 进行采样 {..., -2*dx, -dx, 0, +dx, +2*dx, .. .} 其中 dx 是直方图箱宽度，并且范围足够宽，足以捕获大部分区域。