高效地执行按行分布测试

Question

我有一个矩阵，其中每一行是从分布中抽取的样本。我想使用ks.test进行分布的滚动比较，并在每种情况下保存测试统计量。从概念上实现这个最简单的方式是通过循环：

set.seed(1942)
mt <- rbind(rnorm(5), rnorm(5), rnorm(5), rnorm(5))

results <- matrix(as.numeric(rep(NA, nrow(mt))))

for (i in 2 : nrow(mt)) {

  results[i] <- ks.test(x = mt[i - 1, ], y = mt[i, ])$statistic

}

然而，我的真实数据有大约400列和300,000行的单个示例，而且我有很多这样的示例。因此，我想让它变快一些。Kolmogorov-Smirnov检验在数学上并不是很复杂，所以如果答案是“用Rcpp实现它”，那么我会勉强接受，但我会感到有些惊讶——在R中，对于单个对，计算速度已经非常快了。

我尝试过的方法，但无法使其正常工作：dplyr使用rowwise/do/lag、zoo使用rollapply(这是我用来生成分布的方法)，以及在循环中填充data.table(编辑：这个方法可以运行，但仍然很慢)。

- Ajar

---

你是否真正在使用KernSmooth包？ks.test函数在stats包中。 - davechilders

你说得对！我确实使用了KernSmooth，但不是用于这个函数——我是用它来生成分布的。我会进行编辑。 - Ajar

Answer 1

在Rcpp中进行快速而简单的实现

// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h> 

double KS(arma::colvec x, arma::colvec y) {
  int n = x.n_rows;
  arma::colvec w = join_cols(x, y);
  arma::uvec z = arma::sort_index(w);
  w.fill(-1); w.elem( find(z <= n-1) ).ones();
  return max(abs(cumsum(w)))/n;
}
// [[Rcpp::export]]
Rcpp::NumericVector K_S(arma::mat mt) {
  int n = mt.n_cols; 
  Rcpp::NumericVector results(n);
  for (int i=1; i<n;i++) {
    arma::colvec x=mt.col(i-1);
    arma::colvec y=mt.col(i);
    results[i] = KS(x, y);
    }
  return results;
}

对于大小为(400, 30000)的矩阵，它可以在1秒内完成。

system.time(K_S(t(mt)))[3]
#elapsed 
#   0.98 

结果似乎是准确的。

set.seed(1942)
mt <- matrix(rnorm(400*30000), nrow=30000)
results <- rep(0, nrow(mt))
for (i in 2 : nrow(mt)) {
  results[i] <- ks.test(x = mt[i - 1, ], y = mt[i, ])$statistic
}
result <- K_S(t(mt))
all.equal(result, results)
#[1] TRUE

- Khashaa

---

那很快。我要测试一下！ - Ajar

太快了，真是太厉害了。相比之下，我在大约两个小时后停止了我的rollapplyr()解决方案（此时它已经生成了几乎所有的结果，但仍在运行）。它的结果与ks.test()匹配吗？ - Alex A.

我没有检查准确性，因此使用了标识符“dirty”。 - Khashaa

不完全相同，但非常接近：all.equal(results.ks2, results.cpp[2:280007]) [1] "Mean relative difference: 7.642923e-05"。而且在我的实际数据上，它比ks.test2快大约9倍。 - Ajar

鉴于性能和可接受的准确性，我认为这很可能是您最好的解决方案，@Ajar。 - Alex A.

确实。感谢大家的精彩发言！ - Ajar

Answer 2

加速的一个途径是编写一个较小版本的ks.test，其功能更少。下面的ks.test2比ks.test更为严格。例如，它假设您没有缺失值，并且总是希望获得与双侧检验相关的统计量。

ks.test2 <- function(x, y){

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

  max(abs(z))

}

请确认输出结果与ks.test一致。

set.seed(999)
x <- rnorm(400)
y <- rnorm(400)

ks.test(x, y)$statistic

    D 
0.045

ks.test2(x, y)

[1] 0.045

现在确定来自较小函数的节省：

library(microbenchmark)

microbenchmark(
  ks.test(x, y),
  ks.test2(x, y)
  )

Unit: microseconds
           expr      min       lq      mean   median        uq      max neval cld
  ks.test(x, y) 1030.238 1070.303 1347.3296 1227.207 1313.8490 6338.918   100   b
 ks.test2(x, y)  709.719  730.048  832.9532  833.861  888.5305 1281.284   100  a 

- davechilders

---

我很想看到使用这个函数代替ks.test()的rollapplyr()解决方案的基准测试结果。当前基准测试完成后，我会进行测试。 - Alex A.

我也非常感兴趣！我目前正在测试一些这些答案。 - Ajar

Answer 3

我能够使用rollapplyr()和ks.test()计算成对的Kruskal-Wallis统计量。

results <- rollapplyr(data = big,
                      width = 2,
                      FUN = function(x) ks.test(x[1, ], x[2, ])$statistic,
                      by.column = FALSE)

这样可以得到预期的结果，但对于您的数据集来说速度很慢。非常慢。这可能是因为ks.test()在每次迭代时计算的不仅仅是统计量；它还需要计算p值并进行大量的错误检查。

事实上，如果我们像下面这样模拟一个大型数据集：

big <- NULL
for (i in 1:400) {
    big <- cbind(big, rnorm(300000))
}

rollapplyr() 的方法需要很长时间; 大约 2 小时后我停止了执行，此时它已经计算出了几乎所有（但不是全部）的结果。

看起来，虽然 rollapplyr() 可能比 for 循环更快，但它在性能方面可能不是最佳解决方案。

- Alex A.

Answer 4

这里有一个使用dplyr的解决方案，可以得到与循环相同的结果。我对于它是否比循环更快存疑，但或许它可以作为解决方案的第一步。

require(dplyr)
mt %>% 
  as.data.frame %>%
  mutate_each(funs(lag)) %>%
  cbind(mt) %>%
  slice(-1) %>%
  rowwise %>%
  do({
    x = unlist(.)
    n <- length(x)
    data.frame(ks = ks.test(head(x, n/2), tail(x, n/2))$statistic)
  }) %>%
  unlist %>%
  c(NA, .) %>%
  matrix

- shadow