Julia に翻訳--138 Brown-Forsythe 検定，三群以上の場合の等分散性の検定

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

Brown-Forsythe 検定（三群以上の場合の等分散性の検定）
http://aoki2.si.gunma-u.ac.jp/R/Brown-Forsythe-test.html

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

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

R の split() と sapply() の間単版を Julia で書いてみた。

==========#

using Statistics, Rmath

function brownforsythetest(x, group)
    d = split2(x, group)
    df1 = length(d) - 1
    ni = sapply(d, length)
    mi = sapply(d, mean)
    ui = sapply(d, var)
    ci = (1 .- ni ./ sum(ni)) .* ui
    FBF = sum(ni .* (mi .- mean(x)) .^ 2) / sum(ci)
    C = ci ./ sum(ci)
    df2 = 1 / sum(C .^ 2 / (ni .- 1))
    p = pf(FBF, df1, df2, false)
    println("F = $FBF,  df1 = $df1,  df2 = $df2,  p value = $p")
    Dict(:F => FBF, :df1 => df1, :df2 => df2, :pvalue => p)
end

function split2(x, g)
    index = sort(unique(g))
    y = Array{Array{Float64,1},1}(undef, length(index))
    [y[i] = x[g .== i] for i in index]
end

sapply(x, FUNC) = map(FUNC, x)

function rep(x; each::Int64)
    vcat([repeat([x[i]], each) for i = 1:length(x)]...)
end

x = [3, 3, 4, 2, 5, 2, 3, 4, 8, 8, 5, 6];
g = rep(1:3, each = 4);
brownforsythetest(x, g)
# F = 10.854545454545452,  df1 = 2,  df2 = 7.606873428331936,  p value = 0.00590981731767559

Julia に翻訳--137 Levene 検定，三群以上の場合の等分散性の検定

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

Levene 検定（三群以上の場合の等分散性の検定）
http://aoki2.si.gunma-u.ac.jp/R/levene-test.html

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

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

==========#

using Statistics, Rmath

function levenetest(x, group; method = "mean") # or method = "median"
    n = length(x)
    ni = tapply(x, group, length)
    k = length(ni)
    FUNC = method == "mean" ? mean : median
    x = abs.(x .- tapply(x, group, FUNC)[group])
    sw = sum((ni .- 1) .* tapply(x, group, var))
    dfw = n - k
    dfb = k - 1
    f = ((var(x) * (n - 1) - sw) / dfb) / (sw / dfw)
    P = pf(f, dfb, dfw, false)
    println("等分散性の検定（$method で調整）")
    println("F = $f,  dfb = $dfb,  dfw = $dfw,  p value = $P")
    Dict(:F => f, :dfb => dfb, :dfw => dfw, :pvalue => P, :method => method)
end

function tapply(x, g, func)
    indices = sort(unique(g))
    [func(x[g .== i]) for i in indices]
end

function rep(x; each::Int64)
    vcat([repeat([x[i]], each) for i = 1:length(x)]...)
end

x = [3, 3, 4, 2, 5, 2, 3, 4, 8, 8, 5, 6]
g = rep(1:3, each = 4)

levenetest(x, g)
# 等分散性の検定（mean で調整）
# F = 2.0999999999999996,  dfb = 2,  dfw = 9,  p value = 0.17844673330519398

levenetest(x, g, method = "median")
# 等分散性の検定（median で調整）
# F = 1.9090909090909094,  dfb = 2,  dfw = 9,  p value = 0.203644351599495

Julia に翻訳--136 一元配置分散分析，exact test

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

一元配置分散分析(exact test)
http://aoki2.si.gunma-u.ac.jp/R/exact-oneway-test.html

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

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

==========#

using Statistics, Rmath, Printf

