7
$\begingroup$

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.

Improve this question
$\endgroup$
3
  • 2
    $\begingroup$ 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. $\endgroup$
    – whuber ♦
    Commented Nov 4, 2015 at 13:47
  • $\begingroup$ 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... $\endgroup$
    – Xi'an
    Commented Nov 4, 2015 at 14:13
  • $\begingroup$ @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. $\endgroup$
    – CBowman
    Commented Nov 4, 2015 at 17:41

1 Answer 1

Reset to default
1
$\begingroup$

A Monte Carlo approach offers an easy solution to this problem.

First we need to generate a sample from the PDF - this can be done through interpolation and rejection sampling. Rather than keeping the coordinates of each sample, we need to store the probability density of each sample.

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}$.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.