Julia に翻訳--141 Breslow-Day 検定

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

Breslow-Day 検定
http://aoki2.si.gunma-u.ac.jp/R/BD-test.html

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

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

@. なんか書かなくてもよいようにしてほしいなあ

==========#

using Rmath

function bdtest(m)
    d1, d2, k = size(m)
    Nk = vec(sum(m, dims=[1,2])) # vex() がないと，1x1x8 の3次元行列になる
    psiMH = sum(m[1, 1, :] .* m[2, 2, :] ./ Nk) / sum(m[2, 1, :] .* m[1, 2, :] ./ Nk)
    nk1 = m[1, 1, :] .+ m[1, 2, :]
    nk2 = m[2, 1, :] .+ m[2, 2, :]
    xkt = m[1, 1, :] .+ m[2, 1, :]
    a1 = psiMH - 1
    @. b1 = -psiMH * (nk1 + xkt) - nk2 + xkt
    @. c1 = psiMH * nk1 * xkt
    @. e = (-b1 - sqrt(b1 ^ 2 - 4 * a1 * c1)) / (2 * a1)
    @. v = 1 / (1 / e + 1 / (nk1 - e) + 1 / (xkt - e) + 1 / (nk2 - xkt + e))
    chisqBD = sum((m[1, 1, :] .- e) .^ 2 ./ v) # この行は @. を使うとエラーになる
    df = k - 1
    p = pchisq(chisqBD, df, false)
    println("chisq. = $chisqBD,  df = $df,  p value = $p")
    Dict(:chisq => chisqBD, :df => df, :pvalue => p)
end

m = [1, 8, 0, 10,
     2, 47, 0, 60,
     3, 29, 1, 13,
     3, 17, 0, 7,
     13, 45, 1, 45,
     13, 26, 0, 18,
     8, 15, 0, 10,
     11, 4, 0, 4 ];
m = reshape(m, 2, 2, :)

bdtest(m)
# chisq. = 8.943202421541203,  df = 7,  p value = 0.25675965693687036