Imagine we have $k$ Gaussian RVs $$ X_i \sim N(\mu_i, \sigma_i^2) \text{ for } i=1, \ldots, k $$ and we sample from each of them independently to produce a vector, $\vec{x} = (x_1, \ldots, x_k)$.
For one of the Gaussian RVs, say $X_j$, I am interested in computing the probability that it exceeds all others, i.e. $$ \Pr\left\{ \cap_{i\not= j} \, X_j > X_i \right\}. $$
I know I can use Monte Carlo sampling to estimate this probability. But are there any closed-form analytical methods or approximations?
The distribution of the maximum of several normal random variables is approximately known
One method you could use to approximate this would be to use some results from extreme value theory. As is mentioned in this answer to Is there an analytical expression for the distribution of the max of a normal k sample? the distribution of the maximum of several normal random variables is approximated by a Gumbel distribution. Similar answers can also be found for Asymptotic distribution of maximum order statistic of IID random normals. While I think these may be tangible for uncorrelated random normals, you get several awkward issues when you introduce correlation which I am not sure can be trivially overcome (Distribution of the maximum of two correlated normal variables).
This reduces the problem down to a single integral
Given you have an expression for the distribution of the maximum, it is then relatively easy to write the probability as a two dimensional integral with the joint distribution (using appropriate integration limits corresponding to $X_j > \max(\{X_i\}_{i\neq j})$). Given the normal should be independent to the Gumbel the joint distribution should be separable and this should collapse nicely to just a single integral. This integral is likely quite ugly, but given the functions are all fairly nice it might be feasible to find where the integrand is maximal and just approximate this region to approximate the probability. While this doesn't solve the problem, it should hopefully reduce it to a simpler problem.
Bounding the error is your biggest challenge if you wish to compare it to Monte Carlo
The biggest issue will be bounding the error from using the Gumbel distribution, and then on top of that having to likely approximate the resultant integral.