Distribution of the a random variable defined on the index of a set of independent random variables Scenario 1: Let $X_i,i\in I$ where $I$ is the index set. For simplicity let $I$ be the finite set $I=\{1,2,\ldots,n\}$. Let each $X_i$ be independent and distributed as
$$
X_i\sim N(\delta,1),
$$
and let $\mathbf{X}=(X_1,X_2,\ldots,X_n)$.
Define the  index estimator 
$$
\gamma(\mathbf{X})=\text{argmax}_{i\in I}{\mathbf{(X)}}
$$
The distribution of $\gamma(\mathbf{X})$ is uniform. This is easy to see because each $X_i$ is equally likely to be the maximum since they are all independent and identically distributed.
Scenario 2: In this scenario, $X_i$ are still independent, however, not identically distributed. Each $X_i$ has the distribution 
$$
X_i\sim N(\delta_i,1),
$$
and $\gamma(\mathbf X)$  has the same definition, i.e., 
$$
\gamma(\mathbf{X})=\text{argmax}_{i\in I}{\mathbf{(X)}}.
$$
My question is what is the distribution of $\gamma(\mathbf{X})$ in this scenario?
I started with 
\begin{eqnarray}
\Pr(\gamma(\mathbf X)=i)&=&\Pr(X_i=\max\{X_1,X_2,\ldots,X_n\})\\
                        &=& \Pr(X_i\geq\max\{X_k:k\neq i,k=1,2,\ldots,n\})\\
\end{eqnarray}
And define $Y=\max\{X_k:k\neq i,k=1,2,\ldots,n\}$ which implies 
$$
\Pr(Y\leq y)=\prod_{k \neq i}\Pr(X_k\leq y)=\prod_{k \neq i}\Phi(y)
$$
Then my problem simplifies to 
$$
\Pr(\gamma(\mathbf X)=i)=\Pr(X_i-Y\geq 0).
$$
It is the distribution of $X_i-Y$ that makes life complicated for me. However, if you can proceed in another way it might be easier. 
 A: For $n = 2$, $p(\gamma(x) = 1) = p(x_{1} - x_{2} > 0)$. Let $z_{1,2} = x_{1} - x_{2}$. Obviously, $z_{1,2} \sim N(\delta_{1} - \delta_{2}, 2)$. Hence, $p(\gamma(x) = 1) = 1 - \Phi((\delta_{1} - \delta_{2}) / \sqrt{2})$.
For $n > 2$, we use multivariate normal distribution. Take $n = 3$ for example. $p(\gamma(x) = 3) = p(x_{3} - x_{2} > 0, x_{3} - x_{1} > 0)$. Let $z_{3,2} = x_{3} - x_{2}$ and $z_{3,1} = x_{3} - x_{1}$. Of course, the vector $(z_{3,2}, z_{3,1})$ follows bivariate normal:
\begin{equation*}
 (z_{3,2}, z_{3,1})' \sim 
 N
 \left(
  (\delta_{3} - \delta_{2}, \delta_{3} - \delta_{1})',
  \begin{pmatrix}
   2 & 1 \\
   1 & 2
  \end{pmatrix}
 \right)
\end{equation*}
Then the probability of $p(\gamma(x) = 3) = p(z_{3,2} > 0, z_{3,1} > 0)$ can be calculated from the CDF of the above bivariate normal distribution.
The general case of $n > 2$ is similar. You will define a $n-1$-dimensional multivariate normal, whose covariance matrix has diagonal elements of $2$ and off-diagonal elements of $1$.
