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