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.
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). $\endgroup$