8

我有一个矩阵，其中每一行都是来自分布的样本。我想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。

我尝试过但无法正常工作的方法：dplyrusing rowwise/do/lag、zoousing rollapply（这是我用来生成分布的方法）和data.table循环填充 a （编辑：这个有效，但仍然很慢）。

4

4 回答 4

7

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;
}

对于 size 的矩阵(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
于 2015-04-24T18:17:09.027 回答
3

加速的一个来源是编写一个较小的版本，它的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
于 2015-04-24T17:27:38.087 回答
2

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

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循环更快，但就性能而言，它不太可能是最佳的整体解决方案。

于 2015-04-24T17:33:01.113 回答
1

这是一个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
于 2015-04-24T17:29:15.910 回答