Sample from an estimated probability density in R

sample_from_density.R

## Sample from an estimated probability density
sample.density <- function(dd, size, jitter=FALSE) {
  x <- sample(dd$x, size=size, prob=dd$y, replace=TRUE)
  if (dd$n < 2) jitter=FALSE
  if (jitter) {
    w <- mean(dd$x[2:length(dd$x)] - dd$x[1:(length(dd$x)-1)]) / 2.0
    x <- x + runif(size,-w,w)
  }
  return(x)
}

## pretend this is measured data and we don't know the PDF
dat <- c(rnorm(40), rnorm(40,3,2), runif(20,0,2))
dd <- density(dat)
plot(dd)
for (i in 1:100) {
  s1 <- sample.density(dd, size=1000, jitter=TRUE)
  lines(density(s1), col="#33669980")
}
lines(dd, col="blue")