r - 有效地执行行式分布测试

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在一个循环中(编辑:这个工作，但它仍然很慢)。

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