Computing the probability $P(\exists X\in\{X_1,\ldots, X_N\}:X>\text{max}\{Y_1,\ldots, Y_M\})$ I am trying to somewhat generalize a question, which has been asked in one way or another a several times here on StackExchange**. However, I have not managed to find an answer to the below problem.
Suppose, $X_1,\ldots,X_N\sim\mathcal{N}(\mu_X, \sigma_X)$ and $Y_1,\ldots,Y_M\sim\mathcal{N}(\mu_Y, \sigma_Y)$, which are independently drawn. We are interested in the following probability:
$$
P(\text{there exists at least one }X_i\text{ such that } X_i \text{ greater than max}\{Y_1,\ldots, Y_M\})=P(\exists X\in\{X_1,\ldots, X_N\}:X>\text{max}\{Y_1,\ldots, Y_M\})\tag{1}
$$
I have tried to decompose this problem into smaller problems by writing the "there exists" part as
$$
P(\exists X\in\{X_1,\ldots, X_N\}:X>\text{max}\{Y_1,\ldots, Y_M\})=P(X_1>\text{max}\{Y_1,\ldots, Y_M\}\cup\ldots\cup X_N>\text{max}\{Y_1,\ldots, Y_M\}),
$$
but the RHS does not contain any information about the partial order within the set of $\{X_i\}$. Trying to add this information to the RHS makes it overly complicated and I am not even sure if this would work out. Is there a known closed form expression or a numerical way to calculate the value in Eq.(1)?
**  $P(X_1 < \min(X_i,\ldots, X_n))$ across different normal random variables
**  What is $P(X_1>X_2 , X_1>X_3,... , X_1>X_n)$?
** What is a method to calculate precisely $P(Y \geq X, Y\leq Z)$, given three independent random variables $X, Y$, and $Z$
 A: First of all, you can simplify this condition to $X_{(N)} > Y_{(M)}$ using the notation for order statistics, which will simpify your question a fair bit.  Assuming the values are all independent, these two random variables will also be independent, with respective distribution functions:$^\dagger$
$$F_{X_{(N)}}(x) = \Phi \Big( \frac{x-\mu_X}{\sigma_X} \Big)^N
\quad \quad \quad \quad \quad 
F_{Y_{(M)}}(y) = \Phi \Big( \frac{y-\mu_Y}{\sigma_Y} \Big)^M,$$
and corresponding densities:$^\dagger$
$$\begin{align}
f_{X_{(N)}}(x) = \frac{N}{\sigma_X} \cdot \phi \Big( \frac{x-\mu_X}{\sigma_X} \Big) \cdot \Phi \Big( \frac{x-\mu_X}{\sigma_X} \Big)^{N-1} \\[18pt]
f_{Y_{(M)}}(y) = \frac{M}{\sigma_Y} \cdot \phi \Big( \frac{y-\mu_Y}{\sigma_Y} \Big) \cdot \Phi \Big( \frac{y-\mu_Y}{\sigma_Y} \Big)^{M-1}. \\[6pt]
\end{align}$$
Thus, we can write the probability of interest as:
$$\begin{align}
\mathbb{P}(X_{(N)} > Y_{(M)}) 
&= \int \limits_\mathbb{R} f_{X_{(N)}}(y) \cdot \mathbb{P}(Y_{(M)} \leqslant y) \ dy \\[6pt]
&= \int \limits_\mathbb{R} f_{X_{(N)}}(y) \cdot F_{Y_{(M)}}(y) \ dy \\[6pt]
&=  \frac{N}{\sigma_X} \cdot \int \limits_\mathbb{R} \phi \Big( \frac{y-\mu_X}{\sigma_X} \Big) \cdot \Phi \Big( \frac{y-\mu_X}{\sigma_X} \Big)^{N-1} \cdot \Phi \Big( \frac{y-\mu_Y}{\sigma_Y} \Big)^M \ dy. \\[6pt]
\end{align}$$
There is no closed form expression for this integral, but it can be evaluated numerically to obtain the probability of interest to you.

$^\dagger$ Here we use $\Phi$ and $\phi$ to denote the cumulative distribution function and density function for the standard normal distribution.
