# Calculating the integral of a PDF inside a closed contour of constant density

I'm working with some two-dimensional probability distributions which have emerged from Bayesian inference work I'm doing. These PDFs are stored on regularly spaced Cartesian grids.

I feel like it would be quite common for one to wish to find a contour of constant probability density which encapsulates a certain fraction of the total probability - say one or two sigma worth, or perhaps 95% ect.

From what I can tell, to evaluate this contour we need to be able to calculate the integral of the PDF in the region bounded by a contour. I've written a scheme myself based on splitting the area between two closed contours into a series of triangles to evaluate such an integral, and it seems to work fine.

However - this is such a common problem I feel like there must be a more elegant way of doing this, perhaps exploiting some vector calculus?

If anyone knows something about this I'd be grateful. Feel free to assume that the Cartesian grids are dense enough that interpolation may be used to determine the value of the PDF at an arbitrary point, and that the distributions are uni-modal.

• Since the data are stored on grids, it's hard to see what a "more elegant" method could possibly be. Moreover, the raster format indicates these PDFs might only have been approximated (perhaps with some kind of stochastic algorithm), suggesting that any special effort to achieve high accuracy in this calculation would be wasted. The situation appears to be identical to the one discussed at Integrating kernel density estimator in 2D. – whuber Nov 4 '15 at 13:47
• This is a pretty hard numerical problem that I [in a related research project] concluded could not be solved nicely. So I use instead simulation and a rough convexification for finding HPD regions... – Xi'an Nov 4 '15 at 14:13
• @whuber It's not really a case of desiring high accuracy, It was more about simplicity - the scheme I'd originally come up would have been a nightmare to explain on paper. However I've thought about it more now, and I think I have a much nicer approach based on transforming the integration problem such that one of the integrals is a closed loop integral around the contours. I'll post it as an answer once it's tested. – CBowman Nov 4 '15 at 17:41

Let $z$ be the density of the contour inside which we want to calculate the integral. If $N$ is the total number of samples, and $m$ is the number of samples which have a density value greater than $z$, then
$$\idotsint\limits_{P(\underline{x}) \ge z} P(\underline{x}) \; \mathrm{d}\underline{x} \approx \frac{m}{N}$$
The number of samples can be increased to achieve whatever accuracy is required, as the error in the calculation is proportional to $1 / \sqrt{N}$.