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!
2 Answers
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")
-
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, 2013 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.
qqplot
merely returns both arrays in sorted order), but the plot makes the location ofDmax
visually evident. If you need it numerically, you can obtain it by copying the first page (or so) of code forks.test
. $\endgroup$