Julia に翻訳--126 テューキー法，平均値の線形比較

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

テューキーの方法による平均値の線形比較
http://aoki2.si.gunma-u.ac.jp/R/tukey2.html

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

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

Rmath の ptukey()，qtukey() にバグがある。
lower_tail = false を指定すると NaN を返す。
回避策としては，
ptukey(q, nmeans, df, false)
を
1 - ptukey(q, nmeans, df)
とする。

qtukey(p, nmeans, df, false)
を
qtukey(1-p, nmeans, df)
とする。

==========#

using Rmath

function tlc(n, m, u, g1, g2; conflevel = 0.95)
    all(length(n) == length(m) == length(u) &&
            all(n .> 1) &&
            all(u .> 0) &&
            all(floor.(n) .== n) &&
            all(n .== n[1]) &&
            all(floor.(g1) .== g1) &&
            all(floor.(g2) .== g2)) || error("parameter error")
    ng = length(n)
    nc = sum(n)
    dfw = nc - ng
    Vw = sum(u .* (n .- 1)) / dfw
    n1 = length(g1)
    n2 = length(g2)
    g0 = collect(1:ng)
    g0 = [i for i in g0 if !(i in g1)]
    g0 = [i for i in g0 if !(i in g2)]
    n0 = ng - n1 - n2
    weight = vcat([repeat([i], j)
                   for (i, j) in zip([1 / n1, - 1 / n2, 0], [n1, n2, n0])]...)
    weight = weight[sortperm(vcat([g1, g2, g0]...))]
    theta = sum(weight .* m)
    sq = sqrt(Vw / n[1])
    q = qtukey(conflevel, ng, dfw)
    p = 1 - ptukey(abs(theta / sq), ng, dfw)
    confint = theta .- [1, - 1] * q * sq
    println("theta = $theta, a = $ng, df = $dfw, pvalue = $p")
    println("$(100conflevel) percent confidence interval: $confint")
    Dict(:theta => theta, :a => ng, :df => dfw,
         :pvalue => p, :confint => confint)
end

n = repeat([20], 5)
m = [35.6, 32.7, 29.4, 27.9, 25.7]
u = [5.6, 4.6, 4.9, 3.7, 2.5]

tlc(n, m, u, 1:3, 4:5)
# theta = 5.766666666666664, a = 5, df = dfw, pvalue = 4.640164918967571e-10
# 95.0 percent confidence interval: [3.951633188914924, 7.581700144418404]

tlc(n, m, u, 1:2, 4:5)
# theta = 7.350000000000007, a = 5, df = dfw, pvalue = 4.6348025417586314e-10
# 95.0 percent confidence interval: [5.534966522248267, 9.165033477751747]