Contour plot where levels include quantiles instead of levels of constant density

QuantileLevelSetsContours.R

library(MASS)
xdat = rnorm(1000)
ydat = rnorm(1000)

# https://stackoverflow.com/questions/16225530/contours-of-percentiles-on-level-plot
dens = kde2d(xdat, ydat, n=200); ## estimate the z counts

prob = c(0.95, 0.5)
dx = diff(dens$x[1:2])
dy = diff(dens$y[1:2])
sz = sort(dens$z)
c1 = cumsum(sz) * dx * dy
levs = sapply(prob, function(x) {
  approx(c1, sz, xout = 1 - x)$y
})

plot(xdat,ydat, pch=19, col="red", cex.lab=2,cex.main=1.0,cex.axis=1.75)    
contour(dens, levels=levs, labels=prob, add=T,lwd=2)