Julia に翻訳--106 二要因の分散分析，SAB タイプ，RBFpq デザイン，被検者内計画

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

二要因の分散分析（SAB タイプ；RBFpq デザイン；被検者内計画）
http://aoki2.si.gunma-u.ac.jp/R/SAB.html

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

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

R の apply(X, MARGIN, FUN) と Julia の FUN(X, dims) の関係
reshape は不要な場合もある

x = reshape(1:24, (2, 3, 4))            # == x = array(1:24, dim=c(2,3,4))
reshape(mean(x, dims=1), 3, 4)          # == apply(x, 2:3, mean)
reshape(mean(x, dims=2), 2, 4)          # == apply(x, c(1, 3), mean)
reshape(mean(x, dims=3), 2, 3)          # == apply(x, 1:2, mean)
reshape(mean(x, dims=[1,2]), 4, 1, 1)   # == apply(x, 3, mean)
reshape(mean(x, dims=[2,3]), 2, 1, 1)   # == apply(x, 1, mean)
reshape(mean(x, dims=[1,3]), 3, 1, 1)   # == apply(x, 2, mean)
reshape(mean(x, dims=[1,2,3]), 1, 1, 1) # == mean(x)

==========#

using Rmath, Statistics, Printf

function sab(data)
    Nb, Na, N = size(data)
    Xbar = vec(mean(data, dims=3))
    SD = vec(std(data, dims=3)) .* sqrt((N - 1) / N)
    v1 = mean(data)
    v2 = sum((mean(data, dims=[1,3]) .- v1) .^ 2) * Nb * N
    v3 = sum((mean(data, dims=[2,3]) .- v1) .^ 2) .* Na .* N
    v4 = sum((mean(data, dims=[3]) .- v1) .^ 2) * N
    v42 = v4-v2-v3
    v5 = sum(SD .^ 2) * N
    v6 = sum(data) .^ 2 / (N * Na * Nb)
    v62 = sum(sum(data, dims=[1,2]) .^ 2) / (Na * Nb) - v6
    v7 = sum(sum(data, dims=1) .^ 2) / Nb - v6 - v62 - v2
    v8 = sum(sum(data, dims=2) .^ 2) / Na - v6 - v62 - v3
    v9 = v5 - v62 - v7 - v8
    SS = [v62, v2, v7, v3, v8, v42, v9]
    dfs = N - 1
    dfa = Na - 1
    dfb = Nb - 1
    df = [dfs, dfa, dfs * dfa, dfb, dfs * dfb, dfa * dfb, dfs * dfa * dfb]
    MS = SS ./ df
    F = fill(NaN, 7)
    P = fill(NaN, 7)
    F[[2, 4, 6]] = MS[[2, 4, 6]] ./ MS[[3, 5, 7]]
    P[[2, 4, 6]] = pf.(F[[2, 4, 6]], df[[2, 4, 6]], df[[3, 5, 7]], false)
    @printf("%5s %11s  %3s  %11s  %9s  %7s\n",
            "", "SS", "df", "MS", "F value", "P value")
    name = ["S", "A", "SxA", "B", "SxB", "AxB", "SxAxB"]
    for i = 1:7
        if i % 2 == 1
            @printf("%5s %11.8g  %3d  %11.8g\n", name[i], SS[i], df[i], MS[i])
        else
            @printf("%5s %11.8g  %3d  %11.8g  %9.5f  %7.5f\n",
                    name[i], SS[i], df[i], MS[i], F[i], P[i])
        end
    end
    Dict(:name => name, :SS => SS, :df => df, :MS => MS, :F => F, :P => P)
end

data = [4, 3, 6, 3, 5, 5, 4, 4, 6, 4, 6, 4, 5, 4, 4, 3, 7, 4];
data = reshape(data, 3, 2, 3)
sab(data)
#                SS   df           MS    F value  P value
#     S  0.33333333    2   0.16666667
#     A 0.055555556    1  0.055555556    1.00000  0.42265
#   SxA  0.11111111    2  0.055555556
#     B           4    2            2    1.71429  0.28994
#   SxB   4.6666667    4    1.1666667
#   AxB   11.111111    2    5.5555556   10.00000  0.02778
# SxAxB   2.2222222    4   0.55555556

Julia の小ネタ--012 R の apply(X, MARGIN, FUN) と Julia の FUN(X, dims) の関係

R の apply(X, MARGIN, FUN) と Julia の FUN(X, dims) の関係
reshape は不要な場合もある

x = reshape(1:24, (2, 3, 4))            # == x = array(1:24, dim=c(2,3,4))
reshape(mean(x, dims=1), 3, 4)          # == apply(x, 2:3, mean)
reshape(mean(x, dims=2), 2, 4)          # == apply(x, c(1, 3), mean)
reshape(mean(x, dims=3), 2, 3)          # == apply(x, 1:2, mean)
reshape(mean(x, dims=[1,2]), 4, 1, 1)   # == apply(x, 3, mean)
reshape(mean(x, dims=[2,3]), 2, 1, 1)   # == apply(x, 1, mean)
reshape(mean(x, dims=[1,3]), 3, 1, 1)   # == apply(x, 2, mean)
reshape(mean(x, dims=[1,2,3]), 1, 1, 1) # == mean(x)

