I'm coming from [this question][1] in case anybody wants to follow the trail.

Basically I have a data set $\Omega$ composed of $N$ objects where each object has a given number of measured values attached to it (two in this case):

$$\Omega = o_1[x_1, y_1], o_2[x_2, y_2], ..., o_N[x_N, y_N]$$

I need a way to determine the probability of a **new** object $p[x_p, y_p]$ of belonging to $\Omega$ so I was advised in that question to obtain a probability density $\hat{f}$ through a kernel density estimator, which I believe I already have. 

Since my goal is to obtain the probability of this new object ($p[x_p, y_p]$) of belonging to this 2D data set $\Omega$, I was told to integrate the pdf $\hat{f}$ over "_values of the support for which the density is less than the one you observed_". The "observed" density is $\hat{f}$ evaluated in the new object $p$, ie: $\hat{f}(x_p, y_p)$. So I need to solve the equation:


$$\iint_{x, y:\hat{f}(x, y) < \hat{f}(x_p, y_p)} \hat{f}(x,y)\,dx\,dy$$


Here comes the tricky part, my 2D data set and my pdf (I calculated it through `python`'s [stats.gaussian_kde][2] module) look like this:

![enter image description here][3]

where the red dot represents the new object $p[x_p, y_p]$ plotted over my existing data set.

So the question is: how can I calculate the above integral for the limits $x, y:\hat{f}(x, y) < \hat{f}(x_p, y_p)$ when the pdf looks like that?

It is very likely that I could have misunderstood something in the previous question (my training in statistics is almost null), so if anything makes no sense please do let me know.

  [1]: http://stats.stackexchange.com/questions/63263/probability-of-an-object-of-belonging-to-a-group-of-objects
  [2]: http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.gaussian_kde.html
  [3]: https://i.sstatic.net/RnIDI.png