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

}

ただし、私の実際のデータには、1 つの例で最大 400 列と最大 300,000 行があり、多くの例があります。だから早くしてほしい。コルモゴロフ-スミルノフ検定は数学的にそれほど複雑ではないので、答えが「で実装する」である場合、Rcpp私はそれをしぶしぶ受け入れますが、少し驚いています.1つのペアで計算するのはすでに非常に高速です. Rで。

私が試したがうまくいかなかった方法: dplyrusing rowwise/do/lag、zoousing rollapply(これは私がディストリビューションを生成するために使用するものです)、およびdata.tablea をループに投入します (編集: これは機能しますが、それでも遅いです)。

4

4 に答える 4

3

スピードアップの 1 つのソースは、より少ない機能を備えた のより小さなバージョンを作成すること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 時間後に実行を停止しました。この時点で、ほぼすべての (すべてではない) 結果が計算されました。

whileはループ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 に答える