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

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)

• You probably already realize but for continuous distributions such a test is distribution-free (i.e. its distribution under the null doesn't depend on the common distribution you're drawing from). As a result you can assume X and Y are both drawn from any convenient distribution (such as the standard uniform, say), for obtaining critical values, if that makes it easier to think about. I think the test statistic will turn out to be the negative hypergeometric. Jul 7, 2017 at 18:13
• If you want to determine whether the tails differ, you can usually find more-robust nonparametric tests by choosing a small integer $r$ and counting $y$, the number of $Y$ values falling within the top $r$ values of the combined batch of data. Significance is determined by $r$, $y$, $n_Y$, and $n_X$. You select $r$ depending on how much of the tail you are concerned with (before examining the data). This is an example of a "nonparametric slippage test."
– whuber
Jul 7, 2017 at 19:11

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.

• Excellent, thank you @Glen_b for the very complete and well explained answer! Jul 7, 2017 at 19:32
• As you mentioned, the choose() function in R will fail at n's as large as I'm using. I tried chooseZ() from the gmp package. It works, but it is slower than your solution with lchoose. exp(lchoose(...) - lchoose(...)) is very very fast. Jul 7, 2017 at 19:54
• Indeed; I implemented it just that way when I was checking I hadn't made an error in my reasoning/algebra: function(s,nx,ny) exp(lchoose(nx+ny-1-s,ny-s)-lchoose(nx+ny,ny)) ... if you get up high enough in $n$s or out very very far in the tail it's possible you may even need to take some precautions in relation to how you sum the results (such as summing from the smaller terms to the larger ones), but it would need to be a fair bit more extreme than your example situation. .... NB I have made some small additions at the end of my answer. Jul 7, 2017 at 23:57