从置换数据集中生成大量向量的更快方法 [R]

Question

设置 出于模拟的目的，我生成了一个 B=2000 元素列表，每个元素都是排列过程的输出，在该排列过程中，我首先排列 @ 的行987654323@ 矩阵和对于每一列，我计算了第一和第二个 100 行之间的 Kolmogorov-Smirnov 检验统计量（您可以将前 100 行视为来自一个组和后 100 行作为另一个组的数据）。

问题此过程需要很长时间（大约 30-40 分钟）来生成列表。有更快的方法吗？以后我想把B 增加到更大的值。

代码

B=2000
n.row=200; n.col=8000
#Generate sample data
samp.dat = matrix(rnorm(n.row*n.col),nrow=n.row)

perm.KS.list = NULL
for (b in 1:B){
    #permute the rows
    perm.dat.tmp = samp.dat[sample(nrow(samp.dat)),]
    #Compute the permutation-based test statistics
    perm.KS.list[[b]]= apply(perm.dat.tmp,2,function(y) ks.test.stat(y[1:100],y[101:200]))
}


#Modified KS-test function (from base package)
ks.test.stat <- function(x,y){
  x <- x[!is.na(x)]
  n <- length(x)
  y <- y[!is.na(y)]
  n.x <- as.double(n)
  n.y <- length(y)
  w <- c(x, y)
  z <- cumsum(ifelse(order(w) <= n.x, 1/n.x, -1/n.y))
  z <- z[c(which(diff(sort(w)) != 0), n.x + n.y)] #exclude ties
  STATISTIC <- max(abs(z))
  return(STATISTIC)
}

Answer 1

1:B 循环有几个地方需要优化，但我同意真正的消费者是内部函数。因为您正在模拟表现良好的引导样本，所以您可以做出两个简化的假设，而一般 base 函数不能：

1. 没有缺失值。这避免了 is.na() 调整
2. 两侧（即x & y）的元素个数相同，因此无需分别计算。而不是在循环中拆分y，然后将它们重新加入函数中（进入w），只需将它们放在一起即可。平衡的边还允许简化，例如删除ifelse() 子句。它会产生一堆 0/1，然后用整数运算重新调整为 -1/1。

功能减少，节省约25%的时间。我在cumsum() 中添加了整数，而不是双精度数。

ks.test.stat.balanced <- function(w){
  n     <- as.integer(length(w) * .5)
  # z   <- cumsum(ifelse(order(w) <= n, 1L, -1L)) / n
  z     <- cumsum((order(w) <= n)*2L - 1L) / n
  # z   <- z[c(which(diff(sort(w)) != 0), n + n)] #exclude ties
  return( max(abs(z)) )
}

与你的高斯 rng 不应该经常发生关系，diff(sort(.)) 非常昂贵。如果您愿意取消该保护，则时间将减少约 65%。

如果您将z 的等式移动到abs()，它会比所有这些代表节省一点时间。我把它单独放在上面，这样更容易阅读。

编辑如果模拟不平衡，我建议你：

1. 还是把is.na拒之门外，
2. 还是通过w，
3. 还是尽量保留在integer，不是numeric，而是
4. 现在包括两个组大小的参数n1 和n2。

此外，在 cumsum() 之前预先计算 1/n 进行实验，以避免大量昂贵的除法。尝试考虑其他数学方法来从内部循环中提取计算，以减少它发生的频率。