Julia に翻訳--105 二元配置分散分析

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

二元配置分散分析
http://aoki2.si.gunma-u.ac.jp/R/twoway-anova.html

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

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

元のプログラムにある二次元配列 mab は不要だ。

==========#

using Statistics, Rmath, Printf

function twowayanova(x, a, b)
    indexa, indexb, tbl = table(a, b)
    n = length(x)
    dft = n-1
    na, nb = size(tbl)
    dfa = na - 1
    dfb = nb - 1
    gm = mean(x)
    MSt = var(x)
    sst = dft * MSt
    ma = [mean(x[a .== i]) for i in indexa]
    mb = [mean(x[b .== j]) for j in indexb]
    if any(tbl .== 0)
        error("いくつかの水準の組合せにおいて，データがないセルがあります")
    elseif all(tbl .== 1)
        ssa = nb * sum((ma .- gm) .^ 2)
        ssb = na * sum((mb .- gm) .^ 2)
        sse = sst-ssa-ssb
        dfe = dft-dfa-dfb
        ss = [ssa, ssb, sse, sst]
        df = [dfa, dfb, dfe, dft]
        ms = ss ./ df
        f = ms ./ ms[3]
        p = pf.(f, df, df[3], false)
        f[3:4] .= NaN
        p[3:4] .= NaN
        result = hcat(ss, df, ms, f, p)
        @printf("%12s  %6s  %10s  %10s  %10s\n", "SS", "df", "MS", "F value", "P value")
        @printf("a %10.7g  %6d  %10.7g  %10.7f  %10.7f\n", ssa, dfa, ms[1], f[1], p[1])
        @printf("b %10.7g  %6d  %10.7g  %10.7f  %10.7f\n", ssb, dfb, ms[2], f[2], p[2])
        @printf("e %10.7g  %6d  %10.7g\n", sse, dfe, ms[3])
        @printf("T %10.7g  %6d  %10.7g\n", sst, dft, ms[4])
        Dict(:result => result)
    elseif all(tbl .>= 2) && hirei(tbl)
        sse = 0
        for i in 1:na, j in 1:nb
            selected = x[(a .== indexa[i]) .& (b .== indexb[j])]
            sse += sum((selected .- mean(selected)) .^ 2)
        end
        ssa = sum(sum(tbl, dims=2) .* (ma .- gm) .^ 2)
        ssb = sum(transpose(sum(tbl, dims=1)) .* (mb .- gm) .^ 2)
        ss = [ssa, ssb, sst-ssa-ssb-sse, sse, sst]
        dfab = dfa * dfb
        df = [dfa, dfb, dfab, dft-dfa-dfb-dfab, dft]
        ms = ss ./ df
        f = ms ./ ms[4]
        p = pf.(f, df, df[4], false)
        f[4:5] .= NaN
        p[4:5] .= NaN
        result = hcat(ss, df, ms, f, p)
        names = ["a", "b", "a*b", "e", "T"]
        output("Model-I", names, ss, df, ms, f, p)
        f[1:2] = ms[1:2] ./ ms[3]
        p[1:2] = pf.(f[1:2], df[1:2], df[3], false)
        result2 = hcat(ss, df, ms, f, p)
        println()
        output("Model-II", names, ss, df, ms, f, p)
        Dict(:model1 => result, :model2 => result2)
    else
        d1 = makedummy(a, indexa);
        d2 = makedummy(b, indexb);
        d3 = crossdummy(d1, d2);
        rabab = mreg(hcat(d3, x))
        rab = mreg(hcat(d1, d2, x))
        ra = mreg(hcat(d1, x))
        rb = mreg(hcat(d2, x))
        ss = [rab[1]-rb[1], rab[1]-ra[1], rabab[1]-rab[1], rabab[2], rabab[3]]
        df = [dfa, dfb, dfa*dfb, n-na*nb, n-1]
        ms = ss ./ df
        f = ms ./ ms[4]
        p = pf.(f, df, df[4], false)
        f[4:5] .= NaN
        p[4:5] .= NaN
        result = hcat(ss, df, ms, f, p)
        names = ["a", "b", "a*b", "e", "T"]
        output(names, ss, df, ms, f, p)
        Dict(:result => result)
    end
end

function hirei(tbl)
    for i in 2:size(tbl, 1)
        temp = tbl[i, :] ./ tbl[1, :]
        all(temp .== temp[1]) || return false
    end
    return true
end

function output(model, names, ss, df, ms, f, p)
    println("$model")
    @printf("%14s  %6s  %10s  %10s  %10s\n", "SS", "df", "MS", "F value", "P value")
    for i in 1:3
        @printf("%4s%10.7g  %6d  %10.7g  %10.7f  %10.7f\n",
                names[i], ss[i], df[i], ms[i], f[i], p[i])
    end
    for i in 4:5
        @printf("%4s%10.7g  %6d  %10.7g\n", names[i], ss[i], df[i], ms[i])
    end
