What is the distribution of a bivariate normal component conditional on the max of the other component? Let $n$ be a large integer, and consider two independent multivariate Gaussian $n$-vectors $x, z$ with $x\sim\mathcal{N}\left(0,I\right),$ and $z\sim\mathcal{N}\left(0,\sigma^2 I\right)$. Let $y=x+z$. Conditional on observing $y$, consider the reindexing $(i)$ such that $y_{(1)} \le y_{(2)} \le \cdots \le y_{(n-1)} \le y_{(n)}$. What can we say about $x_{(n)}$?


*

*Do we know the conditional distribution of $x_{(n)}$?

*Failing that, can we compute or estimate the conditional expected value of $x_{(n)}$?

*In each of these, what are the asymptotics of $x_{(n)}$ as $n\to \infty$?


background: this is like an estimation after selection procedure where the $x$ are some unobserved population parameter that you observe with noise $z$, and the indices are over different 'models'. You select the model with the largest estimated parameter.
 A: $Y_{(n)}$ is the largest of $n$ points drawn from a bivariate Normal distribution.  Its $X$ coordinate is more or less likely to be an extreme value among the corresponding $X$ coordinates provided $X$ and $Y$ are strongly correlated.  Even then, the variability in $X$ limits the degree to which this can happen: a manifestation of the regression to the mean phenomenon.

This scatterplot of a sample of size $n=128$ highlights the point with largest $Y$ value, $Y_{(n)}.$ The colored lines show the coordinates.  Note that its x-coordinate, corresponding to the position of the vertical blue line, is not the largest among the x coordinates.  Regression theory explains that this x coordinate $X_{(n)}$ will vary in a Normal fashion from the value predicted from $Y_{(n)}.$  (Note the reversal in the usual roles of $X$ and $Y:$ the setup of the question requires us to predict $X$ from $Y$ rather than the other way around.)

