Como extrair amostras aleatórias de uma distribuição estimada não paramétrica?

Eu tenho uma amostra de 100 pontos que são contínuos e unidimensionais. Estimei sua densidade não paramétrica usando métodos de kernel. Como posso tirar amostras aleatórias dessa distribuição estimada?

Uma estimativa de densidade de kernel é uma distribuição de mistura; para cada observação, há um núcleo. Se o kernel é uma densidade escalada, isso leva a um algoritmo simples para amostragem da estimativa de densidade do kernel:

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

Se (por exemplo) você usou um kernel Gaussiano, sua estimativa de densidade é uma mistura de 100 normais, cada uma centrada em um dos pontos de amostra e todas com desvio padrão igual à largura de banda estimada. Para desenhar uma amostra, você pode apenas amostrar com um dos seus pontos de amostra (por exemplo, x i ) e depois amostrar a partir de um N ( μ = x i , σ = h ) . Em 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)

$M$

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

Se, por algum motivo, você não conseguir extrair do seu kernel (por exemplo, o kernel não é uma densidade), tente com amostragem importante ou MCMC . Por exemplo, usando amostragem de importância:

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

PS Com meus agradecimentos a Glen_b, que contribuiu para a resposta.