end

function mreg(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[1:end-1] ./ prop[end])))
    St = SS[end, end]
    (Sr, St - Sr, St)  # 回帰変動，誤差変動，全変動
end

function makedummy(x, y)
    n = length(x)
    k = length(y)
    d = zeros(n, k)
    for (i, j) = enumerate(indexin(x, y))
        d[i, j] = 1
    end
    d[:, 2:end]
end

function crossdummy(a, b)
    n1, k = size(a)
    n2, l = size(b)
    n1 == n2 || error("data error")
    d = zeros(n1, k * l)
    for i = 1:l
        for j = 1:k
            d[:, j + (i - 1) * k] = a[:, j] .* b[:, i]
        end
    end
    hcat(a, b, d)
end

function table(x, y)
    indicesx = sort(unique(x))
    indicesy = sort(unique(y))
    counts = zeros(Int, length(indicesx), length(indicesy))
    for (i, j) in zip(indexin(x, indicesx), indexin(y, indicesy))
        counts[i, j] += 1
    end
    return indicesx, indicesy, counts
end

a = [1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3];
b = repeat(1:4, 3);
x = [9, 17, 12, 16, 1, 21, 16, 11, 7, 19, 6, 9];
twowayanova(x, a, b)
#           SS      df          MS     F value     P value
# a       21.5       2       10.75   0.6592845   0.5510300
# b   268.6667       3    89.55556   5.4923339   0.0371924
# e   97.83333       6    16.30556
# T        388      11    35.27273

x = [24.8, 23.9, 24.1, 28.8, 22.6, 28.0, 26.4, 27.4, 29.4, 30.0, 28.7, 
    29.2, 25.0, 26.6, 27.9, 28.5, 27.1, 25.2, 27.9, 29.2, 26.7, 25.0, 
    29.9, 29.4, 27.5, 32.5, 29.5, 26.3, 28.2, 31.8, 30.2, 31.7, 29.2, 
    30.3, 30.9, 28.4, 29.8, 26.7, 30.7, 28.5, 26.5, 26.7, 31.7, 27.2, 
    25.5, 29.9, 27.3, 28.8, 28.5, 25.7, 28.7, 31.3, 29.4, 29.8, 30.3, 
    29.6, 31.7, 33.6, 32.0, 34.3];
a = [1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 1, 1, 1, 2, 2, 2, 3, 3, 
    3, 4, 4, 4, 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 1, 1, 1, 2, 2, 
    2, 3, 3, 3, 4, 4, 4, 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4];
b = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 
    2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 
    4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5];
twowayanova(x, a, b)
#  Model-I
#            SS      df          MS     F value     P value
#   a  51.39733       3    17.13244   4.9370667   0.0052039
#   b  106.2873       4    26.57183   7.6572211   0.0001115
# a*b  52.85267      12    4.404389   1.2692154   0.2738926
#   e  138.8067      40    3.470167
#   T   349.344      59    5.921085
#
#  Model-II
#            SS      df          MS     F value     P value
#   a  51.39733       3    17.13244   3.8898573   0.0373908
#   b  106.2873       4    26.57183   6.0330352   0.0067194
# a*b  52.85267      12    4.404389   1.2692154   0.2738926
#   e  138.8067      40    3.470167
#   T   349.344      59    5.921085

a = [1, 1, 1, 1, 2, 2, 2, 2, 2, 2]
b = [1, 1, 2, 2, 1, 1, 2, 2, 2, 2]
x = [17, 16, 25, 22, 18, 26, 34, 30, 34, 30]
twowayanova(x, a, b)
#            SS      df          MS     F value     P value
# a    121.4405       1    121.4405  13.7479784   0.0099953
# b    177.1905       1    177.1905  20.0592992   0.0041979
# a*b  5.142857       1    5.142857   0.5822102   0.4743687
# e          53       6    8.833333
# T       415.6       9    46.17778

Julia に翻訳--104 クラスカル・ウォリス検定と多重比較

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

クラスカル・ウォリス検定（plus 多重比較）
http://aoki2.si.gunma-u.ac.jp/R/kruskal-wallis.html

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

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

==========#

using StatsBase, Rmath, Combinatorics, Printf

