# How to draw random samples from a non-parametric estimated distribution?

I have a sample of 100 points which are continuous and one-dimensional. I estimated its non-parametric density using kernel methods. How can I draw random samples from this estimated distribution?

A kernel density estimate is a mixture distribution; for every observation, there's a kernel. If the kernel is a scaled density, this leads to a simple algorithm for sampling from the kernel density estimate:

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


If (for example) you used a Gaussian kernel, your density estimate is a mixture of 100 normals, each centred at one of your sample points and all having standard deviation $h$ equal to the estimated bandwidth. To draw a sample you can just sample with replacement one of your sample points (say $x_i$) and then sample from a $N(\mu = x_i, \sigma = h)$. In 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)  Strictly speaking, given that the mixture's components are equally weighted, you could avoid the sampling with replacement part and simply draw a sample a size$M\$ from each components of the mixture:

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


If for some reason you can't draw from your kernel (ex. your kernel is not a density), you can try with importance sampling or MCMC. For example, using importance sampling:

# 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. With my thanks to Glen_b who contributed to the answer.

• Sorry I went straight into importance sampling, and then I realized that usually sampling is simpler than that. I added your initial explanation to my answers. Many thanks – Matteo Fasiolo Jan 20 '14 at 10:34
• @ Matteo Fasiolo - Do you have any reference to a paper I can cite for this method. – Pallavi May 14 at 23:43