Issue with Kolmogorov - Smirnov test implementation in Matlab?

Question

I observed a strange behavior of the 2 sample Kolmogorov - Smirnov test in Matlab and I'm not sure if I missing the obvious or if there is indeed an issue with implementation of the test?

Below is a code example. I estimate the power of the test (counting how often H_0 is rejected) as a function of the number of observations (from 10 to 50). I take the observations from two normal distributions which are shifted by 1.5 units in their mean and then apply the built-in KS-Test ("kstest2").

The resulting figure doesn't look right, but the power of the test depends a lot on the sample size in a non-trivial sense. For example: With 14 observations the rejection rate is 51%, with 15 observations it is 38% and with 16 observations it is again 47%. This issue remains when I increase the sample size and I suspect it is related how the p-values are approximated within the kstest2 function. However, according to the Matlab description, the asymptotic assumption should be good for n > 8.

Do yo have any ideas on this?

n_samples = 10000;
n_obs     = 10:50;  
h_vec     = NaN(length(n_obs),n_samples);

for effi = 1:length(n_obs)
    for casi = 1:n_samples
        
      vec1 = normrnd(0,1,n_obs(effi),1);
      vec2 = normrnd(1,1.5,n_obs(effi),1); 
            
      h_KS = kstest2(vec1,vec2); % KS Test
     
      h_vec(effi,casi) = h_KS; 
      
    end
end 

plot(n_obs,100*sum(h_vec,2)./n_samples)

Answer 1

Comment continued. I used simplest possible R code to get power of K-S tests comparing $\mathsf{Norm}(\mu=0, \sigma=1)$ with $\mathsf{Norm}(\mu=1.5, \sigma=1)$ for sample sizes $n = 12, 13,\dots, 17.$ Power values showed jumps: 81.5%, 88.4%, 82.5%, 88.6%, 92.8%, 95.6%. So jumps of the kind you noted for small $n$ are not just an artifact of Matlab.

set.seed(2021)
pv.12 = replicate(10^5, 
         ks.test(rnorm(12,0,1),rnorm(12,1.5,1))$p.val)
mean(pv.12 <= .05)
[1] 0.81534
pv.13 = replicate(10^5, 
         ks.test(rnorm(13,0,1),rnorm(13,1.5,1))$p.val)
mean(pv.13 <= .05)
[1] 0.88418
pv.14 = replicate(10^5, 
         ks.test(rnorm(14,0,1),rnorm(14,1.5,1))$p.val)
mean(pv.14 <= .05)
[1] 0.82486
pv.15 = replicate(10^5, 
         ks.test(rnorm(15,0,1),rnorm(15,1.5,1))$p.val)
mean(pv.15 <= .05)
[1] 0.88563
pv.16 = replicate(10^5, 
         ks.test(rnorm(16,0,1),rnorm(16,1.5,1))$p.val)
mean(pv.16 <= .05)
[1] 0.92805
pv.17 = replicate(10^5, 
         ks.test(rnorm(17,0,1),rnorm(17,1.5,1))$p.val)
mean(pv.17 <= .05)
[1] 0.95564

I think the reason is the discreteness of the K-S test statistic for small sample sizes. It is the largest vertical discrepancy between ECDFs of the two samples compared. With small sample sizes there are not many 'jumps' in the ECDFs, hence not many possible values of the test statistic $D.$

set.seed(413)
x1 = rnorm(14, 0, 1);  x2 = rnorm(14, 1.5, 1)
ks.test(x1, x2)

    Two-sample Kolmogorov-Smirnov test

data:  x1 and x2
D = 0.64286, p-value = 0.0049
alternative hypothesis: two-sided

plot(ecdf(x1), col="blue", xlim=c(min(x1),max(x2)), 
     main="ECDFs of Two Normal Samples")
 lines(ecdf(x2), col="brown")

Note: It may also be worth mentioning that because of the discreteness of the test statistic for small n, the true significance level of a K-S test at the "5%" level may not be near 5%, and that the distribution of the P-value under $H_0$ is hardly uniform. Thus unusual behavior of power levels is not surprising.

pv.14 = replicate(10^5, 
 ks.test(rnorm(14,0,1),rnorm(14,0,1))$p.val)
mean(pv.14 <= .05)
[1] 0.01823   # signif level not 5%

hdr = "P-value of K-S Test Under Null Hypothesis"
hist(pv.14, prob=T, col="skyblue2", main=hdr)

By contrast, the pooled t test is well behaved when comparing normal samples with the same parameters.

pvt.14 = replicate(10^5, 
          t.test(rnorm(14,0,1),rnorm(14,0,1),var.eq=T)$p.val)
mean(pvt.14 <= .05)
[1] 0.04969  # signif level aprx 5%

hdr = "P-value of Pooled t Test Under Null Hypothesis"
hist(pvt.14, prob=T, col="skyblue2", main=hdr)

Answer 2

@Dave: Thanks, but that's clear and doesn't explain why there are jumps (both up and down) in the power when the sample size is varied.

This issue also doesn't occur with other tests in Matlab. If you replace "kstest2" with "ttest2" in the script above, the results look perfectly fine.
x-axis: Number of observations. y-axis: Rejection rate of H_0 (in %).