Question: I'm reading an excellent paper (here) in which the authors begin with a dataset of effect estimates from randomized control trials. In theory, these numbers are a convolution between an unknown distribution and a standard normal. That is to say, each number can be thought of as a draw from the unknown distribution plus some white noise. They claim they can recover the unknown distribution by deconvolving the data with a standard gaussian. I'd like to do this myself in a toy example, but an struggling to obtain sensible results. In the code below I:  Draw 1e5 random numbers from a gamma distribution To each draw, I add white noise I compute a histogram of these new draws, and let the "signal" be the height of the histogram at the midpoint of the bind edges (as defined by numby and matplotlib)  My code to generate the data is as follows 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  Now that I have the signal, I'd like to deconvolve it with a gaussian. The result should hopefully be a gamma distribution (or something close to the gamma density I've used above). 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?

Answer: Yes, when independent random variables X and E are added, the pdf of the sum (X + E) is the convolution of X's pdf and E's pdf. So, if we know the distribution of E and have a histogram estimating the pdf of (X + E), it is possible to obtain a non-parametric estimate of the distribution of X by deconvolution. An important practical limitation to be aware of is that deconvolution is often ill posed, meaning a small change in the input may result in a large change in the deconvolved result. That is a problem here: the histogram is an imperfect estimate of exact pdf of (X + E), and exact deconvolution like scipy.signal.deconvolution will unreasonably amplify those imperfections. To deal with this, it is necessary to use a regularized deconvolution method. There are many ways to do this. A simple classical method is Wiener deconvolution. Here are plots showing results with Wiener deconvolution (code listing below):  The top plot shows our initial estimate of (X + E) from the histogram. The following three plots show the result of Wiener deconvolution with different regularization strengths:  (2nd plot, regularization = 1e-5) If the regularization is too weak, nasty oscillations develop in the output. (3rd plot, regularization = 0.001) It seems that around 0.001 is about right for this problem. (4th plot, regularization = 0.1) If the regularization is too strong, then the deconvolution doesn't do much.  A defect is that the Wiener deconvolution creates negative values, but of course, a pdf should always be nonnegative. A better (but more complicated) method would be to minimize the mean square error objective in the Wiener deconvolution derivation subject to a nonnegativity constraint. And there are other possible techniques beyond this to finesse the regularization for better results. Code: # 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?  The key bit for this is you want to sample the pdf of E over a sequence of points {..., -2*dx, -dx, 0, +dx, +2*dx, ...} where dx is the histogram bin width, and over a range wide enough that most of the area is captured.