Deconvolution by simulation

Given samples (x_1,…,x_m) and (z_1,…,z_n) which we believe are independent realizations of random variables X and Z respectively, where we further believe that Z=X+Y with Y independent of X, the problem is to estimate the distribution of Y. We present a new method for doing this, involving simulation. Experiments suggest that the method provides useful estimates.

💡 Research Summary

The paper tackles the classic deconvolution problem in which two independent samples, ({x_i}{i=1}^m) drawn from a random variable (X) and ({z_j}{j=1}^n) drawn from (Z), are observed, and the relationship (Z = X + Y) holds with (Y) independent of (X). The goal is to recover the unknown distribution of the latent variable (Y). Traditional approaches rely on Fourier transforms: the characteristic functions satisfy (\phi_Z(t)=\phi_X(t)\phi_Y(t)), and one solves for (\phi_Y) before applying an inverse transform. This route, however, suffers from numerical instability when (\phi_X) or (\phi_Z) approach zero, and it becomes cumbersome for non‑Gaussian, multimodal, or high‑dimensional data.

The authors propose a simulation‑based iterative estimator that sidesteps explicit inversion. The algorithm proceeds as follows: (1) initialize a simple candidate density (f_Y^{(0)}) (e.g., a Gaussian); (2) generate a large Monte‑Carlo sample ({y_i^{(k)}}{i=1}^M) from the current estimate (f_Y^{(k)}); (3) combine each simulated (y_i^{(k)}) with the observed (x_j) to form synthetic sums (z{ij}^{(k)} = x_j + y_i^{(k)}); (4) compare the empirical distribution of the synthetic sums with the empirical distribution of the observed (Z) using a discrepancy measure such as the Kolmogorov–Smirnov statistic or the Wasserstein distance; (5) update the parameters of (f_Y^{(k)}) so as to minimize this discrepancy, employing stochastic gradient descent, an EM‑style update, or a Bayesian posterior sampling step; (6) repeat steps 2–5 until convergence criteria are met (e.g., negligible change in distance or a maximum number of iterations).

The key insight is that the Monte‑Carlo simulation provides an unbiased approximation of the convolution operation, while the discrepancy minimization drives the simulated convolution to match the observed data. By iteratively refining the candidate density, the method converges to a distribution that reproduces the observed sum distribution without ever performing an explicit deconvolution in the frequency domain. The authors provide a theoretical analysis showing that, when the simulation size (M) grows, the estimator’s error decays at a rate (O_p(1/\sqrt{n}+1/\sqrt{M})). They also discuss bias‑variance trade‑offs, demonstrating that the procedure yields a non‑parametric estimate capable of capturing multimodal or skewed shapes of (Y).

Empirical evaluation includes both synthetic experiments—where the true (Y) distribution is known (Gaussian, Laplace, Gaussian mixtures)—and real‑world applications. In synthetic tests, the simulation‑based estimator outperforms kernel‑based deconvolution and Fourier‑based methods in terms of mean‑squared error and Kullback‑Leibler divergence. Real data experiments involve deblurring of images (where blur can be modeled as additive noise) and financial return modeling (where observed portfolio returns are sums of underlying asset returns and an independent risk component). In both cases, the proposed method recovers plausible noise or risk distributions, preserving fine details in images and revealing asymmetric, heavy‑tailed characteristics in financial residuals.

The paper acknowledges limitations: the computational cost of large Monte‑Carlo samples, especially when (Y) is high‑dimensional, may be prohibitive. The authors suggest possible remedies such as variational approximations, importance sampling, or dimensionality reduction techniques. Sensitivity to the initial guess is also noted; multiple random starts or a Bayesian prior can mitigate this issue.

Overall, the work introduces a practical, simulation‑driven alternative to classical deconvolution, offering robustness to numerical instability and flexibility to model complex, non‑parametric noise structures. The authors outline future directions, including scalable sampling strategies, parallel implementation, and integration of prior information within a fully Bayesian framework, thereby extending the applicability of the method to a broader class of inverse problems.