Number of simulated statistics more extreme than most extreme real data statistic

Question

I have a distribution of $n_{obs}$ "real" data observations drawn from a normal distribution $X \sim N(\mu,\sigma^2)$, and a number $Q$ of simulated realizations of the same distribution, from which I have drawn $n_{obs}$ observations from each realization.

I am interested in knowing the distribution of the number of simulated observations greater than the maximum real data value.

My attempt:

We can take $\mu=0,\sigma^2=1$ for simplicity's sake.

The probability of a sim test statistic $x$ having a value larger than a value $M$ would be the tail distribution for $x=M$:

$P_G (x>M) = 1-\Phi(M)$

where $\Phi$ is the normal distribution CDF. The number of sims greater than $M$ would then be obtained by multiplying this probability by the number of sim statistics, i.e., $n_{obs} \cdot Q$.

Per probabilityislogic's answer to this question, $M$ above is given by the $n_{obs}$'th (the highest) order statistic, the distribution for which is given by a "beta-F" compound distribution. For the $i=N$ case in their notation, this comes out as:

$ P_M (x_{n_{obs}}) = n_{obs} \cdot f({x_{n_{obs}}})\Phi(x_{n_{obs}})^{n_{obs}-1}$

where $f$ is the normal distribution PDF here. I suppose I'd then want to marginalize over $x_{n_{obs}}$:

$P_{N_{gtr}}(x,n_{obs}) = \int P_G (x|x_{n_{obs}})P_M (x_{n_{obs}})dx_{n_{obs}}$

$ = \int (1-\Phi(x_{n_{obs}}))\cdot n \cdot f (x_{n_{obs}})\Phi(x_{n_{obs}})^{n-1}$

...which appears to be a horrific integral, though I can get Matlab to spit out a finite number with specified inputs. But Mathematica is failing to give me a symbolic representation, and from my numerical experiments, I'm suspecting I may have gone wrong in (at least) this last step.

I'm an astronomer who is a bit out of my depth with this... But would like to understand. It's been a while since I've properly studied formal statistics like this and I likely have some fundamental misconception(s). Approaching this in a brute-force empirical manner starting from pure Gaussian random variables in Matlab reveals that this produces a uniform distribution in the $n_{obs}=1$ case, and becomes exponentially distributed for the $n_{obs}\gg1$ case. I'd like to be able to reproduce or understand that limiting behavior, if possible.

If I numerically sample from $P_M$ and feed the results into $P_G$, I obtain something quite close to the brute-force numerical result, though with some strange scalings, so this suggests I'm at least in the neighborhood of what I want.

Answer 1

Depending on the details (if realisations are without error) then you can simplify this question as drawing two samples of size $n$ and $m= qn$, and express the distribution for the number of times that a variable from the second sample is larger than the maximum of the first sample.

You can also see this as first drawing a single sample of size $(q+1)n$, ordering them, and randomly assigning $n$ of the ranks to the first sample and $qn$ to the second sample.

The approaches below assume a continuous distribution (such that a probability for ties is zero).

---

Approach 1:

You can consider this assignment as an urn model.

• The probability that the highest rank is assigned to the first sample is $\frac{n}{(q+1)n}$
• The probability that the second highest rank is assigned to the first sample, given that the highest rank is not is $\frac{n}{(q+1)n-1}$
• The probability that the third highest rank is assigned to the the first sample, given that the higher rank are not is $\frac{n}{(q+1)n-2}$
• and so on

Then the probability of $k$ items from the second sample being higher than the maximum of the first sample is equal to the probability, is the probability of the event that the $k$ highest ranks are not being assigned to the first sample and the $k-1$ highest rank is being assigned to the first sample.

$$P(K=k) = \frac{n}{(q+1)n-k-1} \prod_{j=0}^{k} \left(1-\frac{n}{(q+1)n-j}\right)$$

---

Approach 2:

The assignment of the ranks is like drawing $n$ samples out of a discrete distribution of $(q+1)n$ numbers. This draw is without replacement, but for large $q$ we can approximate it as with replacement. Then the distribution is like $f(k)F(k)^{n-1}$, a geometric distribution, which simplifies the product from approach 1.

---

Approach 3:

The distribution of the quantile of the highest number from the first sample is beta distributed.

Given the quantile of the highest number from the first sample, the numbers with a higher quantile from the second sample that are higher, is binomial distributed.

Combining those two you have a beta-binomial distribution for the number of cases from the second sample that are higher than the maximum of the first sample.