For $X_i \sim \text{N}(\mu_i, \sigma_i)$, is there a fast method of (roughly) computing
$$\text{P}\left(X_i = \text{max}(X_1, \dots, X_n)\right), \forall i$$
for a set of $n$ tuples of $(\mu_i, \sigma_i)$, $i \in 1, \dots, n$? Or in words, for all $i$, what is the probability that $X_i$ is the maximum in the sample of $\{X_1, \dots, X_N\}$
I believe the exact answer is given by
$$ \text{P}\left(X_i = \text{max}(X_i, \dots, X_n)\right) = \frac{1}{\sqrt{2\pi\sigma_i}} \int_{-\infty}^{\infty} \exp\left(\frac{-(x-\mu_i)^2}{2\sigma_i}\right) \prod_{k \neq i} \Phi\left(\frac{x-\mu_k}{\sigma_k}\right) \text{dx} $$
where $\Phi(x)$ is the cumulative normal distribution function.
I need to repeat this calculation for $M$ different sets of $n$ tuples. I am wondering if there is perhaps a formula which could give an approximation that could possibly be vectorised for fast computation when $n$ is known in advance. Ideally, avoiding having to perform $M\times n$ numerical integrations.
One thought I had was to train a model.
Given that I can calculate the exact values in advance for a set of $n$ tuples, perhaps I could train a machine learning model which might be able to give me a fast but adequate approximation when $n$ is fixed? I could set it up as a multiple response model where the input data consists of rows of vectors made from concatenating $n$ randomly generated $(\mu_i, \sigma_i)$ tuples and the responses are calculated using the formula above. I can generate as much data as I need.
We have
\begin{align} P\left(X_i = \max(X_1, \dots, X_n)\right) &= E\left(\mathbb{I}_{\{X_i = \max(X_1, \dots, X_n)\}} \right)\\ &= E\left(E(\mathbb{I}_{\{X_i = \max(X_1, \dots, X_n)\}}|X_i) \right)\\ &= E\left(P(X_i = \max(X_1, \dots, X_n)|X_i) \right)\\ &= E\left(P(\bigcap_{k \ne i} \{X_k \le X_i\} |X_i) \right)\\ &= E\left( \prod_{k \ne i} P(X_k \le X_i |X_i) \right) \tag{1}\\ &= E\left( \prod_{k \ne i} \Phi\left(\frac{X_i-\mu_k}{\sigma_k}\right) \right) \\ &= \frac{1}{\sqrt{2\pi\sigma_i}}\int_\mathbb{R}\left(\exp\left(\frac{-(x-\mu_i)^2}{2\sigma_i}\right) \prod_{k \neq i} \Phi\left(\frac{x-\mu_k}{\sigma_k}\right)\right)dx \end{align}
At $(1)$, we need the assumption that $(X_i)_{i=1,...,n}$ are independant.