Julia に翻訳--143 クロンバックのα信頼性係数

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

クロンバックのα信頼性係数
http://aoki2.si.gunma-u.ac.jp/R/alpha.html

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

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

出力にはこだわらない。
必要なら，戻り値を使って。

==========#

using Statistics

function alpha(x; detail = false)
    nr, k = size(x)
    (k > 1 && nr > 1) || error("parameter error")
    α = alpha0(x)
    println("α = $α")
    if detail == true && k >= 3
        R2 = 1 .- 1 ./ diag(inv(cor(x)))
        α2 = similar(R2)
        cor2 = similar(R2)
        for i in 1:k
            x2 = x[:, [i != j for j in 1:k]]
            α2[i] = alpha0(x2)
            cor2[i] = cor(x[:, i], vec(sum(x2, dims=2)))
        end
        println("それぞれの変数を除いたときの α\n$α2")
        println("それぞれの変数とその変数を除いたときの合計値との相関係数\n$cor2")
        println("それぞれの変数をその変数以外の変数で予測したときの決定係数\n$R2")
        Dict(:alpha => α, :alpha2 => α2, :R2 => R2, :cor2 => cor2)
    end
end

function alpha0(x)
    k = size(x, 2)
    VarCovMat = cov(x)
    Sy2 = sum(VarCovMat)
    Sj2 = sum(diag(VarCovMat))
    k / (k - 1) * (1 - Sj2 / Sy2)
end

diag(a) = [a[i, i] for i in 1:size(a, 1)]

x = [
    49 44 37 54
    36 36 36 42
    42 51 45 35
    47 67 54 40
    54 59 68 54
    51 55 59 67
    45 36 48 46
    72 49 50 58];

alpha(x)
alpha(x, detail = true)
# α = 0.7403365794138637
# それぞれの変数を除いたときの α
# [0.6657399977072106, 0.7307558528428093, 0.592778807780536, 0.72184]
# それぞれの変数とその変数を除いたときの合計値との相関係数
# [0.5599374682949105, 0.4444039960775981, 0.6801937600489253, 0.4589177292647371]
# それぞれの変数をその変数以外の変数で予測したときの決定係数
# [0.4298971517659005, 0.551857280000196, 0.6071655584759322, 0.4821633678144336]