Hexagon binning plot of bivariate normal data

library(MASS);

library(hexbin);

mu<-c(-1,1);

Sigma<-matrix(0,2,2);

Sigma[1,]<-c(1,0.5); Sigma[2,]<-c(0.5,1);

x<-mvrnorm(10000, mu, Sigma);

plot(hexbin(x[,1], x[,2]), colorcut=seq(0,1,length.out=64),

colramp = function(n) rev(rainbow(64)), legend = 0, xlab="",ylab="" );