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

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.

• "simulated realizations of the same distribution" what sort of simulation? Like a bootstrap distribution created from the data? Feb 17 at 0:09
• Can I simplify this question as drawing two samples of size $n$ and $m$ from a uniform distribution, 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? What's the deal with the $Q$ number of realisations? Is it that a realisation can have some error. Is a sample from a realisation different from a sample from the original distribution? Feb 17 at 0:12
• @SextusEmpiricus Yes, for the purposes of this question they can be treated as two samples of different sizes drawn from the same distribution. Perhaps I shouldn't have even used the word "simulation" here - the context is that they're statistics computed on astrophysical simulations of real phenomena, but that's probably ancillary with respect to my question. $Q$ is just an ingredient to turn the "probability of exceeding largest real" to "number exceeding largest real", so is also not a focus. Feb 17 at 1:05
• "probability of exceeding largest real" that made me think of a third approach. The largest real is beta distributed. Conditional on the largest real, the number of cases that exceed the largest real is binomial distributed. So you can see it as a compound distribution. Feb 17 at 7:50

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.