function kruskalwallis(data, group)
    idx, ni = table(group)
    n = sum(ni)
    r = tiedrank(data)
    Ri = [sum(r[group .== i]) for i in idx]
    chisq = 12 * sum(Ri .^ 2 ./ ni) / (n * (n + 1)) - 3 * (n + 1)
    if length(unique(data)) != n
        idx2, tie = table(data)
        chisq /= 1 - sum(tie .^ 3 - tie) / (n ^ 3 - n)
    end
    a = length(ni)
    df = a - 1
    pvalue = pchisq(chisq, df, false)
    println("chisq. = $chisq,  df = $df,  p value = $pvalue")
    Dict()
    amean = (n + 1) / 2
    Rmean = Ri ./ ni
    V = sum((r .- amean) .^ 2) / (n - 1)
    pair = combinations(1:a, 2)
    S = [(Rmean[i]- Rmean[j]) ^ 2 / (V / ni[i] + V / ni[j]) for (i, j) in pair]
    P = pchisq.(S, df, false)
    pair2 = [(i, j) for (i, j) in pair]
    @printf("%2s: %2s  %7s  %9s\n", "i", "j", "Chi sq.", "p value")
    for i in 1:a
        @printf("%2d: %2d  %7.3f  %9.7f\n", pair2[i][1], pair2[i][2], S[i], P[i])
    end
    Dict(:chisq => chisq, :df => df, :pvalue => pvalue, :S => S, :P => P)
end

function table(x)
    indices = sort(unique(x))
    counts = zeros(Int, length(indices))
    for i = 1:length(indices)
        counts[i] = count(isequal(indices[i]), x)
    end
    return indices, counts
end

data = [ 3.42, 3.84, 3.96, 3.76,
         3.17, 3.63, 3.47, 3.44, 3.39,
         3.64, 3.72, 3.91 ];

group = [1, 1, 1, 1,  2, 2, 2, 2, 2,  3, 3, 3];

kruskalwallis(data, group)
# chisq. = 5.54871794871795,  df = 2,  p value = 0.06238945711580609
#  i:  j  Chi sq.    p value
#  1:  2    4.104  0.1284601
#  1:  3    0.004  0.9981702
#  2:  3    3.703  0.1570357

Julia に翻訳--103 分割表データからクラスカル・ウォリス検定

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

クラスカル・ウォリス検定（分割表データから）
http://aoki2.si.gunma-u.ac.jp/R/kw3.html

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

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

sum(tbl, dims= 1) と sum(tbl, dims=2) が，列ベクトルと行ベクトルの区別があるのは，いいようは悪いような。

==========#

using Rmath

function kw4test(tbl)
    nr, nc = size(tbl)
    ni = sum(tbl, dims=2)
    n = sum(tbl)
    tj = vec(sum(tbl, dims=1))
    rj = vcat(0, cumsum(tj[1:end - 1])) .+ (tj .+ 1) ./ 2
    Sx = 12 * sum((tbl * rj) .^ 2 ./ ni) / (n * (n + 1)) - 3 * (n + 1)
    S0 = Sx / (1 - sum(tj .^ 3 .- tj) / (n ^ 3 - n))
    df = nr - 1
    pvalue = pchisq(S0, df, false)
    println("chisq. = $S0,  df = $df,  p value = $pvalue")
    Dict(:chisq => S0, :df => df, :pvalue => pvalue)
end

tbl = [10 5 1
      4 7 3
      2 4 9];

kw4test(d)
# chisq. = 12.417346395306067,  df = 2,  p value = 0.0020119050943707304

Julia に翻訳--102 フリードマン検定と多重比較

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

フリードマン検定（plus 多重比較）
http://aoki2.si.gunma-u.ac.jp/R/friedman.html

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

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

内包表記より，for ループの方がわかりやすいかな。

==========#

using StatsBase, Rmath, Combinatorics, Printf

function friedman(dat)
    row, col = size(dat)
    df = col - 1
    o = similar(dat, Float64)
    ties = 0
    for i = 1:row
        rank = tiedrank(dat[i, :])
        tie = [count(isequal(s), rank) for s in Set(rank)]
        ties += sum(tie .^ 3 .- tie)
        o[i, :] = rank
    end
    R = sum(o, dims=1)
    chi = 12 * sum((R .- row * (col + 1) / 2) .^ 2) / (row * col * (col + 1) - ties / (col - 1))
    pvalue = pchisq(chi, df, false)
    println("chisq. = $chi,  df = $df,  p value = $pvalue")
    Rm = R / row
    V = sum((o .- (col + 1) / 2) .^ 2)
    S = []
    P = []
    @printf("%2s: %2s  %7s  %9s\n", "i", "j", "Chi sq.", "p value")
    for (i, j) in combinations(1:col, 2)
        s = row^2*df*(Rm[i] - Rm[j])^2/(2*V)
        append!(S, s)
        p = pchisq(s, df, false)
        append!(P, p)
        @printf("%2d: %2d  %7.3f  %9.7f\n", i, j, s, p)
    end
   Dict(:chisq => chi, :df => df, :pvalue => pvalue, :S => S, :P => P)
end

dat = [  5 60 35 62 76
        24 44 74 63 76
        56 57 70 74 79
        44 51 55 23 84
         8 68 50 24 64
        32 66 45 63 46
        25 38 70 58 77
        48 24 40 80 72];

