Restaurar la original distribución de observaciones ruidosas

Question

Hay un conocido conjunto de pares $(y_i, \sigma_i)$ tal que

$y_i = x_i + \sigma_i N_i$

$N_i \sim \mathbf{N}(0,1) $ % todo $i$

$x_i \sim \rho$ % todo $i$

donde $y$ se observa valor, $x$ es el verdadero valor.

¿Cómo restaurar la distribución $\rho$? No es necesariamente normal.

Answer 1

Si el $\sigma_i$ es conocido como un constante entonces esto es solo una variación en la estimación de densidad de kernel. Si no son constantes (pero conocido al menos en una constante), entonces sería una forma de ponderación de estimación de densidad de kernel.

Si el $\sigma_i$ no son conocidos, pero se cree que están relacionados con a $x_i$ a continuación, se hace más complicado, pero usted debe ser capaz de encontrar una estimación razonable el uso de un kernel diferente. Por ejemplo, si es mayor $x_i$ significa un mayor $\sigma_i$, en un determinado $y_i$ serían más propensos a ser de una $x_i$ que es mayor que los más pequeños, por lo que un no-simétrica núcleo sería apropiado. El núcleo sería una convolución de la normal y de la relación entre el$x_i$$\sigma_i$.

Las ideas detrás de registro de spline estimación de densidad de mayo se aplican aquí también. Usted podría utilizar la idea para estimar la densidad de $x$ dada su relación anterior y de máxima verosimilitud.

También hay Bayesiano métodos que utilizan una mezcla de distribuciones de dirichlet para representar y desconocidos de la densidad.

Editar

Aquí hay algunos ejemplos de código R que utiliza el kernel de la estimación de la densidad de intentar reconstruir la densidad original:

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')

Puede suavizar la estimación de la densidad por la multiplicación de los valores de sig por una constante, más constante, más suave.

Y aquí está el código que hace una máxima probabilidad de ajuste a los parámetros al asumir una distribución particular:

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' )

Si usted necesita un no-paramétrico de estimación de la densidad, entonces usted podría utilizar el logspline la estimación de la densidad o de la mezcla de Dirichlets, o de otro tipo no paramétrico de estimaciones en lugar de la paramétrica de la densidad en la anterior.