# Computing the covariance of the truncated normal distribution

In the paper On Moments of Folded and Truncated Multivariate Normal Distributions on page 17, one can find the explicit expression for low order moments of the truncated multivariate normal distribution with $$X \mid X > \alpha$$. However, I do not get the right results when comparing it to some code in Python and R.

Let $$\mu = \begin{bmatrix}0 & 0\end{bmatrix}^T$$ be the mean, $$\Sigma = \begin{bmatrix}1 & 0.1\\0.1 & 1\end{bmatrix}$$ the covariance and $$\alpha = \begin{bmatrix}0 & 0\end{bmatrix}^T$$ the truncation points.

The mean of the univariate truncated normal is here

$$E[X_i] = \mu + \sigma\frac{\phi(\frac{-\mu}{\sigma})}{\Phi(\frac{\mu}{\sigma})} = \frac{\phi(0)}{\Phi(0)}$$

Using the formula from the paper while setting $$\mu_i = 0$$ and $$\alpha_i = 0$$ yields

$$E[X_iX_j] = \rho_{ij} + \frac{(1 - \rho_{ij}^2)\phi_2(0, 0 ; 0, \rho_{ij})}{\Phi_2(0, 0 ; 0, \rho_{ij})}$$

where $$\rho_{ij} = \frac{0.1}{1}$$ is the correlation coefficient, $$\phi_2$$ is the bivariate PDF and $$\Phi_2$$ the bivariate CDF.

Then the covariance of $$X_1$$ and $$X_2$$ is

\begin{align} \text{Cov}(X_1, X_2) &= E[X_1X_2] - E[X_1]E[X_2]\\ &= E[X_1X_2] - \left(\frac{\phi(0)}{\Phi(0)}\right)^2 \end{align}

from scipy.stats import norm
from scipy.stats import multivariate_normal as mvn

E_X_i = (norm.pdf(0)/norm.cdf(0))
print(E_X_i) # result = 0.7978845608028654
E_X_i_X_j = 0.1 + (1 - 0.1 ** 2) * mvn.pdf(0, mean=0, cov=0.1)/mvn.cdf(0, mean=0, cov=0.1)
cov = E_X_i_X_j - E_X_i ** 2
print(cov) # result = 1.9612814244323769


This result is wrong because in R the library tmvtnorm outputs:

install.packages("tmvtnorm")
library(tmvtnorm)

mtmvnorm(mean=c(0, 0), sigma=matrix(c(1, 0.1, 0.1, 1), 2, 2), lower=c(0, 0), upper=c(Inf, Inf))
$tmean  0.8250601 0.8250601$tvar
[,1]       [,2]
[1,] 0.37882152 0.01473302
[2,] 0.01473302 0.37882152


The right answer is $$\text{Cov}(X_i, X_j) = 0.01473302$$. Also $$E[X_i] = 0.8250601$$ is slightly off (but numerical stability could be the issue).

Furthermore, I tried integrating the truncated distribution here, but with different results.

Where lies the problem?

It took me some time, but I found the solution.

1. tmvtnorm is actually numerical stable. I just did not consider that $$\mu_i$$ also changes in the bivariate case. For $$\mu = \begin{bmatrix}0 & 0\end{bmatrix}^T$$, $$\Sigma_{ii} = 1$$ and $$X_i \mid X_i > 0$$, the new $$\mu$$ is:

$$\mu_i = \frac{\rho\phi(0)\Phi(0) + \phi(0)\Phi(0)}{\Phi_2(\begin{bmatrix}0 & 0\end{bmatrix}^T; \begin{bmatrix}0 & 0\end{bmatrix}^T, \Sigma)}$$

1. I did not enter the parameters correctly. It should be:
mvn.cdf(np.array([0, 0]), mean=np.array([0, 0]), cov=np.array([[1, 0.1], [0.1, 1]]))

1. The formula in the paper contains typos. For the case $$\mu = \begin{bmatrix}0 & 0\end{bmatrix}^T$$, $$\Sigma_{ii} = 1$$ and $$X_i \mid X_i > 0$$, the formula should be:

$$E[X_iX_j] = \rho_{ij} + \frac{\sqrt{1 - \rho_{ij}^2}\phi^2(0)}{\Phi_2(\begin{bmatrix}0 & 0\end{bmatrix}^T; \begin{bmatrix}0 & 0\end{bmatrix}^T, \Sigma)}$$

Then the corrected code is:

from scipy.stats import norm
from scipy.stats import multivariate_normal as mvn
import numpy as np

E_X_i = (0.1 * norm.pdf(0) * norm.cdf(0) + norm.pdf(0) * norm.cdf(0))/mvn.cdf(np.array([0, 0]), mean=np.array([0, 0]), cov=np.array([[1, 0.1], [0.1, 1]]))
print(E_X_i) # result = 0.8250601203844766
E_X_i_X_j = 0.1 + norm.pdf(0) ** 2 * np.sqrt(1 - 0.1 ** 2) / mvn.cdf(np.array([0, 0]), mean=np.array([0, 0]), cov=np.array([[1, 0.1], [0.1, 1]]))
cov = E_X_i_X_j - E_X_i ** 2
print(cov) # result = 0.014733023153836067


This result is like tmvtnorm 0.01473302. Hence, my result is correct.