Julia に翻訳--140 Chow 検定

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

Chow 検定
http://aoki2.si.gunma-u.ac.jp/R/Chow.html

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

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

==========#

using Statistics, Rmath

chow = function(dat1, dat2)
    ess12 = ess(dat1) + ess(dat2)
    essc = ess(vcat(dat1, dat2))
    nr, df1 = size(dat1)
    df2 = nr + size(dat2, 1) - 2 * df1
    f = (essc - ess12) * df2 / (df1 * ess12)
    p = pf(f, df1, df2, false)
    println("F = $f,  df1 = $df1,  df2 = $df2,  p value = $p")
    Dict(:F => f, :df1 => df1, :df2 => df2, :pvalue => p)
end

function ess(dat) # 回帰分析を行い，誤差変動を返す関数
    r = cor(dat)
    betas = r[1:end-1, 1:end-1] \ r[1:end-1, end]
    SS = cov(dat, corrected=false) .* size(dat, 1)
    prop = [SS[i, i] for i in 1:size(SS, 1)]
    Sr = sum(SS[1:end-1, end] .* (betas .* sqrt.(prop[end] ./ prop[1:end-1])))
    St = SS[end, end]
    St - Sr
end

dat1 = [
    1.2 1.9 0.9
    1.6 2.7 1.3
    3.5 3.7 2.0
    4.0 3.1 1.8
    5.6 3.5 2.2
    5.7 7.5 3.5
    6.7 1.2 1.9
    7.5 3.7 2.7
    8.5 0.6 2.1
    9.7 5.1 3.6
    ];

dat2 = [
    1.4 1.3 0.5
    1.5 2.3 1.3
    3.1 3.2 2.5
    4.4 3.6 1.1
    5.1 3.1 2.8
    5.2 7.3 3.3
    6.5 1.5 1.3
    7.8 3.2 2.2
    8.1 0.1 2.8
    9.5 5.6 3.9
    ];

chow(dat1, dat2)
# F = 0.07154752544573709,  df1 = 3,  df2 = 14,  p value = 0.9742318542171906