kd.val <- function(q, kd) {
    # Return the height of a kernel density at the x-value of the density 
    # object closest to the desired quantile
    x <- max(kd$x[kd$x < q])
    x.idx <- which(kd$x == x)
    return(kd$y[x.idx])
}

# Simulate sampling distribution of mean of 10 obs from an Exp(.2) distr.
set.seed(125468)
expmat <- matrix(rexp(100e6, .2), ncol=10)
expmeans <- rowMeans(expmat)

# Plot the exponential and normal densities, and their 95% CIs
lambda  <- .2
mu <- 1 / lambda
se <- sd(expmeans)
xlim <- mu + c(-4, 4)*se
kd.exp <- density(expmeans) 
exp.ci <- quantile(expmeans, c(.025, .975))
norm.ci <- qnorm(c(.025, .975), mu, se)
kd.lcl.y <- kd.val(exp.ci[1], kd.exp)
kd.ucl.y <- kd.val(exp.ci[2], kd.exp)

plot(kd.exp, xlim=xlim, lwd=2, main="Sampling Distributions of the Mean of 10 Observations, and Central 95% Intervals, 
     of an Exponential(.2) Random Variable
     and a Normal Random Variable With the Same Mean and Variance as the Exponential", xlab="X", col="darkblue", bty="l", cex.main=1.2)

curve(dnorm(x, mu, se), add=TRUE, lwd=1, col="blue", lty=2)

lines(x=c(norm.ci[1], norm.ci[1]), y=c(0, dnorm(norm.ci[1], mu, se)), 
      lwd=1, lty=2, col="blue")
lines(x=c(norm.ci[2], norm.ci[2]), y=c(0, dnorm(norm.ci[2], mu, se)), 
      lwd=1, lty=2, col="blue")

lines(x=c(exp.ci[1], exp.ci[1]), y=c(0, kd.lcl.y), lwd=2, col="darkblue")
lines(x=c(exp.ci[2], exp.ci[2]), y=c(0, kd.ucl.y), lwd=2, col="darkblue")

legend(x=8, y=.25, c("Exponential", "Normal"), col=c("darkblue", "blue"),
       lwd=c(2, 1), lty=c(1, 2), cex=1.2)