#==========
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
※コメント投稿者のブログIDはブログ作成者のみに通知されます