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

Tengo una muestra de 100 puntos que son continuos y unidimensionales. Calculé su densidad no paramétrica utilizando métodos del núcleo. ¿Cómo puedo extraer muestras aleatorias de esta distribución estimada?

Una estimación de la densidad del grano es una distribución de mezcla; por cada observación, hay un núcleo. Si el núcleo es una densidad escalada, esto lleva a un algoritmo simple para el muestreo de la estimación de densidad del núcleo:

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) usó un núcleo gaussiano, su estimación de densidad es una mezcla de 100 normales, cada una centrada en uno de sus puntos de muestra y todas tienen una desviación estándar igual al ancho de banda estimado. Para extraer una muestra, puede simplemente tomar una muestra con el reemplazo de uno de sus puntos de muestra (digamos x i ) y luego tomar muestras de una N ( μ = x i , σ = h ) . En R:$h$$x_i$$N(\mu = x_i, \sigma = h)$

# 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 los componentes de la mezcla están igualmente ponderados, puede evitar el muestreo con una pieza de repuesto y simplemente extraer una muestra de un tamaño de cada componente de la mezcla:$M$

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

Si por alguna razón no puede extraer de su núcleo (por ejemplo, su núcleo no es una densidad), puede intentar con muestreo importante o MCMC . Por ejemplo, usando muestreo de importancia:

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

PD: Gracias a Glen_b, que contribuyó a la respuesta.