function exactonewaytest(x, y=[]; permutation = true)
    function found()
        hh = performtest(um)
        if hh <= pvalue + EPSILON
            nprod = sum(perm_table[rt .+ 1]) - sum(perm_table[um .+ 1])
            nntrue += exp(nprod - nntrue2 * log_expmax)
            while nntrue >= EXPMAX
                nntrue /= EXPMAX
                nntrue2 += 1
            end
        end
        ntables += 1
    end

    function search(x, y)
        if y == 1
            found()
        elseif x == 1
            t = um[1, 1] - um[y, 1]
            if t >= 0
                um[1, 1] = t
                search(nc, y - 1)
                um[1, 1] += um[y, 1]
            end
        else
            search(x - 1, y)
            while um[y, 1] != 0 && um[1, x] != 0
                um[y, x] += 1
                um[y, 1] -= 1
                um[1, x] -= 1
                search(x - 1, y)
            end
            um[y, 1] += um[y, x]
            um[1, x] += um[y, x]
            um[y, x] = 0
        end
    end

    function permutationtest()
        denom2 = 0
        denom = perm_table[n + 1] - sum(perm_table[ct .+ 1])
        while denom > log_expmax
            denom -= log_expmax
            denom2 += 1
        end
        denom = exp(denom)
        um[:, 1] = rt
        um[1, :] = ct
        search(nc, nr)
        pvalue = nntrue / denom * EXPMAX ^ (nntrue2 - denom2)
        @printf("並べ替え検定による P 値 = % .10g\n", pvalue)
        @printf("査察した分割表の個数は % s 個\n", ntables)
        return pvalue
    end

    function performtest(u)
        x = Array{Array{Float64,1},1}(undef, nr)
        for i in 1:nr
            x[i] = score[u[i, :] .> 0]
        end
        onewaytest(x)
    end

    function simpleonewaytest()
        @printf("観察値による一元配置分散分析の P 値 = % g\n", pvalue)
    end

    if length(y) == 0
        ni = map(length, x)
        y = rep(1:length(ni), ni)
        indexy, score, t = table(y, vcat(x...))
    else
        indexy, score, t = table(y, x)
    end
    EPSILON = 1e-10
    EXPMAX = 1e100
    log_expmax = log(EXPMAX)
    nr, nc = size(t)
    um = zeros(Int, nr, nc)
    rt = sum(t, dims=2)
    ct = transpose(sum(t, dims=1))
    n = sum(t)
    q = cumsum(vcat(0, ct[1:nc-1])) .+ (ct .+ 1) .* 0.5
    half = (n + 1) * 0.5
    pvalue = performtest(t)
    simpleonewaytest()
    if permutation
        perm_table = cumsum(vcat(0, log.(1:n + 1)))
        ntables = nntrue = nntrue2 = 0
        permutationtest()
    end
end

function onewaytest(x)
    n = map(length, x)
    m = map(mean, x)
    v = map(var, x)
    w = n ./ v
    sumw = sum(w)
    tmp = sum((1 .- (w / sumw)) .^ 2 ./ (n .- 1)) / (nr ^ 2 - 1)
    m0 = sum(w .* m) / sumw
    F = sum(w .* (m .- m0) .^ 2) / ((nr - 1) * (1 + 2 * (nr - 2) * tmp))
    df1 = nr - 1
    df2 = 1 / (3 * tmp)
    pf(F, df1, df2, false)
end

function rep(x, n::Array{Int64,1})
    length(x) == length(n) || error("length(x) wasn't length(n)")
    vcat([repeat([x[i]], n[i]) for i = 1:length(x)]...)
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

x = [[36.7, 52.4, 65.8], [45.7, 61.9, 65.3], [52.6, 76.6, 81.3]];
exactonewaytest(x)
# 観察値による一元配置分散分析の P 値 =  0.440037
# 並べ替え検定による P 値 =  0.4107142857
# 査察した分割表の個数は 1680 個

 

Julia に翻訳--135 クラスカル・ウォリス検定，exact test

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

クラスカル・ウォリス検定(exact test)
http://aoki2.si.gunma-u.ac.jp/R/exact-kw.html

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

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

global は不要

==========#

using Rmath, Printf

