Probability that the same r.v. generates the rth order statistic in one noise-added set, and the sth order statistic in another noise-added set

Question

(Note: The title is confusing, as I have no idea if a name / short description exists for the setting below. I'm open to pointers and/or suggestions.)

Setting

Let $X_1, ..., X_N \overset{i.i.d.}{\sim} \mathcal{N}(0, \sigma^2_X)$ for large but fixed $N > 50$. We generate two new sets $Y_1, ..., Y_N$ and $Z_1, ..., Z_N$ by adding two different levels of random noise to $X_i$:

$$\begin{align} Y_i & = X_i + \epsilon_{Yi} , \quad\epsilon_{Yi} \overset{i.i.d.}{\sim} \mathcal{N}(0, \sigma^2_Y)\\ Z_i & = X_i + \epsilon_{Zi} , \quad\epsilon_{Zi} \overset{i.i.d.}{\sim} \mathcal{N}(0, \sigma^2_Z), \end{align}$$ where $X_i$, $\epsilon_{Yi}$ and $\epsilon_{Zi}$ are mutually independent $\forall i$.

We then rank $Y_i$ to get the order statistics $Y_{(1)}, ..., Y_{(N)}$ (where $Y_{(1)}$ represents the smallest of $Y_i$), and rank $Z_i$ to get a different set of order statistics $Z_{(1)}, ..., Z_{(N)}$. The order statistics here can correspond to, say, the result of ranking a set of items under two different noisy estimate of the items' score.

Note there is a 1-1 relationship between the sets $\{X_1, ..., X_N\}$ and $\{Y_{(1)}, ..., Y_{(N)}\}$, and another 1-1 relationship between the sets $\{X_1, ..., X_N\}$ and $\{Z_{(1)}, ..., Z_{(N)}\}$. As a result, two order statistics $Y_{(r)}$ (the $r^{\textrm{th}}$ ranked $Y_n$) and $Z_{(s)}$ (the $s^{\textrm{th}}$ ranked $Z_n$) might actually be generated by the same $X_i$. The question is, how likely is that the case?

Question(s)

The question in its general form reads:

What is the probability, in the setting above, that $Y_{(r)}$ and $Z_{(s)}$ are generated by the same $X_i$, for any given $r, s \leq N$?

I am looking for an explicit formula that may or may not take $N$, $r$, $s$, $\sigma^2_X$, $\sigma^2_Y$, and $\sigma^2_Z$ into account.

Other questions that are perhaps too tiny and coupled to warrant their own CrossValidated question:

1. Is there a name for the "add noise - rank separately - match back" scenario described above?
   I know once you get the ranks $Y_{(r)}$/$Z_{(s)}$, the corresponding $X_{Y[r]}$/$X_{Z[r]}$ is known as the concomitant of the order statistics. However, the concomitants in the literature usually go the other way round, i.e. $X_i$ is ranked, and the behaviour of the induced ranks for $Y_i$/$Z_i$ is studied.
2. Are there anything in the literature that deals with the scenario described?
3. What if I replace the distributions for $X_i$, $\epsilon_{Yi}$, and $\epsilon_{Zi}$ with general distributions $F(\cdot)$, $G(\cdot)$, $H(\cdot)$ respectively, with mean and variance $(\mu_X, \sigma^2_X)$, $(\mu_Y, \sigma^2_Y)$, and $(\mu_Z, \sigma^2_Z)$?

---

Attempts

A natural first answer to the main question is to set the probability as $\frac{1}{N}$. This assumes that $Y_{(r)}$ is independent to $Z_{(s)}$ $\forall r, s$, and $Z_{(s)}$ can point to any $X_i$, one (out of $N$) of which happens to have also generated $Y_{(r)}$. This attempt is clearly wrong as $Y_{(r)}$ and $Z_{(s)}$ are not independent --- see e.g. David and Nagaraja (2004) which discuss the covariance between these quantities.

Some intuitive cases why the above is wrong: $Y_{(N)}$ and $Z_{(N)}$ are not always, but still quite likely to generated by the same $X_i$: If $Y_i$ is the largest amongst its set, chances are that the corresponding $X_i$ is quite large as well, and thus the corresponding $Z_i$ stands a higher-than-random chance to become the largest amongst its set.

On the other hand, one will have a hard time convincing others that $Y_{(1)}$ and $Z_{(N)}$ are generated by the same $X_i$ --- you will need $X_i$ to be somewhere in the middle, $\epsilon_{Yi}$ to be very negative, and $\epsilon_{Zi}$ to be very positive for that to happen. The latter two rarely happen together.

Initial simulations

I've also run some simulation to get a feel of how the probabilities behave. This require reformulating the main problem a bit as:

Fix $r < N$, what is the distribution across all $s \in \{1, ..., N\}$, where $X_i$ generated $Y_{(r)}$ and $Z_i = X_i + \epsilon_{Zi}$ is the $s^{\textrm{th}}$ ranked item?

We can get the answer to the original question above by picking out the probability corresponding to a particular $s$.

The following is a Python snippet that generates an empirical probability mass function for the reformulated problem. I've taken one off the ranks to enable the comparison with beta binomial distributions (see below).

import numpy as np
from scipy.stats import rankdata

N = 100; r = 1; n_runs = 100000
s = []

for run in range(0, n_runs):
    Xi = np.random.normal(0, 1, N)
    Yi = Xi + np.random.normal(0, 0.5, N)
    Zi = Xi + np.random.normal(0, 0.4, N)

    Ir = np.argwhere(rankdata(Yi) == r).flatten()[0]
    s.append(rankdata(Zi)[Ir])

s = np.array(s) - 1

The following two figures shows the empirical PMF for $r=1$ and $r=50$ respectively. To me, they look quite like the beta binomial distributions, which might make sense as we are dealing with order statistics. I am able to fit the parameters for the distribution with reasonable accuracy, but I am more interested in how the parameters come about.

Answer 1

As Matt pointed out in the comments - it is indeed a hard problem. The answer below is a more credible attempt that (in my opinion) better approximates the quantity.

A general prior for the distribution

The reformulated question

Fix $r < N$, what is the distribution across all $s \in \{1, ..., N\}$, where $X_i$ generated $Y_{(r)}$ and $Z_i = X_i + \epsilon_{Zi}$ is the $s^{\textrm{th}}$ ranked item?

can be formalised as the probability that exactly $s-1$ other (assumed independent) $Z_k$ are less than $Z_i$, which makes $Z_i$ the $s^\textrm{th}$ ranked item:

$$\mathbb{P}\left(\sum_{k=1, k\neq i}^{N} \mathbb{I}_{\{Z_k < Z_i\}} = s - 1\right).$$

Let's call the sum of the indicator variables $C$ (for count).

To get the probability mass function of $C$, we first recognise that $Z_i$ is a continuous random variable with probability density $f_{Z_i}(z)$. Also for any realisation $z$, the probability that an independent $Z_k$ is less than $z$ is simply ${p = F_{Z_k}(z)}$, where $F_{Z_k}$ is the cumulative density function of $Z_k$.

Hence, the resultant probability, let's call it $P$, is a random variable (as opposed to a fixed value) that results from transforming $Z_i$ using $F_{Z_k}$, with probability density \begin{align} f_{P}(p) = f_{Z_k}(z) \left|\frac{\textrm{d}z}{\textrm{d}p}\right| = \frac{f_{Z_i}(F_{Z_k}^{-1}(p))}{f_{Z_k}(F_{Z_k}^{-1}(p))} . \end{align} The distribution can then be used as a prior for $C$, which has a binomial likelihood.

Simplifying the prior to something we can work with

Now the problem here is that the prior above is likely to be intractable, even under the normal assumptions as specified in the original question. This is because $Z_k$ is normal, but $Z_i$ is almost certainly not as it carries ranking information from $Y_{(r)}$.

In order to have something that we can actually work with, we approximate the prior using beta distributions. Beta distributions are a natural choice to me as they are closely related to order statistics, and moreover is a conjugate prior to binomial likelihood, which eases the computation of the probability masses for $C$. The plots in the question seems to support such hypothesis as well.

The key idea is to find the beta distribution parameters using method of moments. We first need the mean and variance of $Z_i$ (before transforming it with $F_{Z_k}$), then the mean and variance of $P$ (post transformation, denoted $\mu_{P}$ and $\sigma^2_{P}$ respectively), and finally fit the beta parameters.

The mean and variance of $Z_i$ (before transforming it with $F_{Z_k}$) is:

$$ \mathbb{E}(Z_i) \approx \frac{\sigma^2_X}{\sqrt{\sigma^2_X + \sigma^2_Y}} \Phi^{-1}\left(\frac{r - \alpha}{N - 2\alpha + 1}\right) , \quad \alpha \approx 0.4$$

$$ \textrm{Var}(Z_i) \approx \frac{\sigma^2_1 \sigma^2_X}{\sigma^2_X + \sigma^2_1} \, + \frac{\sigma^4_X}{\sigma^2_X + \sigma^2_1} \frac{r(N-r+1)}{(N+1)^2 (N+2)} \frac{1}{\big(\phi\big(\Phi^{-1}\big(\frac{r}{N+1}\big)\big)\big)^2} + \sigma^2_2 $$ (see this question, this question, [1] and [2] for detailed derivations).

The expected value and variance of ${P = F_{Z_k}(Z_i)}$ can then be approximated using Taylor series expansion: \begin{align} \mu_{P} \approx \, & F_{Z_k}\left(\mathbb{E}(Z_i)\right) + \frac{1}{2} f'_{Z_k}\left(\mathbb{E}(Z_i)\right) \cdot \textrm{Var}(Z_i) + \, ... \,,\\ \sigma^2_{P} \approx \,& \left(f_{Z_k}(\mathbb{E}(Z_i))\right)^2 \cdot \textrm{Var}(Z_i) + \frac{1}{4}\left(f'_{Z_k}(\mathbb{E}(Z_i))\right)^2 \cdot \textrm{Var}\left(\left(Z_i - \mathbb{E}(Z_i)\right)^2\right) + \, ... \,, \end{align}

where $f_{Z_k}$ is the PDF of $Z_k$ (normal w/ mean 0, variance $\sigma^2_X + \sigma^2_Z$), and $f'_{Z_k}$ is the derivative of the PDF.

Finally, the beta distribution parameters can be fit as follow (see this question): \begin{align} \alpha_{P} = \left(\frac{1 - \mu_{P}}{\sigma^2_{P}} - \frac{1}{\mu_{P}} \right) \mu_{P}^2 \;, \quad \beta_{P} = \alpha_{P}\left(\frac{1}{\mu_{P}} - 1\right) \;. \end{align}

---

[1] H. A. David and H. N. Nagaraja (2004) Order statistics. Encyclopedia of Statistical Sciences.
[2] C. H. B. Liu and B. P. Chamberlain (2019) What is the value of experimentation & measurement? ICDM 2019.