E(X1 | X2 > X3) for (X1,X2,X3) multivariate normal

Question

I'd like a closed form solution for $E(X1 | X2 > X3)$ where $(X1, X2, X3)$ is multivariate normal with possibly arbitrary mean vector and covariance matrix.

The conditional distribution $f(X1 | X2, X3)$ is easy to get in closed form, but I am not sure if there is a closed-form solution to handling conditioning on the event $X2 > X3$. Any help would be appreciated!

Answer 1

After thinking a little bit, I was able to solve this. The solution is fairly straightforward. Suppose that $(X1, X2, X3)$ is multivariate normal with the following mean and covariance:

$$ \mu = \begin{bmatrix} \mu_1 \\ \mu_2 \\ \mu_3 \end{bmatrix} \qquad \Sigma = \begin{bmatrix} V_{11} & V_{12} & V_{13} \\ V_{21} & V_{22} & V_{23} \\ V_{31} & V_{32} & V_{33} \end{bmatrix}. $$

We can define $D \equiv X_2 - X_3$. Trivially, $E(D) = \mu_2 - \mu_3$, $V(D) = V_{22} + V_{33} - 2V_{23}$, and $Cov(X_1, D) = V_{12} - V_{23}$.

In other words, $(X1, D)$ is multivariate normal with the following mean and covariance:

$$ \mu^* = \begin{bmatrix} \mu_1 \\ \mu_2 - \mu_3 \end{bmatrix} \qquad \Sigma^* = \begin{bmatrix} V_{11} & V_{12} - V_{23} \\ V_{12} - V_{23} & V_{22} + V_{23} - 2V_{23} \end{bmatrix}. $$

For convenience, let's just rename these parameters:

$$ \mu^* = \begin{bmatrix} \mu_1 \\ \mu_D \end{bmatrix} \qquad \Sigma^* = \begin{bmatrix} V_{11} & V_{1,D} \\ V_{1,D} & V_{D} \end{bmatrix}. $$

Now, we want $E(X_1 \mid D > 0)$, which can be solved using the answer to this question. The solution is:

$$E(X_1 \mid D > 0) = \mu_1 + \frac{V_{1,D}}{\sqrt{V_{D}}} \left[\frac{\phi\left(\frac{-\mu_D}{\sqrt{V_d}}\right)}{1 - \Phi\left(\frac{-\mu_D}{\sqrt{V_d}}\right)}\right]$$

Answer 2

Either random sampling or numerical integration should work fine for this. Here is an implementation using Mathematica:

(* Set parameters *)
Σ = {{1, 0.5, 0.7}, {0.5, 2, 0.7}, {0.7, 0.7, 1}};
μ = {3, 2, 1};
(* Random sampling )
SeedRandom[12345];
x1x2x3 = RandomVariate[MultinormalDistribution[μ, Σ], 1000000];
x1Givenx2x3 = Select[x1x2x3, #[[2]] > #[[3]] &][[All, 1]];
Mean[x1Givenx2x3]
( 2.939435703127434 )
se = StandardDeviation[x1Givenx2x3]/Sqrt[Length[x1Givenx2x3]]
( 0.0011216623753498563 *)
(* Numerical integration )
m = (NIntegrate[x1  PDF[MultinormalDistribution[μ, Σ], {x1, x2, x3}],
      {x1, -∞, ∞}, {x2, -∞, ∞}, {x3, -∞, x2}] // Quiet)/
   NIntegrate[PDF[MultinormalDistribution[μ, Σ], {x1, x2, x3}],
    {x1, -∞, ∞}, {x2, -∞, ∞}, {x3, -∞, x2}] // Quiet
( 2.9412419726547325 *)

Answer 3

Let the mean vector and covariance matrix of $(X_1, X_2, X_3)$ be denoted as $\mu$ and $\Sigma$, respectively, where:

$$\mu = \left[ \begin{array}{c} \mu_1 \\ \mu_2 \\ \mu_3 \end{array} \right]$$

and

$$\Sigma = \left[ \begin{array}{ccc} \sigma_{11} & \sigma_{12} & \sigma_{13} \\ \sigma_{21} & \sigma_{22} & \sigma_{23} \\ \sigma_{31} & \sigma_{32} & \sigma_{33} \end{array} \right].$$

The conditional distribution of $X_1$ given $X_2$ and $X_3$ is also normally distributed, with its mean and variance given by the properties of the multivariate normal distribution. However, the condition $X_2 > X_3$ complicates the calculation, turning it into a problem involving conditional distributions based on an inequality.

For multivariate normals, the conditional expectation $E(X_1 | X_2, X_3)$ can be directly computed if $X_2$ and $X_3$ were given as constants. However, given the condition $X_2 > X_3$, we have to integrate over the joint distribution of $X_2$ and $X_3$ where $X_2 > X_3$, weighted by the conditional density of $X_1$ given $X_2$ and $X_3$.

The closed form solution for $E(X_1 | X_2 > X_3)$ is not straightforward because it depends on the joint distribution of $X_2$ and $X_3$, and involves integrating over a region defined by $X_2 > X_3$. In general, the solution might involve numerical methods or approximations because the integral defining $E(X_1 | X_2 > X_3)$ does not usually simplify nicely for arbitrary covariance matrices.

But in the most general case, where $\Sigma$ allows for arbitrary correlations among $X_1$, $X_2$, and $X_3$, finding a closed-form solution for $E(X_1 | X_2 > X_3)$ requires integrating over the conditional density of $X_1$ given $X_2$ and $X_3$, subject to $X_2 > X_3$, which is typically done using numerical methods rather than analytical expressions.

To approximate $E(X_1 | X_2 > X_3)$ for a multivariate normal distribution with mean vector $\mu = \{3, 2, 1\}$ and covariance matrix $$\Sigma = \left[ \begin{array}{ccc} 1 & 0.5 & 0.7 \\ 0.5 & 2 & 0.7 \\ 0.7 & 0.7 & 1 \end{array} \right]$$

we can use Monte Carlo simulation. Monte Carlo simulation involves generating samples from the multivariate normal distribution and filtering for the condition $X_2 > X_3$. This method is straightforward and relies on the law of large numbers to approximate the expected value.

library(MASS) # for mvrnorm to generate multivariate normal distributions
Define the mean vector and covariance matrix
mu <- c(3, 2, 1) # The mean vector μ
Sigma <- matrix(c(1, 0.5, 0.7,    # First row of the covariance matrix
                  0.5, 2, 0.7,    # Second row
                  0.7, 0.7, 1),   # Third row
                nrow = 3, ncol = 3, byrow = TRUE) # Reshape into 3x3 matrix
Generate samples from the multivariate normal distribution
set.seed(123) # For reproducibility
samples <- mvrnorm(n = 10000, mu = mu, Sigma = Sigma)
Filter samples where X2 > X3
filtered_samples <- samples[samples[,2] > samples[,3],]
Compute the approximate expectation of X1 given X2 > X3
E_X1_given_X2_gt_X3 <- mean(filtered_samples[,1])
print(E_X1_given_X2_gt_X3
2.93153