Julia に翻訳--120 比率の多重比較，ライアン法，テューキー法

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

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

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

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

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, Printf

function pmulticomp(n, r; alpha = 0.05, method = "ryan") # or method = "tukey"
    k = length(n)
    all(k == length(r) && all(n .> 0) && all(r .>= 0) && all(r .<= n) &&
        all(floor.(n) .== n) && all(floor.(r) .== r)) || error("parameter error")
    o = sortperm(r ./ n)
    sr = r[o]
    sn = n[o]
    srsn = sr ./ sn
    numsignificant = nsn = 0
    nss = zeros(Int, k * (k - 1) ÷ 2)
    nsb = zeros(Int, k * (k - 1) ÷ 2)
    # m = k; small = 1
    q = NaN # dummy これがないと動かない
    for m = k:-1:2
        for small = 1:k - m + 1
            big = small + m - 1
            if check(small, big, nsn, nss, nsb)
                prop = sum(sr[small:big]) / sum(sn[small:big])
                se = sqrt(prop * (1 - prop) * (1 / sn[small] + 1 / sn[big]))
                diff = srsn[big] - srsn[small]
                if method == "ryan"
                    nominalalpha = 2 * alpha / (k * (m - 1))
                    rd = se * qnorm(nominalalpha / 2, false)
                    z = se == 0 ? 0 : diff / se
                    p = pnorm(z, false) * 2
                    result = p <= nominalalpha
                else
                    if m == k
                        q = qtukey(1-alpha, k, Inf)
                        qq = q
                    else
                        qq = (q + qtukey(1-alpha, m, Inf)) / 2
                    end
                    WSD = se * qq / sqrt(2)
                    result = diff != 0 && diff >= WSD
                end
                if result
                    numsignificant = 1
                    @printf("p[%2i] = %7.5f vs p[%2i] = %7.5f : diff = %7.5f, ",
                        o[small], srsn[small], o[big], srsn[big], diff)
                    if method == "ryan"
                        @printf("RD = %7.5f : P = %7.5f, alpha' = %7.5f\n", rd, p, nominalalpha)
                    else
                        @printf("WSD = %7.5f\n", WSD)
                    end
                else
                    nsn += 1
                    nss[nsn] = small
                    nsb[nsn] = big
                end
            end
        end
    end
    if numsignificant == 0
        println("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

n = [30, 35, 47, 21, 45];
r = [2, 4, 14, 13, 39];

pmulticomp(n, r)
# p[ 1] = 0.06667 vs p[ 5] = 0.86667 : diff = 0.80000, RD = 0.32472 : P = 0.00000, alpha' = 0.00500
# p[ 1] = 0.06667 vs p[ 4] = 0.61905 : diff = 0.55238, RD = 0.33341 : P = 0.00001, alpha' = 0.00667
# p[ 2] = 0.11429 vs p[ 5] = 0.86667 : diff = 0.75238, RD = 0.30528 : P = 0.00000, alpha' = 0.00667
# p[ 1] = 0.06667 vs p[ 3] = 0.29787 : diff = 0.23121, RD = 0.23054 : P = 0.00979, alpha' = 0.01000
# p[ 2] = 0.11429 vs p[ 4] = 0.61905 : diff = 0.50476, RD = 0.32612 : P = 0.00007, alpha' = 0.01000
# p[ 3] = 0.29787 vs p[ 5] = 0.86667 : diff = 0.56879, RD = 0.26479 : P = 0.00000, alpha' = 0.01000
# p[ 3] = 0.29787 vs p[ 4] = 0.61905 : diff = 0.32118, RD = 0.29877 : P = 0.01239, alpha' = 0.02000

pmulticomp(n, r, method="tukey")
# p[ 1] = 0.06667 vs p[ 5] = 0.86667 : diff = 0.80000, WSD = 0.31555
# p[ 1] = 0.06667 vs p[ 4] = 0.61905 : diff = 0.55238, WSD = 0.32546
# p[ 2] = 0.11429 vs p[ 5] = 0.86667 : diff = 0.75238, WSD = 0.29800
# p[ 1] = 0.06667 vs p[ 3] = 0.29787 : diff = 0.23121, WSD = 0.22695
# p[ 2] = 0.11429 vs p[ 4] = 0.61905 : diff = 0.50476, WSD = 0.32104
# p[ 3] = 0.29787 vs p[ 5] = 0.86667 : diff = 0.56879, WSD = 0.26067
# p[ 3] = 0.29787 vs p[ 4] = 0.61905 : diff = 0.32118, WSD = 0.30102