1
$\begingroup$

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
[1] 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?

$\endgroup$

1 Answer 1

1
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.