# Value at $D_\max$ from Kolmogorov-Smirnov test in R

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!

• Have you drawn the q-q plot?
– whuber
May 21 '13 at 18:34
• I have. Can it be extracted from there? May 21 '13 at 18:58
• 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.
– whuber
May 21 '13 at 19:13
• @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. May 21 '13 at 19:37
• How about simulating Dn for one sample KS test? Oct 3 at 22:23

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

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")  • 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. May 22 '13 at 19:07 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.