# Restoring original distribution from noisy observations

There is a known set of pairs $(y_i, \sigma_i)$ such that

$y_i = x_i + \sigma_i N_i$

$N_i \sim \mathbf{N}(0,1)$ for all $i$

$x_i \sim \rho$ for all $i$

where $y$ is observed value, $x$ is true value.

How to restore the distribution $\rho$ ? It is not neccessarily normal.

• When you say "restore", what do you actually mean? Do you know the functional form of $\rho$ and seek to estimate the parameters? Do you seek a nonparametric estimate of $\rho$? Or something else? In addition, how large a sample do you have? – jbowman Jul 29 '14 at 18:28
• It sounds like you're after something like deconvolution, but it's very hard to tell for sure. Responses to jbowman's questions will make for better answers. – Glen_b -Reinstate Monica Jul 29 '14 at 22:52
• I seek a nonparametric estimate of $\rho$ – user31264 Jul 30 '14 at 19:41

If the $\sigma_i$ are a known constant then this is just a variation on kernel density estimation. If they are not constant (but known at least up to a constant) then it would be a form of weighted kernel density estimation.

If the $\sigma_i$ are not known, but believed to be related to $x_i$ then it becomes more complicated, but you should still be able to find a reasonable estimate using a different kernel. For example if larger $x_i$ means a larger $\sigma_i$ then a given $y_i$ would be more likely to be from an $x_i$ that is larger than smaller, so a non-symmetric kernel would be appropriate. The kernel would be a convolution of the normal and the relationship between $x_i$ and $\sigma_i$.

The ideas behind log spline density estimation may apply here as well. You could use that idea to estimate the density of $x$ given your relationship above and maximum likelihood.

There are also Bayesian methods that use a mixture of dirichlet distributions to represent and unknown density.

Edit

Here is some example R code that uses the kernel density estimate to try to reconstruct the original density:

x <- rgamma(100, 3, 1/3)
e <- rnorm(100)
sig <- runif(100, 0.5, 1.5)
y <- x + sig*e

plot( density(y, kernel='gaussian', bw=1, weights=sig/sum(sig)), type='l' )

densgen <- function(y, sig) {
n <- length(y)
function(x) {
sum( dnorm(x, y, sig)/n )
}
}

tmpdens <- Vectorize(densgen(y, sig))



You can smooth the density estimate by multiplying the values of sig by a constant, the bigger the constant, the more smooth.

And here is some code that does a maximum likelihood fit to the parameters when assuming a particular distribution:

library(distr)
library(stats4)

ll <- function(shape, rate) {
if( shape <= 0 || rate <= 0 ) return(Inf)
-sum( sapply( seq_along(y), function(i) {
E <- Norm(0,sig[i])
Y <- X + E
log( d(Y)(y[i]) )
} ) )
}

fit <- mle( ll,  start = list( shape=3, rate=1/3 ) )

curve( dgamma(x, coef(fit), coef(fit)), add=TRUE, col='blue' )


If you need a non-parametric estimate of the density, then you could use the logspline density estimate or mixture of Dirichlets, or other non parametric estimates in place of the parametric density in the above.

• Each $\sigma_i$ is known (and they are all different), but it is not very helpful that the problem is "just a variation on kernel density estimation". – user31264 Jul 30 '14 at 19:49
• @user31264, I added some R code in my example to demonstrate my suggested approaches. – Greg Snow Jul 31 '14 at 16:58
• @user31264, I just came across and article that based on the abstract looks like it may be what you are looking for (I have not read the full article yet). The article is: Deconvolving kernel density estimators; Statistics: A Journal of Theoretical and Applied Statistics. Volume 21, Issue 2, 1990. pages 169-184. – Greg Snow Aug 1 '14 at 16:10