2dhist()

n1 <- 500
n2 <- 300
n3 <- 200
x <- c(rnorm(n1, 0, 0.5), rnorm(n2, 5, 1), rnorm(n3, 8, 2))
y <- c(rnorm(n1, 0, 2), rnorm(n2, 3, 2), rnorm(n3, -3, 1))
library(gregmisc)
plot(x, y)
h2d <- hist2d(x, y, show = FALSE, same.scale = TRUE, nbins = c(10,10))
filled.contour(h2d$x, h2d$y, h2d$counts, nlevels = 10, col = gray((10:0)/10))
persp(h2d$x, h2d$y, h2d$counts, ticktype = "detailed", theta = 60,phi = 30, shade = 0.5, col = "cyan")