How to evaluate a multivariate gaussian condition

Question

I am trying to come up with estimates on the outcomes from a layered-earth inversion calculation whereby the calculation provides a posterior covariance matrix for the best-fitting model parameters given the data. In particular, I want to see how long a particular feature of the model persists from the surface (ie, I am looking for values $m_{i} \ge x_{0}, i=1..M$ . I believe this is a conditional probability question:

1. What is probability layer 1 satisfies condition? $P(m_1 \ge x_0|m_2,...m_M, d_1,d_2,...,d_N)$
2. What is probability layer 2 also satisfies condition? $P(m_2 \ge x_0|m_1 \ge x_0, m_2,...m_M, d_1,d_2,...,d_N)$

... and so on. In this question, my model parameters are correlated through the posterior covariance matrix, and the postulate of the inversion provides me with Gaussian distribution of the parameters. I have been looking around for days how to evaluate an integral of this sort, but the best I could find was for the bivariate case where it was stated that the problem is not closed and must be evaluated numerically. I typically program in Matlab, but am happy to use any method to solve this problem. Is there anyone who can assist with this, please?

Answer 1

In all likelihood, numeric integration is unavoidable.

An effective way to evaluate multivariate integrals is through regular sequences, or quasi Monte Carlo. Methods of this family generate points that fill the space in a very even way, avoiding however putting the points exactly on the grid. The best known method is probably Halton sequence. For a $p$-dimensional version of it, take the first $p$ primes (refer to them as $b_1=2, b_2=3, b_3=5$, etc.), and define reflection $\phi_b(n)$ of a natural number $n$ with respect to the decimal point: if it takes $K$ places to write $n$ as $$ n = \sum_{k=1}^K a_k b^{k-1}, \quad 0\le a_k < b $$ then the reflection is $$ \phi_b(n) = \sum_{k=1}^K a_k b^{-k}. $$ So for instance, in base $3$, we have $11 = 1\cdot 3^2 + 0\cdot 3^1 + 2\cdot 3^0 = 102_3$, which get reflected into $\phi_3(11) = 0.201_3 = 1\cdot3^{-3} + 0\cdot 3^{-2} + 2\cdot 3^{-1} = \frac1{81} + \frac23 = \frac{55}{81} $.

Halton sequence is the sequences of such reflections using the primes: $$\Bigl\{ {\mathbf a}_n = \Bigl(\phi_{b_1}(n), \phi_{b_2}(n), \ldots, \phi_{b_p}(n)\Bigr) : n \in \mathbf{N} \Bigr\}$$

The way it fills up the $p$-dimensional space is that in the $j$-th dimension, among any $b_j$ consecutive numbers, there's exactly one in each of the intervals $[0, 1/b_j), [1/b_j, 2/b_j), \ldots $. Moreover, the regularity properties also hold in the multivariate space, and there is a measure-theoretic sense in which this sequence (or some modifications of it) provide the closest fit to Lebesgue measure of the unit cube.

For Gaussian integration, you would need to take the inverses of the points using the inverse of the standard normal $\Bigl( \Phi^{-1}[\phi_{b_1}(n)], \Phi^{-1}[\phi_{b_2}(n)], \ldots, \Phi^{-1}[\phi_{b_p}(n)]\Bigr) $, which would give you the standard multivariate normal with covariance matrix $I_p$. To get $N(\mu,\Sigma)$ out of it, you need to transform them as $x = \mu + \Sigma^{1/2} z, z \sim N(0,I_p)$, and you can take say Cholesky decomposition of $\Sigma$ as the square root.

I'd be stunned if Matlab, the environment of choice for numeric math, did not have these implemented -- physicists have been working with them for at least three decades, if not longer. In statistics, they were mentioned briefly by Art Owen, and arguably they have greater following in econometrics under the name of GHK algorithm where they are used for numeric integration in choice modeling with more than two alternatives.

Answer 2

There are two ways of thinking about the MVN CDF, which are detailed below.

Mahalanobis Distance

Here, the MVN CDF is the probability that a sample vector, $\vec{x_{0}}$ falls within the ellipsoid defined by your $MVN(\vec{\mu},\Sigma)$ density function. Compute the squared Mahalanobis distance, $d^{2}$ from $\vec{\mu}$ to $\vec{x_{0}}$. Then:

$$F(d^{2}) = P(d^{2} \leqslant \chi_{k}^{2}(\alpha)) = 1-\alpha$$

i.e., the $\operatorname{quantile}(d^{2})=\chi_{k}^{2}(d^{2})$

where $k$ is the dimension of $\vec{\mu}$. Note that this is for the population parameters $\vec{\mu}$ and $\Sigma$. When using sample parameters, see Distribution of an observation-level Mahalanobis distance .

Simultaneous X

The MVN CDF is defined, for a random vector $\vec{x}$, as the probability that every component of $\vec{x}$ is less than a reference vector, $\vec{x_{0}}$:

$$F(\vec{x_{0}}) = P(\vec{x} \leqslant \vec{x_{0}})$$

This can be generalized to the probability of $\vec{x}$ being between two vectors, $\vec{h}$ and $\vec{k}$:

$$F(\vec{h},\vec{k}) = P(\vec{h} \leqslant \vec{x} \leqslant \vec{k})$$

Note that $\vec{h}=-\infty$ and/or $\vec{k}=+\infty$ are possible values. Calculation of this involves the numerical integration the multivariate integral:

$$F_{n}(\vec{h},\vec{k};\vec{\mu} ,\Sigma ) = P(h_{i}<X_{i}<k_{i}; i=1\ldots n) = \int_{h_{1}}^{k_{1}}\int_{h_{2}}^{k_{2}}\cdots \int_{h_{i}}^{k_{i}}f(x_{1},x_{2},\ldots ,x_{i})~dx_{1}dx_{2}\cdots dx_{i}$$

Note that this computes the CDF over a rectangular domain.

The go-to resource for evaluating this is the work of Alan Genz. See his home page for links to journal articles on evaluating this integral and his software page for code in R, MATLAB, and FORTRAN.

MATLAB has the built-in mvncdf() function which computes this integral, with various options for limits.

Special Case

A useful special case for the bivariate normal with correlation coefficient $\rho$ is:

$$F(0,0,\rho) = 1/4 + \frac{1}{2\pi}\arcsin(\rho)$$

This gives you an independent method to determine the probability to check that you are using a software function correctly.

Non-Rectangular Regions

If you need to compute the MVN CDF over some convex region $A$, then see "A Fortran 90 Program for Evaluation of Multivariate Normal and Multivariate t Integrals Over Convex Regions" in Vol. 3, Issue 4 (Feb 1999) of the Journal of Statistical Software.

Should you need to compute it over ellipsoidal regions, see "Numerical Computation of Multivariate Normal and Multivariate t Probabilities over Ellipsoidal Regions" in Vol. 6, Issue 8 (Jun 2001) of the Journal of Statistical Software.

You will find most journal articles for Multivariate Normal CDF will discuss the standardized MVN, where $\vec{\mu} = \vec{0}$ and $\Sigma = I_{n}$. You can standardize your data with the transform $Z=\left | \Sigma \right |^{-1/2}(\vec{x}-\vec{\mu})$.