Comment tirer des échantillons aléatoires d'une distribution estimée non paramétrique?

14

J'ai un échantillon de 100 points qui sont continus et unidimensionnels. J'ai estimé sa densité non paramétrique en utilisant les méthodes du noyau. Comment puis-je tirer des échantillons aléatoires de cette distribution estimée?

lovekesh
la source

Réponses:

21

Une estimation de la densité du noyau est une distribution de mélange; pour chaque observation, il y a un noyau. Si le noyau est une densité échelonnée, cela conduit à un algorithme simple d'échantillonnage à partir de l'estimation de la densité du noyau:

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

hXjeN(μ=Xje,σ=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))

Si pour une raison quelconque vous ne pouvez pas puiser dans votre noyau (ex. Votre noyau n'est pas une densité), vous pouvez essayer avec un échantillonnage d'importance ou MCMC . Par exemple, en utilisant l'échantillonnage d'importance:

# 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 Avec mes remerciements à Glen_b qui a contribué à la réponse.

Matteo Fasiolo
la source
1
Désolé, je suis entré directement dans l'échantillonnage d'importance, puis j'ai réalisé que l'échantillonnage est généralement plus simple que cela. J'ai ajouté votre explication initiale à mes réponses. Merci beaucoup
Matteo Fasiolo
@ Matteo Fasiolo - Avez-vous une référence à un article que je peux citer pour cette méthode.
Pallavi