friedman(dat)
# chisq. = 15.9,  df = 4,  p value = 0.0031563263734228895
#  i:  j  Chi sq.    p value
#  1:  2    3.600  0.4628369
#  1:  3    4.225  0.3764110
#  1:  4    5.625  0.2289584
#  1:  5   15.625  0.0035659
#  2:  3    0.025  0.9999225
#  2:  4    0.225  0.9941270
#  2:  5    4.225  0.3764110
#  3:  4    0.100  0.9987909
#  3:  5    3.600  0.4628369
#  4:  5    2.500  0.6446358

x = [9 17 12 16; 5 21 16 11; 7 19 6 9; 8 11 11 8; 9 8 9 9; 2 4 5 8; 3 8 10 9];
friedman(x)
# chisq. = 6.1875,  df = 3,  p value = 0.10283586999577443
#  i:  j  Chi sq.    p value
#  1:  2    4.687  0.1961632
#  1:  3    3.797  0.2842498
#  1:  4    3.797  0.2842498
#  2:  3    0.047  0.9973385
#  2:  4    0.047  0.9973385
#  3:  4    0.000  1.0000000

x = [9 17 12 16; 1 21 16 11; 7 19 6 9];
friedman(x)
# chisq. = 7.0,  df = 3,  p value = 0.07189777249646513
#  i:  j  Chi sq.    p value
#  1:  2    6.400  0.0936908
#  1:  3    0.400  0.9402425
#  1:  4    1.600  0.6593898
#  2:  3    3.600  0.3080222
#  2:  4    1.600  0.6593898
#  3:  4    0.400  0.9402425

Julia に翻訳--101 乱塊法

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

乱塊法
http://aoki2.si.gunma-u.ac.jp/R/randblk.html

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

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

%w.dG の結果が他の言語とちょっと違う（末尾0が除かれる）のが，ちょっtいや。

==========#

using Statistics, Rmath, Printf

function randblk(dat)
    nr, nc = size(dat)
    gmean = mean(dat)
    tmean = mean(dat, dims=1)
    rmean = mean(dat, dims=2)
    SSt = sum((dat .- gmean) .^ 2)
    SSc = nr * sum((tmean .- gmean) .^ 2)
    SSr = nc * sum((rmean .- gmean) .^ 2)
    SSe = SSt - SSr - SSc
    SS = [SSc, SSr, SSe, SSt]
    dfc = nc - 1
    dfr = nr - 1
    dfe = dfc * dfr
    dft = nc * nr - 1
    df = [dfc, dfr, dfe, dft]
    MS =  SS ./ df
    Fs = MS / MS[3]
    Ps = pf.(Fs, df, dfe, false)
    Fs[3:4] .= NaN
    Ps[3:4] .= NaN
    result = hcat(SS, df, MS, Fs, Ps)
    @printf("%11s  %10s  %4s  %10s  %10s  %10s\n", "Source", "SS", "d.f.", "MS", "F value", "P value")
    @printf("%11s  %10.7g  %4s  %10.7g  %10.7g  %10.5f\n", "Treatment", SSc, dfc, MS[1], Fs[1], Ps[1])
    @printf("%11s  %10.7g  %4s  %10.7g  %10.7g  %10.5f\n", "Replication", SSr, dfr, MS[2], Fs[2], Ps[2])
    @printf("%11s  %10.7g  %4s  %10.7g\n", "Residual", SSe, dfe, MS[3])
    @printf("%11s  %10.7g  %4s  %10.7g\n", "Total", SSt, dft, MS[4])
    Dict(:SS => SS, :df => df, :MS => MS, :Fs => Fs[1:2], :Ps => Ps[1:2])
end

dat = [9  17  12  16 
       1  21  16  11 
       7  19   6   9]

randblk(dat)
#      Source          SS  d.f.          MS     F value     P value
#   Treatment    268.6667     3    89.55556    5.492334     0.03719
# Replication        21.5     2       10.75   0.6592845     0.55103
#    Residual    97.83333     6    16.30556
#       Total         388    11    35.27273

Julia に翻訳--100 Jonckheere 検定

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

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

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

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

R の table() に相当する関数は 2 つくらいあるのだけど，table 関数ほど使い勝手がよくないので，自分で書いた。

==========#

using Rmath