The conditions given in the question imply $(X,Y)$ has a bivariate Normal distribution.  Its correlation coefficient is
$$\rho = \frac{\operatorname{Cov}(X,Y)}{\sqrt{\operatorname{Var}(X)\operatorname{Var}{Y}}} = \frac{\operatorname{Cov}(X,X+Z)}{\sqrt{\operatorname{Var}(X)\operatorname{Var}(X+
Z)}} = \frac{1}{\sqrt{(1)(1+\sigma^2}} = \sqrt{\frac{1}{1+\sigma^2}}.$$
We may analyze all bivariate Normal distributions at once, recognizing that suitable choices of units for $X$ and $Y$ reduce it to the condition where $X$ and $Y$ are standard Normal with correlation $\rho.$  (In particular, $X$ is already standard Normal and so needs no standardization at all.)
Let $\Phi$ be the standard Normal CDF.  Following the analyses at Extreme Value Theory - Show: Normal to Gumbel, the CDF of the maximum of $n$ iid observations $Y_{(n)}$ is $\Phi^n$ which, for large $n$, can be standardized to approximate a Gumbel distribution.  At the same time, the theory of linear regression establishes that the distribution of $X$ is that of $\rho Y + \sqrt{1-\rho^2}\,\varepsilon$ where $\varepsilon$ has a standard Normal distribution (independent of $(X,Y)$).  Therefore 

$X$ is distributed as the linear combination of independent variables with $\Phi^n$ and $\Phi$ distribution functions.  The coefficients are $\rho$ and $\sqrt{1-\rho^2},$ respectively.

Asymptotically, an appropriately scaled version of $Y$ has a Gumbel distribution (id.).  The scaling factor decreases to zero in the limit.  One asymptotic expression for the location factor is
$$\mu(n) = \Phi^{-1}\left(1-\frac{1}{n}\right).$$
Therefore

Asymptotically, $(X - \rho\mu(n))/\sqrt{1-\rho^2}$ has a standard Normal distribution.


The following figures present simulation results that check and illustrate these conclusions.  "$n$" is the sample size.  A correlation coefficient of $\rho=0.8$ was used.  Each of the four (independent) simulations obtained $100,000$ independent realizations of $(X_{(n)}, Y_{(n)}).$  The upper plots are histograms of the realized values of $Y_{(n)}$ and the lower plots are the corresponding histograms of the realized values of $X_{(n)}.$  The overplotted colored curves are the asymptotic Gumbel (red, top) and Normal (blue, bottom) distributions, scaled to fit the data.

It is evident that (a) the Gumbel approximation is good by the time $n=100$ and (b) although the convergence of $X$ to a normal variable is slow, by the time $n=10^{25}$ the difference is undetectable.

The R code that produced these results runs almost instantaneously--demonstrating the relative simplicity of all the calculations--and is provided for further experimentation.
#
# Generate random variates far into the upper tail of the standard Normal
# distribution (Mills' Ratio approximation).
#
qnorm.0 <- function(log.q) {
  f <- function(x) sqrt(-2*log(-sqrt(2*pi) * x * log.q))
  ifelse(log.q < -1e-12, qnorm(exp(log.q)), f(f(f(7))))
  # obj <- uniroot(function(x) sqrt(2*pi) * x * log.q + exp(-x^2/2), c(lower, upper))
  # obj$root
}
#
# Gumbel PDF.
#
dGumbel <- function(x, mu=0, sigma=1) {z <- (x-mu)/sigma; exp(-z - exp(-z))/sigma}
#
# Standardization of Normal maximum distribution for largish n.
#
snorm <- function(n) {
  c(mu = qnorm.0(-1/n), sigma = 1 / sqrt(2*log(n) - log(2*pi))) 
}
#
# Run the simulations, plotting as we go.
#
rho <- 0.8
N <- 1e5          # Simulation size
set.seed(17)      # For reproducibility
par(mfcol=c(2,4))
for (n in 10^c(1, 2, 5, 25)) {
  # Generate and display Y
  y <- qnorm.0(-rexp(N, n))
  hist(y, main=paste("n =", n), breaks=30, freq=FALSE, 
       xlab=expression(y[(n)]), cex.lab=1.25)
  ab <- snorm(n)
  curve(dGumbel(x, ab[1], ab[2]), add=TRUE, col="#e02020", lwd=2)

  # Generate and display X
  x <- rho * y + sqrt(1-rho^2) * rnorm(N, 0, 1)
  tau <- sqrt(1-rho^2)
  hist(x, breaks=30, freq=FALSE, ylim=c(0, sqrt(1/(2*pi)) / tau), 
       xlab=expression(x[(n)]), cex.lab=1.25)
  curve(dnorm(x, rho*ab[1], tau), add=TRUE, col="#2020e0", lwd=2)
}
par(mfrow=c(1,1))
# #
# # Check the fast method (above) by comparison to the straightforward one.
# # The QQ plot lies close to the X=X.1 line, demonstrating equivalence.
# #
# library(MASS)
# n <- 1000
# N <- 1e4
# xy <- aperm(array(mvrnorm(n*N, c(0,0), matrix(c(1,rho,rho,1), 2)), c(N,n,2)), c(2,3,1))
# 
# x <- apply(xy, 3, function(xy) {
#   i <- which.max(xy[, 2])
#   xy[i, 1]
# })
# # hist(x, main=paste0("n =", n))
# 
# x.1 <- rho * qnorm(runif(N)^(1/n)) + sqrt(1-rho^2) * rnorm(N, 0, sigma)
# plot(sort(x), sort(x.1), main="q-q plot", xlab="Direct", ylab="Shortcut")
# abline(c(0,1), col="Red", lwd=2)

A: My comment to @whuber is based on the following reasoning:
Suppose that $x \sim \mathcal{N}\left(\mu_x,\sigma_x^2\right)$,
and $y \left| x \right. \sim \mathcal{N}\left(x,\sigma_y^2\right)$.
Then the stacked vector of $x, y$ is bivariate normal:
$$
\left[x, y\right]^{\top} \sim \mathcal{N}\left(\left[\mu_x, \mu_x\right]^{\top}, 
\left[\begin{array}{cc}
\sigma_x^2 & \sigma_x^2\\
\sigma_x^2 & \sigma_x^2 + \sigma_y^2
\end{array}\right]
\right).
$$
By the conditional normal distribution, 
conditional on $y=a$ for some $a$, $x$ is normal
$$
x \left| \left\{ y = a \right\} \right. \sim \mathcal{N}\left(\mu_x +
\frac{\sigma_x}{\sqrt{\sigma_x^2 + \sigma_y^2}} \rho \left(a - \mu_x\right),
\left(1-\rho^2\right)\sigma_x^2\right),
$$
where $\rho$ is the correlation of $x$ and $y$, which takes value
$$
\rho = \frac{\sigma_x^2}{\sigma_x \sqrt{\sigma_x^2 + \sigma_y^2}} 
= \frac{\sigma_x}{\sqrt{\sigma_x^2 + \sigma_y^2}}.
$$
Thus
$$
x \left| \left\{ y = a \right\} \right. \sim
\mathcal{N}\left(\left(1-\rho^2\right) \mu_x + \rho^2 a,
\left(1-\rho^2\right)\sigma_x^2\right).
$$
Substituting in $\mu_x=0$ and $\sigma_x=1$ gives the desired relation.
