How do you get N pairs of values, which represent a joint probability (2d density) of a much larger pairs of values?
I do MCMC sampling on parameters of a function, and I want to visualize the posterior density of that function by plotting, say, 20 semi-transparent lines that visualizes the density of the function. I know that I can just make a sufficiently large sample to approximate the density (like this), but this could clutter the graph. Rather, I'd imagine something like percentiles would work. E.g. a pair for each 2% change in percentile should accurately portray the density using 50 pairs.
Specifically, I'm sampling from a bayesian model with a two-parameter geometric series. The density is not necessarily monotonically increasing (there can be several peaks). This is a bit complicated, but just to get started, a multivariate normal would be more minimal to work with:
library(MASS)
pairs = mvrnorm(10000, c(1,3), rbind(c(0.2, 0.1), c(0.1, 0.2)))
# Easy to just sample the rows to get an approximate representation:
pairs[sample(nrow(pairs), size=50, replace=FALSE)
# But how to get more certainty that the samples are really representative?
# This would probably start with the density
dens = kde2d(pairs[,1], pairs[,2], n=100)
Here, dens$z is a 100*100 matrix containing a density for each of the combinations of dens$x and dens$y, i.e. the (binned) pairs in pairs. Here's an image(dens) visualization:
