Julia に翻訳--122 平均値の多重比較，ライアン法，テューキー法

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

平均値の多重比較（ライアンの方法，テューキーの方法）
http://aoki2.si.gunma-u.ac.jp/R/m_multi_comp.html

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

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

Rmath の qtukey() にバグがある。

qtukey(q::Number, nmeans::Number, df::Number, nranges::Number=1.0,
    lower_tail::Bool=true, log_p::Bool=false) =
    ccall((:qtukey ,libRmath), Float64,
    (Float64, Float64, Float64, Float64, Int32, Int32),
    p, nranges, nmeans, df, lower_tail, log_p)
end

# 最初の引数 'q::Number' は q ではなく，p である。すなわち，'p::Number'
# しかし，それを直しても lower_tail = false を指定すると NaN を返す。
# 回避策としては，
# qtukey(p, nmeans, df, false)
# を
# qtukey(1-p, nmeans, df)
# とする。

==========#

using Rmath, Roots, Printf

function mmulticomp(n, me, s,; alpha = 0.05, method = "ryan") # or method = "tukey"
    k = length(n)
    all(k == length(me) && k == length(s) && all(n .> 0) && all(floor.(n) .== n) && all(s .> 0)) || error("parameter error")
    o = sortperm(me)
    sn = n[o]
    sm = me[o]
    ss = s[o]
    nt = sum(sn)
    mt = sum(sn .* sm) / nt
    dfw = nt - k
    vw = sum((sn .- 1) .* ss .^ 2) / dfw
    numsignificant = nsn = 0
    nss = zeros(Int, k * (k - 1) ÷ 2)
    nsb = zeros(Int, k * (k - 1) ÷ 2)
    for m = k:-1:2
        for small = 1:k - m + 1
            big = small + m - 1
            if check(small, big, nsn, nss, nsb)
                t0 = (sm[big] - sm[small]) / sqrt(vw * (1 / sn[big] + 1 / sn[small]))
                if method == "ryan"
                    P = pt(t0, dfw, false) * 2
                    nominalalpha = 2 * alpha / (k * (m - 1))
                    result = P <= nominalalpha
                else
                    t0 *= sqrt(2)
                    q = (qtukey(0.95, k, dfw) + qtukey(0.95, m, dfw)) / 2
                    WSD = q * sqrt(vw * (1 / sn[big] + 1 / sn[small]) / 2)
                    P = find_zero(p -> ((qtukey(1 - p, k, dfw) +
                                 qtukey(1 - p, m, dfw)) / 2 - t0), [0, 1], Bisection())
                    result = sm[big] - sm[small] >= WSD
                end
                if result
                    numsignificant = 1
                    @printf("mean[ % 2i] = %7.5f vs mean[ % 2i] = %7.5f : diff = % 7.5f, ",
                        o[small], me[small], o[big], me[big], me[big] - me[small])
                    if method == "ryan"
                        @printf("t = %7.5f : P = %7.5f, alpha' = %7.5f\n", t0, P, nominalalpha)
                    else
                        @printf("WSD = %7.5f : t = %7.5f : P = %7.5f\n", WSD, t0, P)
                    end
                else
                    nsn += 1
                    nss[nsn] = small
                    nsb[nsn] = big
                end
            end
        end
    end
    if numsignificant == 0
        print("Not significant at all")
    end
end

function check(s, b, nsn, nss, nsb)
    if nsn > 1
        for i in 1:nsn
            if (nss[i] <= s <= nsb[i]) && (nss[i] <= b <= nsb[i])
                return false
            end
        end
    end
    return true
end

mmulticomp([8, 11, 22, 6], [135.83, 160.49, 178.35, 188.06], [19.59, 12.28, 15.01, 9.81], method = "ryan")
# mean[  1] = 135.83000 vs mean[  4] = 188.06000 : diff =  52.23000, t = 6.53866 : P = 0.00000, alpha' = 0.00833
# mean[  1] = 135.83000 vs mean[  3] = 178.35000 : diff =  42.52000, t = 6.96308 : P = 0.00000, alpha' = 0.01250
# mean[  2] = 160.49000 vs mean[  4] = 188.06000 : diff =  27.57000, t = 3.67279 : P = 0.00066, alpha' = 0.01250
# mean[  1] = 135.83000 vs mean[  2] = 160.49000 : diff =  24.66000, t = 3.58814 : P = 0.00085, alpha' = 0.02500
# mean[  2] = 160.49000 vs mean[  3] = 178.35000 : diff =  17.86000, t = 3.26998 : P = 0.00212, alpha' = 0.02500

mmulticomp([8, 11, 22, 6], [135.83, 160.49, 178.35, 188.06], [19.59, 12.28, 15.01, 9.81], method = "tukey")
# mean[  1] = 135.83000 vs mean[  4] = 188.06000 : diff =  52.23000, WSD = 21.34697 : t = 9.24706 : P = 0.00000
# mean[  1] = 135.83000 vs mean[  3] = 178.35000 : diff =  42.52000, WSD = 15.57114 : t = 9.84728 : P = 0.00000
# mean[  2] = 160.49000 vs mean[  4] = 188.06000 : diff =  27.57000, WSD = 19.14118 : t = 5.19411 : P = 0.00258
# mean[  1] = 135.83000 vs mean[  2] = 160.49000 : diff =  24.66000, WSD = 16.11328 : t = 5.07440 : P = 0.00195
# mean[  2] = 160.49000 vs mean[  3] = 178.35000 : diff =  17.86000, WSD = 12.80554 : t = 4.62444 : P = 0.00482