Hypothesis Test: Bound for number of observations from Y that exceed max(X) if X=Y in distribution

Question

I believe this question is related to Extreme Value Theory, an area of statistics that I have not studied.

Let

• $X$ and $Y$ be random variables
• the null hypothesis be that $X=Y$ in distribution
• $\{X_i\}_{i=1}^{n_X}$ be a sample of size $n_X$ from $X$ and $\{Y_i\}_{i=1}^{n_Y}$ be a sample of size $n_Y$ from $Y$
• $X_{max} = max(\{X_i\}_{i=1}^{n_X})$
• $S = \sum_{i=1}^{n_Y} I(Y_i > X_{max}) $, where $I$ is the indicator function
• $n_Y >> n_X$ (e.g. $n_Y= 100000$ and $n_X=500$). This isn't a necessary assumption, but it is useful for understanding why I'm using $S$ as a test statistic.

The goal is to determine when $S$ is unusually high under the null hypothesis. Meaning, find $b$ such that we reject the null when $S>b$, where b is a sample statistic and $P(S > b) < \alpha$ under the null. For example the 99th quantile of a binomial of size $n_Y$ and probability $1/(n_X +1)$ (this bound can be verified to be wrong through simulation). Can a valid $b$ be determined? Or more generally:

How do we use $S$ to determine if the tails of $X$ and $Y$ are different?

Possible solutions:

• Determine the distribution of $S$ under the null hypothesis
• Something involving Order Statistics
• Something involving Extreme Value Theory

Answers may certainly include additional assumptions if needed for this to be tractable.

Here is a quick simulation in R (with the wrong bound) showing what I am looking for:

NSIM <- 1000
nX <- 500
nY <- 100000
alpha <- .01

exc <- rep(NA, NSIM)
for(i in 1:NSIM){
   cat('\r', paste0("Sim: ",i,"/", NSIM))
   X <- rnorm(nX)
   Xmax <- max(X)
   Y <- rnorm(nY)
   # b is the wrong bound
   b <- qbinom(p = 1-alpha, size = nY, prob = 1/(nX+1))
   S <-sum(Y > Xmax)
   exc[i] <- S > b
}
### This is the type I error rate
### should be around alpha=.01 if correct bound
sum(exc)/NSIM
### [1] 0.307 (usually like .27-.33)

Answer 1

Let us reason as follows.

If the null hypothesis is true, then in a combined sample, any of the $n_X+n_Y$ observations has the same chance to be labelled with $Y$ as any other observation.

Counting how many $Y$s stick out one end is like we have a deck of $n_X$ red cards and $n_Y$ white cards, and we deal cards off the top of a shuffled deck until we hit the first red card, and we count how many white cards there were before then.

So that would suggest that the pmf under the null is that of a negative hypergeometric distribution for the number of successes until the first failure, where there are $n_Y$ successes and $n_X$ failures.

I think that will boil down to:

$$P(S=s) = \frac{{n_X+n_Y-1-s} \choose {n_Y-s}}{{n_X+n_Y}\choose{n_Y}}$$

A quick plausibility check:

Consider $n_X=2, n_Y=3$. We can compute from the above formula:

  s       0    1    2    3
P(S=s)   0.4  0.3  0.2  0.1

Now let's try a simulation to check it:

nsim=1000000L
table(replicate(nsim,{a=runif(3);b=runif(2);sum(a>max(b))}))/nsim
res
       0        1        2        3 
0.400026 0.300800 0.199712 0.099462 

That looks like we should expect.

For a given $n_X$ and $n_Y$ you can use this to find the smallest value $s_\text{crit}$ that has $P(S\geq s_\text{crit})\leq \alpha$ and then reject for any observed $s$ that is at least that large.

[Note that to do the calculations with large arguments you want a function like R's lchoose (which computes the log of ${n} \choose {x}$, or failing that, at least something like its lgamma (the log of a gamma function).]

Alternatively, you can compute a p-value for some observed $s$ as $P(S\geq s)$.

It may often be more convenient to compute $P(S< s)$ and take its complement.

When $n_Y$ and $n_X$ are both very large you may be able to use a geometric approximation (with $p=\frac{n_X}{n_X+n_Y+1}$). That may at least be useful at least in figuring out about where to sum up to, to find a more accurate critical value from that approximate one.

From the look of it for your example $n$'s and $\alpha=0.01$ that approximation would work well, taking you just a few values past the required quantile; it's easy to take the cumulative sum of the pmf up to there and so have the cdf to good accuracy.

Note that if you choose not the number of $Y$'s past the largest $X$-value, but above say the $k$-th largest $X$ (e.g. the number of $Y$'s past the tenth-highest $X$), that should still be negative hypergeometric.

I would advise you to consider the power properties of this test for distributions that look something like the data you have. Tests that look similar to this one may have great power in some situations but relatively poor power in others. In particular, if the upper tail is heavier than exponential, the power is likely to be quite poor, but if the distribution has a very light upper tail, the power may be quite good. Simulation to check you can have a reasonable chance to reject the null when you think you should be able to would be advisable.