Julia に翻訳--115 対応のあるデータの二つの相関係数の相等性の検定

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

対応のあるデータの二つの相関係数の相等性の検定
http://aoki2.si.gunma-u.ac.jp/R/diff_r2.html

ファイル名: diffr2.jl
関数名:    diffr2(データ行列)
          diffr2(相関係数行列, サンプルサイズ)

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

データ行列を与える場合と，cor と n を与える場合の 2 通りの関数を定義した。

==========#

using Statistics, Rmath

function diffr2(x)
    n, nc = size(x)
    nc == 4 && n != nc || error("required nx4 data matrix")
    diffr2(cor(x), n)
end

function diffr2(r, n)
    z12 = atanh(r[1,2])
    z34 = atanh(r[3,4])
    a1 = r[1, 3] * r[2, 4] + r[1, 4] * r[2, 3]
    a2 = -r[3, 4] * (r[1, 3] * r[2, 3] + r[1, 4] * r[2, 4])
    a3 = -r[1, 2] * (r[1, 3] * r[1, 4] + r[2, 3] * r[2, 4])
    a4 = r[1, 2] * r[3, 4] * (r[1, 3]^2 + r[1, 4]^2 + r[2, 3]^2 + r[2, 4]^2) / 2
    d = (1 - r[1, 2]^2) * (1 - r[3, 4]^2)
    chisq = (n - 3) * (z12 - z34)^2 / (2 - 2 * (a1 + a2 + a3 + a4) / d)
    df = 1
    p = pchisq(chisq, df, false)
    println("chisq. = $chisq,  df = $df,  p value = $p")
    Dict(:chisq => chisq, :df => df, :pvalue => p)
end

include("gendat.jl")
d = gendat(220, [0.439,0.288,0.329,0.354,0.320,0.595]);
cor(d)
# 4×4 Array{Float64,2}:
#  1.0    0.439  0.288  0.329
#  0.439  1.0    0.354  0.32
#  0.288  0.354  1.0    0.595
#  0.329  0.32   0.595  1.0
diffr2(d)
# chisq. = 5.500229788108368,  df = df,  p value = 0.01901397495593562

r = [  1.0    0.439  0.288  0.329
       0.439  1.0    0.354  0.32
       0.288  0.354  1.0    0.595
       0.329  0.32   0.595  1.0   ];
diffr2(r, 220)
# chisq. = 5.500229788108428,  df = 1,  p value = 0.019013974955934914