Numerical computation of the means and covariance in a truncated bivariate normal distribution

Question

How can I compute the means and covariance of a truncated bivariate normal distribution? I am particularly worried about the case when the truncation occurs very far from the mean. Is there a robust numerical evaluation procedure?

By truncated bivariate normal distribution, I mean a density of the form:

$$f(\vec{x}) \propto \exp \left\{ -\left(\vec x - \vec \mu\right)^T \Sigma^{-1} \left(\vec x - \vec \mu\right) \right\}$$

for $\vec a \le \vec x \le \vec b$ (component-wise inequalities), while $f(\vec x) = 0$ othewise. Here $\vec x, \vec \mu, \vec a, \vec b$ are two-dimensional real vectors and $\Sigma$ is a $2\times2$ invertible real symmetric matrix that is positive semi-definite. The implicit proportionality constant guarantees normalization over the truncation rectangle.

I want to compute the means and covariance in this distribution. Explicitly, given that $\vec x = (x_1,x_2),\vec a=(a_1,a_2),\vec b=(b_1,b_2)$, I want to compute:

$$\langle x_1x_2\rangle = \int_{a_1}^{b_1}\mathrm dx_1\int_{a_2}^{b_2}\mathrm dx_2\ x_1x_2 f(x_1,x_2)$$

$$\langle x_i\rangle = \int_{a_1}^{b_1}\mathrm dx_1\int_{a_2}^{b_2}\mathrm dx_2\ f(x_1,x_2) x_i,\qquad i = 1,2$$

$$\langle x_i^2\rangle = \int_{a_1}^{b_1}\mathrm dx_1\int_{a_2}^{b_2}\mathrm dx_2\ f(x_1,x_2) x_i^2,\qquad i = 1,2$$

I presume there are no analytical formulas, so I am looking for an efficient and robust numerical method, probably involving integration. Specifically, I am worried about the scenario where the truncation rectangle is far from the peak $\vec \mu$. In this case a naive integration might fail due to underflow, even though the moments are well defined even in this case.

Update: By re-scaling and translation, we can assume that $\vec \mu = 0$ and that

$$\Sigma=\left(\begin{array}{cc} 1 & -\rho\\ -\rho & 1 \end{array}\right)$$

for some $-1<\rho<1$. An example set of values where the packages I have tried fail is: $\rho = 0.0220$, $\vec a = (724.128, -0.324)$, $\vec b = (2518.364, 0.511)$. I always get NaNs.

Answer 1

The following is a solution strategy for finding an analytical expression that solves the case of finding the means (at least it solves the problem analytically, the expression is not much simplified and computation, of the analytical expression, might still be difficult, yet simpler because more is know about computing exponential functions and error functions)

On page 152 equation 43 of

Murray Geller and Edward W. Ng, 1971, A Table of Integrals of the Error Function. II Additions and Corrections J. Res. Natl. Bur. Stand., 75, pp.149-163

you will find

$$ 2b^2 \int \text{erf}(ax+c) \text{exp}(-b^2x^2)x \,dx = \\ \frac{a}{\sqrt{a^2+b^2}} \text{erf}\left(x \sqrt{(a^2+b^2)} + \frac{ac}{\sqrt{a^2+b^2}} \right) \text{exp} \left( -\frac{b^2c^2}{a^2+b^2} \right) -\text{erf}(ax+c) \text{exp}(-b^2x^2)$$

which can be used to solve your problem (the case of finding the means) in terms of sums and products of exponentials and error functions. (although it is a bit awkward to write it down fully).

---

The mode of solution would be:

• Transform your problem into a case with no correlation $\Sigma_{12} = \Sigma_{21} = 0$.

  Note that the rectangular truncation is now rotated.

• Change the probability, which is a product of two exponential terms, into a marginal probability by integrating one of the variables. Then one of the exponential terms becomes an expression with terms like $\text{erf}(ax+c)$ (the error function as part of the expression for the CDF of a Gaussian and the linear expression for the rotated truncation)

  You will have to do this for different regions since linear functions for the diagonal bounding box are not the same for the entire region over which the marginal probability is evaluated.

• Work out the resulting integrals based on the indefinite integral.

While this method will give you a result that is reduced into simpler terms of error functions and exponentials, it will look like a monster.