Julia に翻訳--123 テューキー法，平均値の対比較

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

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

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

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

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

二重配列（リスト）を平坦化するには，vcat(A...) を使う。

==========#

using Statistics, Combinatorics, Rmath, Printf

function tukey(data, group; method = "Tukey") # or method = "Games-Howell"
    indices, n = table(group)
    a = length(n)
    npairs = binomial(a, 2)
    phie = sum(n) - a
    Mean = [mean(data[group .== g]) for g in indices]
    Variance = [var(data[group .== g]) for g in indices]
    pair = combinations(1:a, 2)
    if method == "Tukey"
        ve = sum((n .- 1) .* Variance) / phie
        t = [abs(Mean[i] - Mean[j]) / sqrt(ve * (1 / n[i] + 1 / n[j])) for (i, j) in pair]
        p = 1 .- ptukey.(t .* sqrt(2), a, phie)#, false)
        @printf(" pair  %7s  %7s\n", "t value", "p value")
        for (i, (j, k)) in enumerate(pair)
            @printf("%2d:%2d  %7.5f  %7.5f\n", j, k, t[i], p[i])
        end
        Dict(:t => t, :pvalue => p, :phi => phie, :v => ve)
    else
        t = [ts0(Mean[i], Variance[i], n[i], Mean[j], Variance[j], n[j]) for (i, j) in pair]
        df = [dfs0(Variance[i], n[i], Variance[j], n[j]) for (i, j) in pair]
        p = 1 .- ptukey.(t .* sqrt(2), a, df) #, false)
        @printf(" pair  %7s  %9s  %7s\n", "t value", "df", "p value")
        for (i, (j, k)) in enumerate(pair)
            @printf("%2d:%2d  %7.5f  %9.5f  %7.5f\n", j, k, t[i], df[i], p[i])
        end
        Dict(:t => t, :df => df, :pvalue => p)
    end
end

ts0(mi, vi, ni, mj, vj, nj) = abs(mi - mj) / sqrt(vi / ni + vj / nj)

dfs0(vi, ni, vj, nj) = (vi / ni + vj / nj) ^ 2 /
                       ((vi / ni)^2 / (ni - 1) + (vj / nj)^2 / (nj - 1))

function table(x) # indices が少ないとき
    indices = sort(unique(x))
    counts = zeros(Int, length(indices))
    for i in indexin(x, indices)
        counts[i] += 1
    end
    return indices, counts
end

data = [
  14., 15, 14, 16, 15, 17, 17,
  17, 16, 17, 16, 15, 18, 19, 15,
  18, 19, 20, 19, 17, 17,
  20, 21, 19, 20, 19, 22, 20,
  19, 20, 19, 17, 17, 17, 18
];

group = [repeat([i], count) for (i, count) in enumerate([7, 8, 6, 7, 7])]
group = vcat(group...)

tukey(data, group)
# pair  t value  p value
# 1: 2  1.85409  0.36302
# 1: 3  4.18754  0.00198
# 1: 4  7.07368  0.00000
# 1: 5  4.07273  0.00269
# 2: 3  2.53703  0.10906
# 2: 4  5.45158  0.00006
# 2: 5  2.35220  0.15676
# 3: 4  2.60863  0.09414
# 3: 5  0.27459  0.99868
# 4: 5  3.00096  0.03978

tukey(data, group, method = "Games-Howell")
# pair  t value         df  p value
# 1: 2  1.72859   12.97639  0.45125
# 1: 3  4.21141   10.84628  0.01043
# 1: 4  7.50517   11.65358  0.00007
# 1: 5  4.08185   11.97449  0.01084
# 2: 3  2.43499   11.69245  0.17285
# 2: 4  5.48706   12.78706  0.00087
# 2: 5  2.24124   12.99984  0.22492
# 3: 4  2.83393   10.14001  0.10023
# 3: 5  0.28228   10.70743  0.99838
# 4: 5  3.26970   11.80872  0.04412