# What is the probability that random variable $x_1$ is maximum of random vector $X=(x_i)$ from a multivariate normal distribution?

Given a $n$-dimensional multivariate normal distribution $X=(x_i) \sim \mathcal{N}(\mu, \Sigma)$ with mean $\mu$ and covariance matrix $\Sigma$, what is the probability that $\forall j\in {1,\ldots,n}:x_1 \geq x_j$?

The question reads to me like the OP was asking when $U = (X,Y,Z)^{\mathrm{T}}$ are jointly normal then what is the probability $P(X \geq Y \mbox{ and } X \geq Z)$?

For that question we could look at the joint distribution of $AU$ where $A$ looks like
$$A=\left[ \begin{array}{ccc} 1 & -1 & 0 \newline 1 & 0 & -1 \end{array}\right]$$
Of course, $AU$ is also jointly normal with mean $A\mu$ and variance-covariance $A\Sigma A^{\mathrm{T}}$, and the desired probability is $P(AU > \mathbf{0}_{n-1})$. We could get this in R with something like

set.seed(1)

Mu <- c(1,2,3)

library(MCMCpack)
S <- rwish(3, diag(3))  # get var-cov matrix

A <- matrix(c(1,-1,0, 1,0,-1), nrow = 2, byrow = TRUE)

newMu <- as.vector(A %*% Mu)
newS <- A %*% S %*% t(A)

library(mvtnorm)
pmvnorm(lower=c(0,0), mean = newMu, sigma = newS)


which is about 0.1446487 on my system. If a person knew something about the matrix $\Sigma$ then (s)he might even be able to write something down that looks like a formula (I haven’t tried, though).