#==========
Julia の修行をするときに,いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。
対応のある代表値の差の検定(ウィルコクソンの符号付順位和検定)
http://aoki2.si.gunma-u.ac.jp/R/wilcox-paired.html
ファイル名: wilcoxonsignedranktest.jl 関数名: wilcoxonsignedranktest
翻訳するときに書いたメモ
R の wilcox.test() に含まれる Wilcoxon Signed Rank Test では,同値がなくても n < 50 のときしか exact な検定を行えないようになっている。
しかし,使用される psignrank() は n < 1039 まで可能である。
このプログラムでは,可能な限り exact な検定を行い,同値がある場合には漸近検定を行う。exact な検定が行われた場合も,漸近検定の結果は表示される。
==========#
using StatsBase, Rmath
function calctie(r)
t = [count(isequal(s), r) for s in Set(r)]
sum(t .^ 3 .- t)
end
function wilcoxonsignedranktest(x, y; correct=true)
x .-= y
n0 = length(x)
x = x[x .!= 0]
n = length(x)
r = tiedrank(abs.(x))
w = sum(r[x .> 0])
ties = calctie(r)
z = w - n * (n + 1)/4
if n < 1039 && ties == 0 && n0 == n
pvalue = z > 0 ? psignrank(w - 1, n, false) : psignrank(w, n)
pvalue = min(2 * pvalue, 1)
println("w = $w, p value = $pvalue (exact)")
end
method = "asymptotic test"
sigma = sqrt(n * (n + 1) * (2 * n + 1)/24 - ties/48)
correction = 0
if correct
correction = 0.5sign(z)
method = method * " with continuity correction"
end
z = (z - correction) / sigma
pvalue = 2 * min(pnorm(z), pnorm(z, false))
println("z = $z, p value = $pvalue ($method)")
end
同値がある場合は漸近検定
x = [269, 230, 365, 282, 295, 212, 346, 207, 308, 257];
y = [273, 213, 383, 282, 297, 213, 351, 208, 294, 238];
wilcoxonsignedranktest(x, y)
# z = 0.0, p value = 1.0 (asymptotic test with continuity correction)
x = [48.1, 51.2, 61.9, 40.7, 49.2, 53.9, 61.0, 33.6, 21.5, 72.3]
y = [58.1, 63.2, 61.8, 39.0, 56.2, 63.9, 71.2, 62.9, 48.1, 66.6]
wilcoxonsignedranktest(x, y)
# z = -2.141909506235002, p value = 0.03220076475250472 (asymptotic test with continuity correction)
同値がない場合は常に exact な検定と漸近検定の両方の結果が表示される
x = [51.4, 42.2, 61.7, 47.9, 50.0, 63.6, 53.0, 34.2, 56.2, 49.9];
y = [43.8, 54.8, 51.6, 54.0, 61.0, 57.4, 67.0, 77.6, 65.9, 62.6];
wilcoxonsignedranktest(x, y)
# w = 10.0, p value = 0.083984375 (exact)
# z = -1.732800450887927, p value = 0.0831311427008039 (asymptotic test with continuity correction)
同値さえなければ,サンプルサイズが 1038 以下ならば exact な検定が行われる。
using Random; Random.seed!(123);
n = 1038
x = randn(n)
y = randn(n)
wilcoxonsignedranktest(x, y)
# w = 280908.0, p value = 0.24279165850242684 (exact)
# z = 1.1683136400265737, p value = 0.24268027537889492 (asymptotic test with continuity correction)