N <- 10
x <- rgamma(N, 1, 1)
plot(density(x), col="red")
rug(x)
for (i in 1:N) {
   lines(seq(-1, 4, 0.01), 0.1 * dnorm(seq(-1, 4, 0.01), x[i], density(x)$bw))
}