Computing confidence region for Gaussian mixture model I have a 2-d Gaussian mixture model and would like to compute a confidence region for it. Our application is that the two dimensions are latitude and longitude; that is, we want to say something like "the model is 95% confident that the true point lies within this (perhaps noncontiguous) region".
In other words, I think what we want is to find the contour which contains 95% (or whatever) of the probability volume in the model.
My questions:


*

*Is there a simple and fast way to do this?

*Are there any libraries which can do this, ideally in Python?


We need the region as a (multi)polygon object which we can then pass around the rest of our software, not a plot.
 A: In general, it is possible for instance to compute what is the probability contained in a ball $\mathcal{B}(c,r)$. I suppose that your gaussian mixture writes
$$p(x) = \sum_{j=1}^K \mathcal{N}(x;\mu_j,\Sigma_j)\mathbb{P}(J=j)$$
There exist elementary pieces of code to make an algorithm to compute $F(c,r)=\int_{\mathcal{B}(c,r)} p(x)dx $ for each $(c,r)\in \mathbb{R}^N \times \mathbb{R}^+$ which I will detail further. First note that since for a fixed $c$, the function $r\mapsto F(c,r)$ is increasing on $\mathbb{R}^+$ then a search by dichotomy can numerically solve the problem of finding $r$ such as, for a given $c$, $F(c,r)=95\%$. More efficient methods such as the secant method exist.
To handle the general case, if the matrices $\Sigma_j$ have no particular form except that they are positive definite, you may have to use [wikipedia link to the generalized non central chi square cumulative function][1]. For this purpose, the C source code will be useful [source code from Robert Davies' website][2]. You'll find the related documentation [on Robert Davies' page][3] and in his paper [The Distribution of a Linear Combination of Chi Squared Random Variables][4] in which he adopts the same symbols.
If the matrices $\Sigma_j$ are each one proportional to the identity matrix, you may use either the already mentioned generalized non central chi square cumulative function or the non central chi square cumulative function which is more common (see [Non Central Chi Square Law][5]). This function is available in MATLAB for instance.
Now, here is how you can use it.
$$p(x) = \sum_{j=1}^K \mathcal{N}(x;\mu_j,\Sigma_j)\mathbb{P}(J=j)$$
is the density of a variable that we call $X_J$ where for each $i \in \{1,...,K\}$, the variable $X_i$ is a gaussian random variable which has a probability density function given by $\mathcal{N}(x;\mu_j,\Sigma_j)$ and where $J$ is a random discrete variable on the set $\{1,...,K\}$ independent of each $X_i$, and which follows the known law $\mathbb{P}(J=j)$.
Since we can decompose the probability
$$\int_{\mathcal{B}(c,r)} p(x)dx = \mathbb{P}\left(X_J \in \mathcal{B}(c,r) \right) = \sum_{j=1}^K \mathbb{P}\left(X_J \in \mathcal{B}(c,r) |J=j\right)\mathbb{P}(J=j)$$
What we need to compute is $\mathbb{P} \left(X_J \in \mathcal{B}(c,r) |J=j\right)=\mathbb{P}\left(\|X_j-c\|^2 \leq r^2 \right)$.
Since $X_j-c$ is a gaussian random variable which follows $\mathcal{N}(x;\mu_j-c,\Sigma_j)$, then $\|X_j-c\|^2$ follows a (generalized non central) chi squared law.
Some derivations must be done to identify the parameters (which we call $\theta_j$) of this law. (I can be more explicit on demand). 
Then the only thing remaining to do is to evaluate the Chi Squared cumulative function (which we call $S$) of this chi squared law at $r$.
Finally : 
$$F(c,r)=\int_{\mathcal{B}(c,r)} p(x)dx =
\sum_{j=1}^K S(r;\theta_j)\mathbb{P}(J=j) $$
Then, one can apply a dichotomy or a secant method to find the better approximation to $r$ which ensusres $\mathcal{B}(c,r)$ to contain 95%.

If you are in the particular case where the matrices $\Sigma_j$ are diagonal, then you can find a rectangular region. By rectangle I mean a domain which is a cartesian product of intervals. You'll need the erf function, which is related to the cumulative function of a gaussian probability density function.

To answer another question posted : The union of the contours contains exacly or more than 95% if each contour contains 95% of probability.
Here is why. Let $E_i$ be the contour such as $\mathbb{P}(X_i \in E_i)=95\%$ and let $\bigcup_{i=1}^K  E_i $ be the union of the contours then
$$ \int_{ \bigcup_{i=1}^K E_i} p(x)dx = \mathbb{P}\left(X_J \in \bigcup_{i=1}^K E_i\right)$$
$$= \sum_{i=1}^K \mathbb{P}\left(X_J \in  \bigcup_{i=1}^K E_i |J=j\right)\mathbb{P}(J=j)$$
in which each term $\mathbb{P}\left(X_J \in \bigcup_{i=1}^K E_i \Big| J=j \right) \geq 95\%$. Finally, because $\mathbb{P}(J=j)$ sums to 1, this last line is a weighted average of values greater than 95%, thus the sum is greater than 95%.
