Julia に翻訳--127 シェッフェ法，平均値の線形比較

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

シェッフェの方法による平均値の線形比較
http://aoki2.si.gunma-u.ac.jp/R/scheffe.html

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

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

かなりスイスイ書けるようになってきた。

==========#

using Rmath

function scheffe(n, m, u, g1, g2; conflevel = 0.95)
    all(length(n) == length(m) == length(u) &&
        all(n .> 1) &&
        all(u .> 0) &&
        all(floor.(n) .== n) &&
        all(floor.(g1) .== g1) &&
        all(floor.(g2) .== g2)) || error("parameter error")
    ng = length(n)
    k1 = ng - 1
    nc = sum(n)
    dfw = nc - ng
    Vw = sum(u .* (n .- 1)) / dfw
    n1 = length(g1)
    n2 = length(g2)
    g0 = collect(1:ng)
    g0 = [i for i in g0 if !(i in g1)]
    g0 = [i for i in g0 if !(i in g2)]
    n0 = ng - n1 - n2
    weight = vcat([repeat([i], j)
                   for (i, j) in zip([1 / n1, - 1 / n2, 0], [n1, n2, n0])]...)
    weight = weight[sortperm(vcat([g1, g2, g0]...))]
    theta = sum(weight .* m)
    Vtheta = Vw * sum(weight .^ 2 ./ n)
    confint = theta .- [1, - 1] *
        sqrt(k1 * qf(1 - conflevel, k1, dfw, false) * Vtheta)
    F0 = theta ^ 2 / k1 / Vtheta
    p = pf(F0, k1, dfw, false)
    println("theta = $theta,  V(theta) = $Vtheta")
    println("F = $F0,  df1 = $k1,  df2 = $dfw,  p value = $p")
    println("$(100conflevel) percent confidence interval = $confint")
    Dict(:theta => theta, :Vtheta => Vtheta, :F => F0,
        :df1 => k1, :df2 => dfw, :pvalue => p, :confint => confint)
end

n = [8, 11, 22, 6];
m = [135.83, 160.49, 178.35, 188.06];
sd = [19.59, 12.28, 15.01, 9.81];

scheffe(n, m, sd .^ 2, 1:2, 3:4)
# theta = -35.04499999999997,  V(theta) = 23.40938365266032
# F = 17.488030202230846,  df1 = 3,  df2 = 43,  p value = 1.435044718014149e-7
# 95.0 percent confidence interval = [-49.1218509474347, -20.968149052565245]

scheffe(n, m, sd .^ 2, 1:3, 4)
# theta = -29.836666666666673,  V(theta) = 42.81362201961475
# F = 6.9310236305159085,  df1 = 3,  df2 = 43,  p value = 0.000658278780876043
# 95.0 percent confidence interval = [-48.87379807609931, -10.799535257234034]