# Calculating 2D Confidence Regions from MCMC Samples

I would like to plot 2D confidence regions (at 1-sigma, 2-sigma) for a model that I've fit to data. I've used PyMC to generate 50k MCMC posterior samples for my model with 6 parameters.

I know the process to create confidence regions is something similar to: 1.) create a histogram of the samples in the 2D space 2.) identify iso-density contours 3.) from a selected start point (eg, the mean) integrate outwards perpendicular to iso-density contours until the desired fraction of sample points are contained in the region.

Is there a convenient function in the numpy/scipy/pymc/pylab/etc world that will create the 2D confidence region plot? Alternatively, where can I find a coded algorithm, or stand-alone tool, that will compute the contours for later plotting?

I once did something like this with pymc, matplotlib, and scipy that you could adapt, the relevant code is in this gist, and the resulting plot looks like this: There is the R function ci2d from gplots that can create 2-dimensional empirical confidence regions: https://rdrr.io/cran/gplots/man/ci2d.html