裏 RjpWiki

文字通り,RjpWiki の裏を行きます
Julia, R, プログラム, コンピュータ・サイエンス, 統計学

Brunner-Munzel の並べ替え検定

2015年02月20日 | ブログラミング

http://oku.edu.mie-u.ac.jp/~okumura/stat/brunner-munzel.html

で紹介されているが,結構時間がかかる。

> system.time(print(okumura(x, y)))
[1] 0.008037645
   ユーザ   システム       経過  
  1898.105     10.102   1923.545

brunner.munzel.test の必要な部分だけ取り出して,2倍強の速度になった。

> system.time(print(brunner.munzel.permutation.test(x, y)$p.value))
[1] 0.008037645
   ユーザ   システム       経過  
   788.666      5.037    794.799

brunner.munzel.permutation.test = function(x, y) {
    BM = function(X) {
        x = xandy[X]
        y = xandy[-X]
        r1 = rank(x)
        r2 = rank(y)
        r = rank(c(x, y))
        m1 = mean(r[1:n1])
        m2 = mean(r[n1 + 1:n2])
        v1 = sum((r[1:n1] - r1 - m1 + (n1 + 1)/2)^2)/(n1 - 1)
        v2 = sum((r[n1 + 1:n2] - r2 - m2 + (n2 + 1)/2)^2)/(n2 - 1)
        (m2 - m1)/sqrt(n1 * v1 + n2 * v2)
    }
    x = na.omit(x)
    y = na.omit(y)
    n1 = length(x)
    n2 = length(y)
    xandy = c(x, y)
    structure(list(method = "Brunner-Munzel Permutation Test",
        p.value = mean(abs(combn(n1+n2, n1, BM)) >= abs(BM(1:n1))),
        data.name = paste(deparse(substitute(x)), "and", deparse(substitute(y)))),
        class = "htest")
}

okumura = function(x, y) {
    bm = brunner.munzel.test(x, y)$statistic
    n1 = length(x)
    n2 = length(y)
    N = n1 + n2
    xandy = c(x, y)
    foo = function(X) {
        brunner.munzel.test(xandy[X], xandy[-X])$statistic
    }
    z = combn(1:N, n1, foo)
    mean(abs(z) >= abs(bm))
}

コメント   この記事についてブログを書く
« こりもせず,単純な暗号 | トップ | How to Transition from Exce... »
最新の画像もっと見る

コメントを投稿

ブログ作成者から承認されるまでコメントは反映されません。

ブログラミング」カテゴリの最新記事