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

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.

• – becko Feb 14 '18 at 13:10
• Exactly what form does this truncation take? For bivariate variables, truncation can be arbitrarily complicated, because generally it refers to restricting the distribution to some measurable subset of the plane. Unless you make severe restrictions on what that subset might be, your only hope is numerical integration; and what form to use depends on the subset. But why worry about truncation far from the mean? That will have almost no effect at all on the low moments of any Normal distribution. – whuber Feb 14 '18 at 14:10
• @whuber Rectangular truncation. I edited to clarify. – becko Feb 14 '18 at 15:04
• What do you mean by the last statement: But why worry about truncation far from the mean? That will have almost no effect at all on the low moments of any Normal distribution.? Far from the mean the exponential gets very small, which can lead to numerical problems I think. At least I had some issues in the univariate case (see related link in my previous comment). – becko Feb 14 '18 at 15:05
• Right: there is almost no probability far from the mean, so removing it will scarcely affect the lower moments. BTW, since your inequalities are componentwise, there doesn't seem to anything bivariate about your question: just compute the means and variances of the truncated marginals (which doesn't require numerical integration). Or do you also want the covariance? – whuber Feb 14 '18 at 15:23

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.

• I do at the moment notice a problem that this method does not (directly)provide the constant of proportionality for the relative expression $f(x) \propto exp(...)$. The used integral above is like $$\int erf(ax+c) exp(-b^2x^2) x dx$$ and you will also need $$\int erf(ax+c) exp(-b^2x^2) dx$$ (without the $x$). – Sextus Empiricus May 28 '18 at 23:41