# 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")

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")

## No comments:

## Post a Comment

Note: Only a member of this blog may post a comment.