7
$\begingroup$

I am comparing two distributions using the Kolmogorov-Smirnov test [ks.test()] in R and would like to know what the numerical value is where $D_\text{max}$ occurs. It's not in the output of ks.test() as far as I can see, so I was wondering if anyone had any ideas about how I could find that out. Thanks!

$\endgroup$
6
  • 1
    $\begingroup$ Have you drawn the q-q plot? $\endgroup$
    – whuber
    May 21 '13 at 18:34
  • $\begingroup$ I have. Can it be extracted from there? $\endgroup$ May 21 '13 at 18:58
  • $\begingroup$ No (qqplot merely returns both arrays in sorted order), but the plot makes the location of Dmax visually evident. If you need it numerically, you can obtain it by copying the first page (or so) of code for ks.test. $\endgroup$
    – whuber
    May 21 '13 at 19:13
  • $\begingroup$ @whuber I appreciate the suggestion. I'm looking into the code for the ks.test and having a hard time interpreting where this value may be. Thanks for pointing me in this direction though. $\endgroup$ May 21 '13 at 19:37
  • $\begingroup$ How about simulating Dn for one sample KS test? $\endgroup$
    – RLOUIS
    Oct 3 at 22:23
8
$\begingroup$

Something like that? Dmax occurs at the value max.at.

set.seed(12345)

x <- rnorm(10000, 5, 5)
y <- rnorm(10000, 7, 6.5)

# remove any missings from the data

x <- x[!is.na(x)]
y <- y[!is.na(y)]

ecdf.x <- ecdf(x)
ecdf.y <- ecdf(y)

plot(ecdf.x, xlim=c(min(c(x,y)), max(c(x,y))), verticals=T, cex.lab=1.2, cex.axis=1.3,
     las=1, col="skyblue4", lwd=2, main="")

plot(ecdf.y, verticals=T, add=T, do.points=FALSE, cex.lab=1.2,
     cex.axis=1.3, col="red", lwd=2)

n.x <- length(x)
n.y <- length(y)

n <- n.x * n.y/(n.x + n.y)
w <- c(x, y)

z <- cumsum(ifelse(order(w) <= n.x, 1/n.x, -1/n.y))

max(abs(z)) # Dmax
[1] 0.1664

ks.test(x,y)$statistic # the same
     D 
0.1664

max.at <- sort(w)[which(abs(z) == max(abs(z)))]
[1] 9.082877

# Draw vertical line

abline(v=max.at, lty=2)

lines(abs(z)~sort(w), col="purple", lwd=2)

legend("topleft", legend=c("x", "y", "|Distance|"), col=c("skyblue4", "red", "purple"), lwd=c(2,2,2), bty="n")

Dmax visualization

$\endgroup$
1
  • 2
    $\begingroup$ Thanks for this piece of code @COOLSerdash. I had to get rid of some duplicates and NA's from my data set to get this to work, but it appears to be working now. $\endgroup$ May 22 '13 at 19:07
2
$\begingroup$

You could also use @COOLSerdash's answer plus environments to make the ks.test function output the value directly, like this:

ks.test.2 <- function(x, y, ..., alternative = c("two.sided", "less", "greater"), 
    exact = NULL) 
{
e <- new.env()
ks.test.2 <- ks.test
environment(ks.test.2) <- e
e$C_pkstwo <- stats:::C_pkstwo
    e$C_psmirnov2x <- stats:::C_psmirnov2x
e$C_pkolmogorov2x <- stats:::C_pkolmogorov2x
    e$return <- function(x){
    w<-get("w", envir=parent.frame())
    z<-get("z", envir=parent.frame())
    x$max.at <- sort(w)[which(abs(z) == max(abs(z)))]
    return(x)
    }
ks.test.2(x, y, ..., alternative = c("two.sided", "less", "greater"), 
    exact = NULL)
}

The function ks.test.2 should behave exactly like ks.test except that now it also returns the desired max.at component.

set.seed(12345)

x <- rnorm(10000, 5, 5)
y <- rnorm(10000, 7, 6.5)

ks.test.2(x,y)$max.at
# [1] 9.082877

This is only for the two-sided alternative, but you could enhance it to deal with the one-sided alternative if desired.

$\endgroup$

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.