> 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))
}
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
