2
$\begingroup$

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)

$\endgroup$
2
  • $\begingroup$ Power depends on sample size. Think about the equation for the t-statistic in a t-test: $t=\dfrac{\bar{x}-\mu_0}{s/\sqrt{n}}$. As $n$ increases, $t$ increases, and power increases. $\endgroup$
    – Dave
    Apr 13, 2021 at 22:42
  • $\begingroup$ Can you say explicitly what two normal distributions you are using? Long time since I used Matlab, but I can't match your 'offset' description with the code. Clarification might encourage more responses. // I did brief simulations in R and found somewhat similar jumps. $\endgroup$
    – BruceET
    Apr 14, 2021 at 5:52

2 Answers 2

1
$\begingroup$

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")

enter image description here

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)

enter image description here

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)

enter image description here

$\endgroup$
-1
$\begingroup$

@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.
KS-Test rejection rate. x-axis: Number of observations. y-axis: Rejection rate of H_0 (in %).

$\endgroup$
1
  • 1
    $\begingroup$ This answer is better as a follow-on comment to @Dave's comment. Please reserve answer space for true answers to the question. $\endgroup$
    – R Carnell
    Apr 13, 2021 at 23:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.