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

I'd like a closed form solution for $$E(X_1 \mid X_2 > X_3)$$ where $$(X_1, X_2, X_3)$$ is multivariate normal with possibly arbitrary mean vector and covariance matrix.

The conditional distribution $$f(X_1 \mid X_2, X_3)$$ 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 $$X_2 > X_3$$. Any help would be appreciated!

• To echo what @MattF. mentioned, there won't be nice closed-form expressions mainly because of the need to integrate a normal CDF. The special cases for where there is a closed-form are likely not very useful. For example if all of the means are $\mu$ and variances are $\sigma^2$ and $X_1$ and $X_2$ are both independent of $X_3$, then the conditional expectation is $\mu+\rho_{12} \sigma/\sqrt{\pi}$. If $X_2$ is independent of both $X_1$ and $X_3$ (with as before all means and variances are equal), the conditional mean is $\mu-\rho_{13} \sigma/\sqrt{\pi}$. There might be other special cases.
– JimB
Commented Feb 14 at 4:55
• Thanks, that's very helpful. I think I could relax "closed form" to "easily evaluated with a computer." For example, if I need to do MCMC for an hour, that is probably too much. But if I all I need to do is evaluate the normal CDF or maybe even a numerical integral of a relatively simple expression that contains a term with the normal CDF, that is ok. Thanks! Commented Feb 14 at 13:18

After thinking a little bit, I was able to solve this. The solution is fairly straightforward. Suppose that $$(X_1, X_2, X_3)$$ 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 $$\operatorname{Cov}(X_1, D) = V_{12} - V_{23}$$.

In other words, $$(X_1, 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]$$

• +1 Very good! (Glad I'm wrong.)
– JimB
Commented Feb 17 at 16:50
• One nice corollary here is that $E(X_1|D>0)>E(X_1)$ iff $X$ and $D$ are positively correlated. Meanwhile the final ratio could also be written with all positive terms as $$\left[\frac{\phi(\frac{\mu_D}{\sqrt{V_d}})}{\Phi(\frac{\mu_D}{\sqrt{V_d})})}\right]$$
– user225256
Commented Feb 26 at 20:36

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 *)


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