Julia に翻訳--099 ファン・デル・ワーデン検定

#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。

ファン・デル・ワーデン検定
http://aoki2.si.gunma-u.ac.jp/R/vdw-test.html

ファイル名: vdwtest.jl　　関数名: vdwtest

翻訳するときに書いたメモ

x = vcat(x, y) を append!(x, y) にすると，呼出元の x が変更されてしまう。

==========#

using Rmath, StatsBase

function vdwtest(x, y)
    n1 = length(x)
    n2 = length(y)
    n = n1 + n2
    x = vcat(x, y) # 呼出元に影響を及ぼさない
    x = qnorm.(tiedrank(x) ./ (n + 1))
    S = abs(sum(x[1:n1]))
    V = n1 * n2 / (n ^ 2 - n) * sum(x .^ 2)
    Z = S / sqrt(V)
    pvalue = pnorm(Z, false) * 2
    println("S = $S,  V = $V,  Z = $Z,  p value = $pvalue")
    Dict(:S => S, :V => V, :Z => Z, :pvalue => pvalue)
end

gr1 = [1.2, 1.9, 2.5, 6.7];
gr2 = [1.5, 3.1, 10.5];
vdwtest(gr1, gr2)
# S = 0.7944989941443014,  V = 1.0741549510392339,  Z = 0.7665842364210947,  p value = 0.44332875078500844

Julia に翻訳--098 マン・ホイットニーの U 検定

#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。

マン・ホイットニーの U 検定
http://aoki2.si.gunma-u.ac.jp/R/u-test.html

ファイル名: utest.jl
関数名:    utest(tbl; correct = true)
          utest(x, y; correct = true)

翻訳するときに書いたメモ

==========#

using StatsBase, Rmath

function utest(tbl; correct = true)
    nr, nc = size(tbl)
    nr == 2 || error("2 x m matrix")
    x = []
    y = []
    for i = 1:nc
        append!(x, repeat([i], tbl[1, i]))
        append!(y, repeat([i], tbl[2, i]))
    end
    utest(x, y; correct)
end

function utest(x, y; correct = true)
    n1 = length(x)
    n2 = length(y)
    n = n1 + n2
    xy = vcat(x, y)
    r = tiedrank(xy)
    U1 = n1 * n2 + n1 * (n1 + 1) / 2 - sum(r[1:n1])
    ties = calctie(r)
    U = min(U1, n1 * n2 - U1)
    V = n1 * n2 *(n ^ 3 - n - ties) / 12 / (n ^ 2 - n)
    E = n1 * n2 / 2
    Z = (abs(U - E) - (correct ? 0.5 : 0)) / sqrt(V)
    P = pnorm(Z, false) * 2
    println("U = $U,  E(U) = $E,  V(U) = $V,  Z = $Z,  p value = $P")
    Dict(:U => U, :EU => E, :VU => V, :Z => Z, :pvalue => P)
end

function calctie(r)
    t = [count(isequal(s), r) for s in Set(r)]
    sum(t .^ 3 .- t)
end

二つのデータベクトルを与えるとき

A = vcat(repeat([1], 9), repeat([2], 12), repeat([3], 6), repeat([4], 3))
B = vcat(repeat([1], 4), repeat([2], 9), repeat([3], 11), repeat([4], 5))
utest(A, B, correct=false)
# U = 310.5,  E(U) = 435.0,  V(U) = 3993.5593220338983,  Z = 1.9701045835900841,  p value = 0.04882638572318423

utest(A, B)
# U = 310.5,  E(U) = 435.0,  V(U) = 3993.5593220338983,  Z = 1.9621925169893206,  p value = 0.04974007481574249

2 x C の集計表を与えるとき

tbl = [10 5 1; 4 7 3]
utest(tbl, correct=false)
# U = 70.0,  E(U) = 112.0,  V(U) = 481.9862068965517,  Z = 1.913074951284128,  p value = 0.05573845794420809

utest(tbl)
# U = 70.0,  E(U) = 112.0,  V(U) = 481.9862068965517,  Z = 1.8903002494831267,  p value = 0.05871781545216507

Julia に翻訳--097 二次データで三群以上の等分散性の検定

#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。

三群以上の等分散性の検定（二次データ）
http://aoki2.si.gunma-u.ac.jp/R/my-bartlett-test.html

ファイル名: mybartletttest.jl　　関数名: mybartletttest

翻訳するときに書いたメモ

==========#

using Rmath

function mybartletttest(n, u)
    (length(n) == length(u) && all(n .> 1) && all(u .> 0) && all(floor.(n) .== n)) || error("argument error")
    ng = length(n)
    temp1 = n .- 1
    temp2 = sum(temp1)
    chisq0 = temp2 * log(sum(temp1 .* u) / temp2) - sum(temp1 .* log.(u))
    co = 1 + (sum(1 ./ temp1) - 1/temp2) / (3 * ng - 3)
    chisq = chisq0 / co
    df = ng - 1
    pvalue = pchisq(chisq, df, false)
    println("chisq. = $chisq,  df = $df,  p value = $pvalue")
    Dict(:chisq => chisq, :df => df, :pvalue => pvalue)
end

n = [11, 25, 22]
u = [23.7, 25.7, 24.1]
mybartletttest(n, u)
# chisq. = 0.03290799098798558,  df = 2,  p value = 0.9836806320912985

Julia に翻訳--096 一元配置分散分析，三群以上の平均値の差の検定，二次データ

#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。

一元配置分散分析（三群以上の平均値の差の検定，二次データ）
http://aoki2.si.gunma-u.ac.jp/R/my-oneway-ANOVA.html

ファイル名: myonewayanova.jl　　関数名: myonewayanova

翻訳するときに書いたメモ

==========#

using Rmath

function myonewayanova(n, m, u; varequal=false)
    ng = length(n)
    dfb = ng - 1
    (ng > 1 && length(m) == ng && length(u) == ng && all(n .> 1) && all(floor.(n) == n) && all(u .> 0)) || error("parameter error")
    if varequal
        nc = sum(n)
        sw = sum(u .* (n .- 1))
        sb = sum(n .* (m .- sum(n .* m)/nc) .^ 2)
        dfw = nc - ng
        msb, msw = sb / dfb, sw / dfw
        f = msb / msw
        pvalue = pf(f, dfb, dfw, false)
    else
        w = n ./ u
        sumw = sum(w)
        m0 = sum(w .* m) / sumw
        temp = sum((1 .- w ./ sumw) .^ 2 ./ (n .- 1)) / (ng^2-1)
        f = sum(w .* (m .- m0) .^ 2) / ((ng - 1) * (1 + 2 * (ng - 2) * temp))
        dfw = 1 / (3 * temp)
        pvalue = pf(f, dfb, dfw, false)
    end
    println("F = $f,  dfb = $dfb,  dfw = $dfw,  p value = $pvalue")
    return Dict(:F => f, :dfb => dfb, :dfw => dfw, :pvalue => pvalue)
end

n = [8, 11, 22, 6];
m = [135.83, 160.49, 178.35, 188.06];
SD = [19.59, 12.28, 15.01, 9.81];
u = SD .^ 2;

myonewayanova(n, m, u, varequal=true)
# F = 20.82824265645924,  dfb = 3,  dfw = 43,  p value = 1.7374840062619836e-8

myonewayanova(n, m, u)
# F = 17.56460925505331,  dfb = 3,  dfw = 16.50402543036444,  p value = 2.191921147555868e-5

Julia に翻訳--095 対応のある代表値の差の検定，ウィルコクソンの符号付順位和検定

#==========
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)