Julia に翻訳--139 分散・共分散行列の同等性の検定

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

分散・共分散行列の同等性の検定
http://aoki2.si.gunma-u.ac.jp/R/eq-cov.html

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

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

関数を入れ子にして一行で書くのもよいが，for ループで書くとわかりやすいというメリットもある。

==========#

using Statistics, LinearAlgebra, Rmath

function boxm(x, gvar)
    nv = size(x, 2)
    nv >= 2 || error("分散共分散行列を計算する変数は2個以上必要です")
    index, ni = table(gvar)
    n = length(gvar)
    g = length(ni)
    sumlogdetSi = 0
    S = zeros(nv, nv)
    for (i, gvari) in enumerate(index)
        y = x[gvar .== gvari, :]
        covy = cov(y)
        sumlogdetSi += (ni[i] - 1) * log(det(covy))
        S .+= (ni[i] - 1) .* covy
    end
    S ./= n - g
    M = (n - g) * log(det(S)) - sumlogdetSi
    df1 = (g - 1) * nv * (nv + 1) / 2
    rho = 1 - (2 * nv ^ 2 + 3 * nv - 1) / (6 * (nv + 1) * (g - 1)) * (sum(1 ./ (ni .- 1)) - 1 / (n - g))
    tau = (nv - 1) * (nv + 2) / (6 * (g - 1)) * (sum(1 ./ (ni .- 1) .^ 2) - 1 / (n - g) ^ 2)
    df2 = (df1 + 2) / abs(tau - (1 - rho) ^ 2)
    gamma = (rho - df1 / df2) / df1
    F = M * gamma
    p = pf(F, df1, df2, false)
    println("M = $M\nF = $F,  df1 = $df1,  df2 = $df2\np value = $p")
    Dict(:M => M, :F => F, :df1 => df1, :df2 => df2, :pvalue => p)
end

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

function rep(x, n::Array{Int64,1})
    length(x) == length(n) || error("length(x) wasn't length(n)")
    vcat([repeat([x[i]], n[i]) for i = 1:length(x)]...)
end

x = [
   2.9 161.7 120.8
   2.3 114.8  85.2
   2.0 128.4  92.0
   3.2 149.2  97.3
   2.7 126.0  81.1
   4.4 133.8 107.6
   4.1 161.3 114.0
   2.1 111.5  77.3
   4.8 198.7 172.9
   3.6 199.3 157.9
   2.0 188.4 152.7
   4.9 183.6 164.2
   3.9 173.5 172.2
   4.4 184.9 163.2 ];
gvar = rep(1:2, [8, 6]);

boxm(x, gvar)
# M = 9.730421574785467
# F = 1.1535477771673381,  df1 = 6.0,  df2 = 795.1946969489453
# p value = 0.32937843690375435