¿Cómo extraer las muestras aleatorias de una distribución no paramétrica estimada?

Question

Tengo una muestra de 100 puntos que son continuas y unidimensional. Calcula su densidad no paramétricos utilizando métodos kernel. ¿Cómo puedo sacar muestras al azar de esta distribución estimada?

Answer 1

Una estimación de densidad de kernel es una mezcla de distribución; para cada observación, hay un núcleo. Si el núcleo es un modelo a escala de densidad, esto conduce a un algoritmo simple para el muestreo a partir de la estimación de densidad de kernel:

repeat nsim times:
  sample (with replacement) a random observation from the data
  sample from the kernel, and add the previously sampled random observation

Si (por ejemplo) que utiliza un núcleo Gaussiano, su estimación de la densidad es una mezcla de 100 normales, cada uno centrado en uno de sus puntos de muestra y todos con desviación estándar $h$ igual a la estimación de ancho de banda. Para dibujar una muestra que usted puede simplemente muestra con reemplazo de uno de sus puntos de muestra (decir $x_i$) y, a continuación, muestra de una $N(\mu = x_i, \sigma = h)$. En R:

# Original distribution is exp(rate = 5)
N = 1000
x <- rexp(N, rate = 5)

hist(x, prob = TRUE)
lines(density(x))

# Store the bandwith of the estimated KDE
bw <- density(x)$bw

# Draw from the sample and then from the kernel
means <- sample(x, N, replace = TRUE)
hist(rnorm(N, mean = means, sd = bw), prob = TRUE)

Estrictamente hablando, dado que la mezcla de los componentes son de igual peso, usted podría evitar el muestreo con reemplazo de parte y simplemente sacar una muestra de un tamaño de $M$ a partir de cada uno de los componentes de la mezcla:

M = 10
hist(rnorm(N * M, mean = x, sd = bw))

Si por alguna razón usted no puede sacar de su núcleo (ex. su núcleo no es una densidad), puedes probar con la importancia de muestreo o MCMC. Por ejemplo, usando la importancia de muestreo:

# Draw from proposal distribution which is normal(mu, sd = 1)
sam <- rnorm(N, mean(x), 1)

# Weight the sample using ratio of target and proposal densities
w <- sapply(sam, function(input) sum(dnorm(input, mean = x, sd = bw)) / 
                                 dnorm(input, mean(x), 1))

# Resample according to the weights to obtain an un-weighted sample
finalSample <- sample(sam, N, replace = TRUE, prob = w)

hist(finalSample, prob = TRUE)

P. S. Con mi agradecimiento a Glen_b que han contribuido a la respuesta.