function jonckheere(x, g; correct=false, alternative = "increasing" #= or "decreasing" =#)
    sgn = alternative == "increasing" ? "<" : ">"
    alternative = "control $sgn= treat_1 $sgn= ... $sgn= treat_n"
    gv, ni = table(g) # ni[1]; names(ni)
    a = length(ni)
    n = sum(ni)
    sumSij = sumTij = 0  # i = 1; j = 2; typeof(gv)
    x = vec(transpose(x))
    for i in 1:a - 1
        di = x[g .== gv[i]]
        for j in i + 1:a
            dj = x[g .== gv[j]]
            if sgn == "<"
                sumTij += sum(di .< dj')
            elseif sgn == ">"
                sumTij += sum(di .> dj')
            end
            sumSij += sum(di .== dj')
        end
    end
    val, tau = table(x)
    V = (n*(n-1)*(2*n+5) - sum(ni .* (ni .- 1) .* (2 .* ni .+ 5)) -
        sum(tau .* (tau .- 1) .* (2 .* tau .+ 5)))/72 +
        sum(ni .* (ni .- 1) .* (ni .- 2)) * sum(tau .* (tau .- 1) .* (tau .- 2)) /
        (36*n*(n-1)*(n-2)) + sum(ni .* (ni .- 1)) * sum(tau .* (tau .- 1)) /
        (8*n*(n-1))
    J = sumTij + sumSij / 2
    E = (n^2 - sum(ni .^ 2)) / 4
    Z = (abs(J - E) - 0.5correct) / sqrt(V)
    pvalue = pnorm(Z, false)
    println("J = $J,  E(J) = $E,  V(J) = $V,  Z = $Z,  p value = $pvalue")
    Dict(:J => J, :E => E, :V => V, :Z => Z, :pvalue => pvalue)
end

function table(x)
    indices = sort(unique(x))
    counts = zeros(Int, length(indices))
    for i = 1:length(indices)
        counts[i] = count(isequal(indices[i]), x)
    end
    return indices, counts
end

x = [153 153 152 156 158 151 151 150 148 157
     158 152 152 152 151 151 157 147 155 146
     153 146 138 152 140 146 156 142 147 153
     137 139 141 141 143 133 147 144 151 156];
g = vcat(repeat([1], 10), repeat([2], 10), repeat([3], 10), repeat([4], 10));

jonckheere(x, g, alternative="decreasing")
# J = 446.5,  E(J) = 300.0,  V(J) = 1705.0775978407557,  Z = 3.547852454888916,  p value = 0.00019419286167098897

jonckheere(x, g, correct=true, alternative="decreasing")
# J = 446.5,  E(J) = 300.0,  V(J) = 1705.0775978407557,  Z = 3.535743743438783,  p value = 0.00020331446444516517

Julia に翻訳--099 ファン・デル・ワーデン検定

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

ファン・デル・ワーデン検定
http://aoki2.si.gunma-u.ac.jp/R/vdw-test.html

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

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

x = vcat(x, y) を append!(x, y) にすると，呼出元の x が変更されてしまう。

==========#

using Rmath, StatsBase

function vdwtest(x, y)
    n1 = length(x)
    n2 = length(y)
    n = n1 + n2
    x = vcat(x, y) # 呼出元に影響を及ぼさない
    x = qnorm.(tiedrank(x) ./ (n + 1))
    S = abs(sum(x[1:n1]))
    V = n1 * n2 / (n ^ 2 - n) * sum(x .^ 2)
    Z = S / sqrt(V)
    pvalue = pnorm(Z, false) * 2
    println("S = $S,  V = $V,  Z = $Z,  p value = $pvalue")
    Dict(:S => S, :V => V, :Z => Z, :pvalue => pvalue)
end

gr1 = [1.2, 1.9, 2.5, 6.7];
gr2 = [1.5, 3.1, 10.5];
vdwtest(gr1, gr2)
# S = 0.7944989941443014,  V = 1.0741549510392339,  Z = 0.7665842364210947,  p value = 0.44332875078500844

Julia に翻訳--098 マン・ホイットニーの U 検定

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

マン・ホイットニーの U 検定
http://aoki2.si.gunma-u.ac.jp/R/u-test.html

ファイル名: utest.jl
関数名:    utest(tbl; correct = true)
          utest(x, y; correct = true)

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

==========#

using StatsBase, Rmath

function utest(tbl; correct = true)
    nr, nc = size(tbl)
    nr == 2 || error("2 x m matrix")
    x = []
    y = []
    for i = 1:nc
        append!(x, repeat([i], tbl[1, i]))
        append!(y, repeat([i], tbl[2, i]))
    end
    utest(x, y; correct)
end

function utest(x, y; correct = true)
    n1 = length(x)
    n2 = length(y)
    n = n1 + n2
    xy = vcat(x, y)
    r = tiedrank(xy)
    U1 = n1 * n2 + n1 * (n1 + 1) / 2 - sum(r[1:n1])
    ties = calctie(r)
    U = min(U1, n1 * n2 - U1)
    V = n1 * n2 *(n ^ 3 - n - ties) / 12 / (n ^ 2 - n)
    E = n1 * n2 / 2
    Z = (abs(U - E) - (correct ? 0.5 : 0)) / sqrt(V)
    P = pnorm(Z, false) * 2
    println("U = $U,  E(U) = $E,  V(U) = $V,  Z = $Z,  p value = $P")
    Dict(:U => U, :EU => E, :VU => V, :Z => Z, :pvalue => P)
end

function calctie(r)
    t = [count(isequal(s), r) for s in Set(r)]
    sum(t .^ 3 .- t)
end

二つのデータベクトルを与えるとき

A = vcat(repeat([1], 9), repeat([2], 12), repeat([3], 6), repeat([4], 3))
B = vcat(repeat([1], 4), repeat([2], 9), repeat([3], 11), repeat([4], 5))
utest(A, B, correct=false)
# U = 310.5,  E(U) = 435.0,  V(U) = 3993.5593220338983,  Z = 1.9701045835900841,  p value = 0.04882638572318423

utest(A, B)
# U = 310.5,  E(U) = 435.0,  V(U) = 3993.5593220338983,  Z = 1.9621925169893206,  p value = 0.04974007481574249

2 x C の集計表を与えるとき

tbl = [10 5 1; 4 7 3]
utest(tbl, correct=false)
# U = 70.0,  E(U) = 112.0,  V(U) = 481.9862068965517,  Z = 1.913074951284128,  p value = 0.05573845794420809

utest(tbl)
# U = 70.0,  E(U) = 112.0,  V(U) = 481.9862068965517,  Z = 1.8903002494831267,  p value = 0.05871781545216507

Julia に翻訳--097 二次データで三群以上の等分散性の検定

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

三群以上の等分散性の検定（二次データ）
http://aoki2.si.gunma-u.ac.jp/R/my-bartlett-test.html

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

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

==========#

using Rmath

function mybartletttest(n, u)
    (length(n) == length(u) && all(n .> 1) && all(u .> 0) && all(floor.(n) .== n)) || error("argument error")
    ng = length(n)
    temp1 = n .- 1
    temp2 = sum(temp1)
    chisq0 = temp2 * log(sum(temp1 .* u) / temp2) - sum(temp1 .* log.(u))
    co = 1 + (sum(1 ./ temp1) - 1/temp2) / (3 * ng - 3)
    chisq = chisq0 / co
    df = ng - 1
    pvalue = pchisq(chisq, df, false)
    println("chisq. = $chisq,  df = $df,  p value = $pvalue")
    Dict(:chisq => chisq, :df => df, :pvalue => pvalue)
end

n = [11, 25, 22]
u = [23.7, 25.7, 24.1]
mybartletttest(n, u)
# chisq. = 0.03290799098798558,  df = 2,  p value = 0.9836806320912985

Julia に翻訳--096 一元配置分散分析，三群以上の平均値の差の検定，二次データ

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

一元配置分散分析（三群以上の平均値の差の検定，二次データ）
http://aoki2.si.gunma-u.ac.jp/R/my-oneway-ANOVA.html

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

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

==========#

using Rmath

function myonewayanova(n, m, u; varequal=false)
    ng = length(n)
    dfb = ng - 1
    (ng > 1 && length(m) == ng && length(u) == ng && all(n .> 1) && all(floor.(n) == n) && all(u .> 0)) || error("parameter error")
    if varequal
        nc = sum(n)
        sw = sum(u .* (n .- 1))
        sb = sum(n .* (m .- sum(n .* m)/nc) .^ 2)
        dfw = nc - ng
        msb, msw = sb / dfb, sw / dfw
        f = msb / msw
        pvalue = pf(f, dfb, dfw, false)
    else
        w = n ./ u
        sumw = sum(w)
        m0 = sum(w .* m) / sumw
        temp = sum((1 .- w ./ sumw) .^ 2 ./ (n .- 1)) / (ng^2-1)
        f = sum(w .* (m .- m0) .^ 2) / ((ng - 1) * (1 + 2 * (ng - 2) * temp))
        dfw = 1 / (3 * temp)
        pvalue = pf(f, dfb, dfw, false)
    end
    println("F = $f,  dfb = $dfb,  dfw = $dfw,  p value = $pvalue")
    return Dict(:F => f, :dfb => dfb, :dfw => dfw, :pvalue => pvalue)
end

n = [8, 11, 22, 6];
m = [135.83, 160.49, 178.35, 188.06];
SD = [19.59, 12.28, 15.01, 9.81];
u = SD .^ 2;

myonewayanova(n, m, u, varequal=true)
# F = 20.82824265645924,  dfb = 3,  dfw = 43,  p value = 1.7374840062619836e-8

myonewayanova(n, m, u)
# F = 17.56460925505331,  dfb = 3,  dfw = 16.50402543036444,  p value = 2.191921147555868e-5

Julia に翻訳--095 対応のある代表値の差の検定，ウィルコクソンの符号付順位和検定

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

対応のある代表値の差の検定（ウィルコクソンの符号付順位和検定）
http://aoki2.si.gunma-u.ac.jp/R/wilcox-paired.html

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

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

R の wilcox.test() に含まれる Wilcoxon Signed Rank Test では，同値がなくても n < 50 のときしか exact な検定を行えないようになっている。
しかし，使用される psignrank() は n < 1039 まで可能である。

このプログラムでは，可能な限り exact な検定を行い，同値がある場合には漸近検定を行う。exact な検定が行われた場合も，漸近検定の結果は表示される。

==========#

using StatsBase, Rmath

function calctie(r)
    t = [count(isequal(s), r) for s in Set(r)]
    sum(t .^ 3 .- t)
end

function wilcoxonsignedranktest(x, y; correct=true)
    x .-= y
    n0 = length(x)
    x = x[x .!= 0]
    n = length(x)
    r = tiedrank(abs.(x))
    w = sum(r[x .> 0])
    ties = calctie(r)
    z = w - n * (n + 1)/4
    if n < 1039 && ties == 0 && n0 == n
        pvalue = z > 0 ? psignrank(w - 1, n, false) : psignrank(w, n)
        pvalue = min(2 * pvalue, 1)
        println("w = $w,  p value = $pvalue (exact)")
    end
    method = "asymptotic test"
    sigma = sqrt(n * (n + 1) * (2 * n + 1)/24 - ties/48)
    correction = 0
    if correct
        correction = 0.5sign(z)
        method = method * " with continuity correction"
    end
    z = (z - correction) / sigma
    pvalue = 2 * min(pnorm(z), pnorm(z, false))
    println("z = $z,  p value = $pvalue ($method)")
end

同値がある場合は漸近検定

x = [269, 230, 365, 282, 295, 212, 346, 207, 308, 257];
y = [273, 213, 383, 282, 297, 213, 351, 208, 294, 238];
wilcoxonsignedranktest(x, y)
# z = 0.0,  p value = 1.0 (asymptotic test with continuity correction)

x = [48.1, 51.2, 61.9, 40.7, 49.2, 53.9, 61.0, 33.6, 21.5, 72.3]
y = [58.1, 63.2, 61.8, 39.0, 56.2, 63.9, 71.2, 62.9, 48.1, 66.6]
wilcoxonsignedranktest(x, y)
# z = -2.141909506235002,  p value = 0.03220076475250472 (asymptotic test with continuity correction)

同値がない場合は常に exact な検定と漸近検定の両方の結果が表示される

x = [51.4, 42.2, 61.7, 47.9, 50.0, 63.6, 53.0, 34.2, 56.2, 49.9];
y = [43.8, 54.8, 51.6, 54.0, 61.0, 57.4, 67.0, 77.6, 65.9, 62.6];
wilcoxonsignedranktest(x, y)
# w = 10.0,  p value = 0.083984375 (exact)
# z = -1.732800450887927,  p value = 0.0831311427008039 (asymptotic test with continuity correction)

同値さえなければ，サンプルサイズが 1038 以下ならば exact な検定が行われる。

using Random; Random.seed!(123);
n = 1038
x = randn(n)
y = randn(n)
wilcoxonsignedranktest(x, y)
# w = 280908.0,  p value = 0.24279165850242684 (exact)
# z = 1.1683136400265737,  p value = 0.24268027537889492 (asymptotic test with continuity correction)

 

Julia に翻訳--094 対応のある平均値の差の検定，原データ，二次データ

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

対応のある平均値の差の検定（原データ，二次データ）
http://aoki2.si.gunma-u.ac.jp/R/paired-t-test.html

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

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

==========#

using Statistics, Rmath

function pairedttest(x, y, ux2=NaN, uy2=NaN, r=NaN, n=NaN)
    if isnan(ux2)
        n = length(x)
        x .-= y
        t = abs(mean(x)) / (std(x) / sqrt(n))
    else
        t = abs(x .- y) / sqrt((ux2 + uy2 - 2 * r * sqrt(ux2 * uy2)) / n)
    end
    df = n - 1
    pvalue = pt(t, df, false)*2
    println("t = $t,  df = $df,  p value = $pvalue")
    Dict(:t => t, :df => df, :pvalue => pvalue)
end

x = [-1.1, -1.1, -0.5, -3.3, 0.0, -1.1, 1.6, 1.6, -0.7]
y = [-1.5, -0.7, -0.8, -2.7, 0.1, -0.4, 0.5, 0.2, -0.3]

二次データから検定する場合

pairedttest(mean(x), mean(y), var(x), var(y), cor(x, y), length(x))
# t = 0.44499415948998505,  df = 8,  p value = 0.6681179327157327

一次データから検定する場合

pairedttest(x, y)
# t = 0.44499415948998505,  df = 8,  p value = 0.6681179327157327

 

Julia に翻訳--093 二群の等分散性の検定，二次データ

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

二群の等分散性の検定（二次データ）
http://aoki2.si.gunma-u.ac.jp/R/my-var-test.html

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

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

単純なので 3 行をまとめて 1 行に書いてみるが，あまりお勧めのプログラミングスタイルとは思えない。

==========#

using Rmath

function myvartest(nx, vx, ny, vy)
    if vx > vy
        F, df1, df2 = vx / vy, nx - 1, ny - 1
    else
        F, df1, df2 = vy / vx, ny - 1, nx - 1
    end
    pvalue = 2 * pf(F, df1, df2, false)
    println("F = $F,  df1 = $df1,  df2 = $df2,  p value = $pvalue")
    Dict(:F => F, :df1 => df1, :df2 => df2, :pvalue => pvalue)
end

myvartest(42, 1.1, 63, 3.2)
# F = 2.909090909090909,  df1 = 62,  df2 = 41,  p value = 0.00044318863143253935