3
$\begingroup$

I have two datasets $S_1$ and $S_2$. I can run KS 2 samples test on these datasets to obtain the value of the KS 2 samples test.

Is there a correct way to bootstrap the KS 2 sample test using the datasets $S_1$ and $S_2$?

In other words, how exactly should sample from $S_1$ and $S_2$ to obtain a bootstrap distribution for the KS 2 sample test.

For example if we have in scipy the following 2 datasets :

from scipy import stats

n1 = 200  # size of first sample
n2 = 300  # size of second sample

S1 = stats.norm.rvs(size=n1, loc=0., scale=1)
S2 = stats.norm.rvs(size=n2, loc=0.5, scale=1.5)
stats.ks_2samp(S1, S2)

Then stats.ks_2samp returns a KS value. I would like to know how could we do bootstrap the samples $S_1$ and $S_2$ to obtain a bootstrap distribution for the ks_2samp values. Any help is appreciated.

$\endgroup$
  • 1
    $\begingroup$ What are you trying to accomplish? Why isn't enough to have just one K-S test comparing the two samples? I could see the purpose of a simulation procedure if you were trying to find the power of the the K-S test in certain circumstances. $\endgroup$ – BruceET Jun 16 at 6:56
  • 2
    $\begingroup$ @BruceET Besides power studies, I have found a few additional reasons to carry out such a procedure. One is where the two-sample test available does not correctly compute the p-value; for instance, maybe it uses a doubtful approximation or it simply refuses to run with datasets of a certain size. The other reason is where either of the datasets has a tie, which makes the KS-test inapplicable--but the KS statistic may still be useful for distinguishing the underlying distributions. $\endgroup$ – whuber Jun 16 at 14:33
  • $\begingroup$ @whuber , Do you have an idea for how to carry how bootstrap in this case? thanks. $\endgroup$ – Steve Jun 16 at 15:56
6
$\begingroup$

The null hypothesis is that the batches of data are independent simple random samples from a common continuous distribution. The (two-sample) Kolmogorov-Smirnov (KS) statistic measures a difference in their empirical distributions. Thus, exactly as in any other case where you have a measure of difference, you can bootstrap it by resampling from the combined dataset.

There's a slight problem: about $1/e\approx 37\%$ of a bootstrap sample will be duplicated, which is bad news for the KS test: many implementations will complain. One workaround is to add a tiny amount of noise to each resampled value. This is legitimate because adding a tiny bit of noise will scarcely change the underlying distribution (provided it has no discrete components).

(Another attractive alternative is to perform a permutation test in which one "bootstrap" sample is obtained without replacement from the combined dataset and is compared to the remainder of the data. If there are no ties the combined dataset, there will never be a problem with the KS statistic. The permutation test and bootstrap perform very nearly the same.)

Below is an example in R, commented and coded to port readily to Python. It works by concatenating the two datasets into an array S and sampling from this array (or, literally, from its indexes using sample.int). Its output documents bootstrapping five situations like the one in the question, giving some indication of the power of the test (that is, its likelihood of rejecting the null). A bootstrap of 500 resamples would be more than enough, but here I show 5000 resamples to create precise results.

Figure

(The code will perform a permutation test by setting do.bootstrap to FALSE. If nothing else is modified, the same random datasets will be generated, enabling direct comparison between the bootstrap and permutation test results.)

ppts <- function(n) (1:n - 1/2)/n
stat <- function(S1, S2, tol=1e-4) {
  sigma <- sd(c(S1,S2)) * tol # (See the text)
  X1 <- S1 + rnorm(length(S1), 0, sigma)
  X2 <- S2 + rnorm(length(S2), 0, sigma)
  ks.test(X1, X2)$statistic
}

do.bootstrap <- FALSE # Alternatively, use a permutation test.

par(mfcol=c(2,5))
set.seed(17)
seeds <- round(runif(5)*2^31)
n1 <- 10
for (i in 1:5) {
  #
  # Generate data.
  #
  set.seed(seeds[i])
  n2 <- round(n1 * 3/2)
  S1 <- rnorm(n1, 0, 1)
  S2 <- rnorm(n2, 0, 1.5)
  #
  # Compare the data graphically.
  #
  qqplot(S1, S2, bty="n", main="QQ Plot")
  abline(0:1)
  #
  # Bootstrap the statistic.
  #
  S <- c(S1, S2)
  boot <- replicate(5000, {
    if (do.bootstrap) {
      i <- sample.int(length(S), length(S), replace=TRUE)
      stat(S[i[1:n1]], S[i[-(1:n1)]])
    } else {
      i <- sample.int(length(S), n1, replace=FALSE)
      stat(S[i], S[-i], tol=0)
    }
  })
  #
  # Display the bootstrap results.
  #
  ks.data <- stat(S1, S2)
  boot <- c(ks.data, boot)
  p.value <- mean(boot >= ks.data)
  hist(boot, freq=FALSE, breaks=30, col="#f0f0f0", border="Gray",
       main=paste0("P-value is ", signif(p.value, 3)))
  abline(v = ks.data, lwd=2, col="Red")
}
par(mfcol=c(1,1))
| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ thanks man, this is pretty helpful. $\endgroup$ – Steve Jun 16 at 18:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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