RG#93: Add countour or heat map plot to XY scatter plot

# data
set.seed(1234)
n <- 10000
X = rnorm (n, 10, 4)
Y = X*1.5 + rnorm (n, 0, 8)

## colour brewing
library(RColorBrewer)
g = 11
my.cols <- rev(brewer.pal(g, "RdYlBu"))
#compute 2D kernel density

# kernel density using MASS 
library(MASS)
z <- kde2d(X, Y, n=50)
plot(X, Y, xlab="X", ylab="Y", pch=19, cex=.3, col = "gray60")
contour(z, drawlabels=FALSE, nlevels=g, col=my.cols, add=TRUE, lwd = 2)
abline(h=mean(Y), v=mean(X), lwd=2, col = "black")
legend("topleft", paste("r=", round(cor(X, Y),2)), bty="n")

## estimate the z counts
prob <- c(.99, .95, .90, .8, .5, .1, 0.05)
dx <- diff(z$x[1:2])
dy <- diff(z$y[1:2])
sz <- sort(z$z)

c1 <- cumsum(sz) * dx * dy

levels <- sapply(prob, function(x) {
              approx(c1, sz, xout = 1 - x)$y })

plot(X,Y, col = "gray80", pch = 19, cex = 0.3)
contour(z, levels= round (levels,7), add=T, col = "red")

# smooth scatter
require(KernSmooth)
smoothScatter(X, Y, nrpoints=.3*n, colramp=colorRampPalette(my.cols), pch=19, cex=.3, col = "green1")