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.
 A: 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))
curve(tmpdens, from=0, to=25, add=FALSE)

curve(dgamma(x,3,1/3), add=TRUE, col='red')

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)
    X <- Gammad(shape,1/rate)
    -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)[1], coef(fit)[2]), 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.
