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!
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]$$
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
