[R] Contour plot (level curves)

[Eric Elguero]

> I have a sample of n values from a bivariate distribution (from a MCMC
> procedure). How could I draw a contour plot of "the joint density" based on
> that sample ?

here is a fast 2D density estimator. Not very sophisticated, but works.
The function assumes that data are in the form of a matrix
with (first) two columns containing x and y coordinates.

To plot the result:
image(dens2d(x)) or contour(dens2d(x))

Play with the h parameter to change the smoothness of the surface.


>dens2d
function(x, nx = 20, ny = 20, margin = 0.05, h = 1)
{
 xrange <- max(x[, 1]) - min(x[, 1])
 yrange <- max(x[, 2]) - min(x[, 2])
 xmin <- min(x[, 1]) - xrange * margin
 xmax <- max(x[, 1]) + xrange * margin
 ymin <- min(x[, 2]) - yrange * margin
 ymax <- max(x[, 2]) + yrange * margin
 xstep <- (xmax - xmin)/(nx - 1)
 ystep <- (ymax - ymin)/(ny - 1)
 xx <- xmin + (0:(nx - 1)) * xstep
 yy <- ymin + (0:(ny - 1)) * ystep
 g <- matrix(0, ncol = nx, nrow = ny)
 n <- dim(x)[[1]]
 for(i in 1:n) {
  coefx <- dnorm(xx - x[i, 1], mean = 0, sd = h)
  coefy <- dnorm(yy - x[i, 2], mean = 0, sd = h)
  g <- g + coefx %*% t(coefy)/n
 }
 return(list(x = xx, y = yy, z = g))
}