function exactkw(x, y = []; exact = true)
    function found()
        hh = sum((um * q) .^ 2 ./ rt)
        if hh >= stat_val || abs(hh - stat_val) <= EPSILON
            nprod = sum(perm_table[rt .+ 1]) - sum(perm_table[um .+ 1])
            nntrue += exp(nprod - nntrue2 * log_expmax)
            while nntrue >= EXPMAX
                nntrue /= EXPMAX
                nntrue2 += 1
            end
        end
        ntables += 1
    end

    function search(x, y)
        if y == 1
            found()
        elseif x == 1
            t = um[1, 1] - um[y, 1]
            if t >= 0
                um[1, 1] = t
                search(nc, y - 1)
                um[1, 1] += um[y, 1]
            end
        else
            search(x - 1, y)
            while um[y, 1] != 0 && um[1, x] != 0
                um[y, x] += 1
                um[y, 1] -= 1
                um[1, x] -= 1
                search(x - 1, y)
            end
            um[y, 1] += um[y, x]
            um[1, x] += um[y, x]
            um[y, x] = 0
        end
    end

    function exacttest()
        denom2 = 0
        denom = perm_table[n + 1] - sum(perm_table[ct .+ 1])
        while denom > log_expmax
            denom -= log_expmax
            denom2 += 1
        end
        denom = exp(denom)
        um[:, 1] = rt
        um[1, :] = ct
        search(nc, nr)
        @printf("正確な P 値 = % .10g\n", nntrue / denom * EXPMAX ^ (nntrue2 - denom2))
        @printf("査察した分割表の個数は % s 個\n", ntables)
    end

    function kwtest(u)
        return sum((u * q) .^ 2 ./ rt)
    end

    function asymptotic()
        chisq = (stat_val * 12 / (n * (n + 1)) - 3 * (n + 1)) / (1 - sum(ct .^ 3 .- ct) / (n ^ 3 - n))
        @printf("カイ二乗値 = % g, 自由度 = % i, P 値 = % g\n", chisq, nr - 1, pchisq(chisq, nr - 1, false))
    end

    if ndims(x) == 2
        t = x
    elseif length(y) == 0
        ni = map(length, x)
        y = rep(1:length(ni), ni)
        indexy, indexx, t = table(y, vcat(x...))
    else
        indexy, indexx, t = table(y, x)
    end
    EPSILON = 1e-10
    EXPMAX = 1e100
    log_expmax = log(EXPMAX)
    nr, nc = size(t)
    rt = sum(t, dims=2)
    ct = transpose(sum(t, dims=1))
    um = zeros(Int, nr, nc)
    n = sum(t)
    q = cumsum(vcat(0, ct[1:nc-1])) .+ (ct .+ 1) .* 0.5
    half = (n + 1) * 0.5
    stat_val = kwtest(t)
    asymptotic()
    if exact
        perm_table = cumsum(vcat(0, log.(1:(n + 1))))
        ntables = nntrue = nntrue2 = 0
        exacttest()
    end
end

function rep(x, n::Array{Int64,1})
    length(x) == length(n) || error("length(x) wasn't length(n)")
    vcat([repeat([x[i]], n[i]) for i = 1:length(x)]...)
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

分割表を与える場合 Array{Int64,2}

x = [5 3 2 1
     4 3 5 2
     2 3 1 2]
exactkw(x)
# カイ二乗値 =  1.32485, 自由度 =  2, P 値 =  0.5156
# 正確な P 値 =  0.5268191237
# 査察した分割表の個数は 24871 個

データベクトルと factor ベクトルを与える場合

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 = rep(1:3, [4, 5, 3]);
exactkw(data, group)
# カイ二乗値 =  5.54872, 自由度 =  2, P 値 =  0.0623895
# 正確な P 値 =  0.0538961039
# 査察した分割表の個数は 27720 個

二重配列で与える場合 Array{Array{Float64,1},1}

x = [[3.42, 3.84, 3.96, 3.76], [3.17, 3.63, 3.47, 3.44, 3.39], [3.64, 3.72, 3.91]]
exactkw(x)
# カイ二乗値 =  5.54872, 自由度 =  2, P 値 =  0.0623895
# 正確な P 値 =  0.0538961039
# 査察した分割表の個数は 27720 個