Julia に翻訳--134 マン・ホイットニーの U 検定，exact test

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

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

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

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

global は不要だった。

==========#

using Rmath, Printf

function exactmw(x, y = NaN; exact = true)
    function found()
        hh = abs(temp2 - sum(um[1, :] .* r))
        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 asymptotic()
        z = stat_val / sqrt(n1n2 / (n * (n - 1)) * (n ^ 3 - n - sum(ct .^ 3 .- ct)) / 12)
        @printf("U = % g, Z = % g, P 値 = % g\n", n1n2 / 2 - stat_val, z, pnorm(z, false) * 2)
    end

    if ndims(x) == 2
        t = x
    else
        nx = length(x)
        ny = length(y)
        indexx, indexy, t = table(rep(1:2, [nx, ny]), vcat(x, y))
    end
    EPSILON = 1e-10
    EXPMAX = 1e100
    log_expmax = log(EXPMAX)
    nr, nc = size(t)
    nr == 2 || error("nrow(table) wasn't 2")
    um = zeros(Int, nr, nc)
    rt = sum(t, dims=2)
    ct = transpose(sum(t, dims=1))
    n1 = rt[1]
    n2 = rt[2]
    n1n2 = n1 * n2
    n = n1 + n2
    r = cumsum(vcat(0, ct[1:nc-1])) .+ (ct .+ 1) ./ 2
    temp = n1n2 + n1 * (n1 + 1) / 2
    temp2 = temp - n1n2 / 2
    stat_val = abs(temp2 - sum(t[1, :] .* r))
    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

x = [5 3 2 1
     2 3 1 2];
exactmw(x)
# U =  33.5, Z =  0.907299, P 値 =  0.364248
# 正確な P 値 =  0.3837421608
# 査察した分割表の個数は 91 個

2 つのデータベクトルが与えられる場合

x = [78, 64, 75, 45, 82];
y = [110, 70, 53, 53];
exactmw(x, y)
[1 0 1 0 1 1 1 0; 0 2 0 1 0 0 0 1]
# U =  9, Z =  0.245976, P 値 =  0.805701
# 正確な P 値 =  0.8492063492
# 査察した分割表の個数は 91 個

Julia に翻訳--133 フィッシャーの正確確率検定

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

フィッシャーの正確確率検定
http://aoki2.si.gunma-u.ac.jp/R/fisher.html

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

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

hybrid は省略した。

==========#

using Rmath, Printf

function fisher(x, y = NaN; exact = true, method = "Fisher") # or method = "Pearson"
    function found()
        if method == "Fisher"
            nprod = temp - sum(perm_table[um .+ 1])
            if nprod <= criterion + EPSILON
                nntrue += exp(nprod - nntrue2 * log_expmax)
                while nntrue >= EXPMAX
                    nntrue /= EXPMAX
                    nntrue2 += 1
                end
            end
        else
            hh = sum((um .- ex) .^ 2 ./ ex)
            if hh >= stat_val || abs(hh - stat_val) <= EPSILON
                nprod = temp - sum(perm_table[um .+ 1])
                nntrue += exp(nprod - nntrue2 * log_expmax)
                while nntrue >= EXPMAX
                    nntrue /= EXPMAX
                    nntrue2 += 1
                end
            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("% s の方法による，正確な P 値 = % .10g\n", method, nntrue / denom * EXPMAX ^ (nntrue2 - denom2))
        @printf("査察した分割表の個数は % s 個\n", ntables)
    end

    function chisqtest(t)
        return sum((t .- ex) .^ 2 ./ ex)
    end
    function asymptotic()
        @printf("カイ二乗値 = % g, 自由度 = % i, P 値 = % g\n", stat_val, df, pchisq(stat_val, df, false))
    end

    if ndims(x) == 2
        t = x
    else
        indexy, indexx, 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 = sum(t, dims=1)
    n = sum(t)
    ex = (rt * ct) ./ n
    stat_val = chisqtest(t)
    df = (nr - 1) * (nc - 1)
    asymptotic()
    if exact
        perm_table = cumsum(vcat(0, log.(1:n + 1)))
        temp = sum(perm_table[rt .+ 1])
        criterion = temp - sum(perm_table[t .+ 1])
        ntables = nntrue = nntrue2 = 0
        exacttest()
    end
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 = [5 3 2 1
     4 3 5 2
     2 3 1 2]

fisher(x)
# カイ二乗値 =  3.39631, 自由度 =  6, P 値 =  0.757711
# Fisher の方法による，正確な P 値 =  0.8091124268
# 査察した分割表の個数は 24871 個

fisher(x, method = "Pearson")
# カイ二乗値 =  3.39631, 自由度 =  6, P 値 =  0.757711
# Pearson の方法による，正確な P 値 =  0.7878188077
# 査察した分割表の個数は 24871 個

2 つのベクトルを与える場合

x = [1, 1, 1, 3, 2, 1, 1, 3, 3, 2, 3, 3, 1, 1, 3, 1, 1, 1, 2, 1,
    1, 2, 2, 3, 1, 1, 2, 3, 3, 2];
y = [1, 2, 2, 3, 2, 1, 2, 2, 2, 3, 3, 2, 2, 1, 2, 1, 2, 2, 3, 1,
    1, 1, 1, 3, 1, 2, 2, 3, 1, 3];

fisher(x, y)
# カイ二乗値 =  9.17504, 自由度 =  4, P 値 =  0.0568702
# Fisher の方法による，正確な P 値 =  0.03320873103
# 査察した分割表の個数は 1353 個

fisher(x, y, method="Pearson")
# カイ二乗値 =  9.17504, 自由度 =  4, P 値 =  0.0568702
# Pearson の方法による，正確な P 値 =  0.05323058543
# 査察した分割表の個数は 1353 個

Julia に翻訳--132 適合度の検定，exact test

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

適合度の検定（exact test）
http://aoki2.si.gunma-u.ac.jp/R/gft.html

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

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

==========#

using Distributions, Printf, SpecialFunctions

function gft(o, p = repeat([1 / length(o)], length(o)))
    function gen_tab(y)
        if y == 1
            x2 = sum((o .- e) .^ 2 ./ e)
            if x2 >= x0 || abs(x2 - x0) <= EPSILON
                w2 += exp(temp - sum(fact[o.+1])) * prod(p .^ o)
            end
        else
            gen_tab(y - 1)
            while o[1] > 0
                o[y] += 1
                o[1] -= 1
                gen_tab(y - 1)
            end
            o[1] += o[y]
            o[y] = 0
        end
    end

    EPSILON = 1e-10
    total = sum(o)
    p ./= sum(p)
    e = p .* total
    x0 = sum((o .- e) .^ 2 ./ e)
    n = length(o)
    p_val = ccdf(Chisq(n - 1), x0)
    @printf("カイ二乗値は %g，自由度は %i，P値は %g\n", x0, n - 1, p_val)
    fact = logfactorial.(0:total)
    temp = fact[total+1]
    w2 = 0.0
    o = zeros(Int64, n)
    o[1] = total
    gen_tab(n)
    @printf("正確なP値は %g\n", w2)
end

gft([1, 2, 3, 4, 5, 6])
# カイ二乗値は 5，自由度は 5，P値は 0.41588
# 正確なP値は 0.467389

gft([10,1,1,1,1,1,2])
# カイ二乗値は 27.8824，自由度は 6，P値は 9.88791e-05
# 正確なP値は 0.000223375

 

Julia に翻訳--131 シャーリー・ウィリアムズ法，代表値の多重比較

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

シャーリー・ウィリアムズの方法による代表値の多重比較
http://aoki2.si.gunma-u.ac.jp/R/Shirley-Williams.html

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

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

# t = (M - sum(r[group .== 1]) / ni[1]) / sqrt(V * (1 / ni[p] + 1 / ni[1]))
の行は多分バグ。（結果に違いはないが，Julia だとエラーになる）

tapply() と rep() を書いてみた。

==========#

using StatsBase, Rmath

function shirleywilliams(data, group; method = "up") # or method = "down"
    func = method == "down" ? minimum : maximum
    ni = tapply(data, group, length)
    a = length(ni)
    s = zeros(a - 1)
    for p in a:-1:2 # p = 3
        select = 1 .<= group .<= p
        r = tiedrank(data[select])
        g = group[select]
        M = func(cumsum(reverse(tapply(r, g, sum))[1:p-1]) ./ cumsum(reverse(ni[2:p])))
        N = sum(ni[1:p])
        V = (sum(r .^ 2) - N * (N + 1) ^ 2 / 4) / (N - 1)
        # t = (M - sum(r[group .== 1]) / ni[1]) / sqrt(V * (1 / ni[p] + 1 / ni[1]))
        t = (M - sum(r[g .== 1]) / ni[1]) / sqrt(V * (1 / ni[p] + 1 / ni[1]))
        s[p - 1] = method == "down" ? -t : t
    end
    t = reverse(s)
    [println("$(a - i + 1):  t = $t") for (i, t) in enumerate(t)]
    return t
end

function table(x) # indices が少ないとき
    indices = sort(unique(x))
    counts = zeros(Int, length(indices))
    for i in indexin(x, indices)
        counts[i] += 1
    end
    return indices, counts
end

function tapply(x, g, func)
    indices = sort(unique(g))
    [func(x[g .== i]) for i in indices]
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 rep(x; each::Int64)
    vcat([repeat([x[i]], each) for i = 1:length(x)]...)
end

function rep(x, n::Int64)
    repeat(x, n)
end

data = [
  13, 23,  8, 17, 25, 34, 18, 26, 10, 28, 18, 21,
  26, 22, 30, 38, 15, 24, 18, 11, 21, 30, 31, 23,
  22, 10, 29, 37, 22, 13, 29, 28, 21, 16, 21, 26,
  26, 34, 30, 45, 17, 19, 27, 18, 36, 24, 25, 31
];
group = rep(1:4, each=12);

shirleywilliams(data, group, method = "up")
# 4:  t = 2.1898989426889655
# 3:  t = 1.0917905788994513
# 2:  t = 1.2434675147149

Julia に翻訳--130 スティール・ドゥワス法，代表値の多重比較

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

スティール・ドゥワスの方法による代表値の多重比較
http://aoki2.si.gunma-u.ac.jp/R/Steel-Dwass.html

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

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

Rmath の ptukey() にバグがある。
lower_tail = false を指定すると NaN を返す。
回避策としては，
ptukey(q, nmeans, df, false)
を
1 - ptukey(q, nmeans, df)
とする。

==========#

using Combinatorics, StatsBase, Rmath, Printf

function steeldwass(data, group)
    index, ni = table(group)
    ng = length(ni)
    k = binomial(ng, 2)
    t = zeros(k)
    pair = fill("", k)
    for (n, (i, j)) in enumerate(combinations(1:ng, 2))
        r = tiedrank(vcat(data[group .== i], data[group .== j]))
        R = sum(r[1:ni[i]])
        N = ni[i] + ni[j]
        E = ni[i] * (N + 1) / 2
        V = ni[i] * ni[j] / (N * (N - 1)) * (sum(r .^ 2) - N * (N + 1) ^ 2 / 4)
        pair[n] = string(i) * ":" * string(j)
        t[n] = abs(R - E) / sqrt(V)
    end
    p = 1 .- ptukey.(t * sqrt(2), ng, Inf)
    @printf(" pair  %10s  %10s\n", "t value", "p value")
    for i in 1:k
        @printf("%5s  %10.6f  %10.6f\n", pair[i], t[i], p[i])
    end
    Dict(:pair => pair, :t => t, :pvalue => p)
end

function table(x) # indices が少ないとき
    indices = sort(unique(x))
    counts = zeros(Int, length(indices))
    for i in indexin(x, indices)
        counts[i] += 1
    end
    return indices, counts
end

data = [
    6.9, 7.5, 8.5, 8.4, 8.1, 8.7, 8.9, 8.2, 7.8, 7.3, 6.8,
    9.6, 9.4, 9.5, 8.5, 9.4, 9.9, 8.7, 8.1, 7.8, 8.8,
    5.7, 6.4, 6.8, 7.8, 7.6, 7.0, 7.7, 7.5, 6.8, 5.9,
    7.6, 8.7, 8.5, 8.5, 9.0, 9.2, 9.3, 8.0, 7.2, 7.9, 7.8
    ];
group = vcat([repeat([i], n) for (i, n) in enumerate([11, 10, 10, 11])]...);
steeldwass(data, group)
# pair     t value     p value
#  1:2    2.680234    0.036960
#  1:3    2.539997    0.053981
#  1:4    1.282642    0.574012
#  2:3    3.746076    0.001031
#  2:4    2.046776    0.170966
#  3:4    3.384456    0.003977

Julia に翻訳--129 スティール法，代表値の多重比較

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

スティールの方法による代表値の多重比較
http://aoki2.si.gunma-u.ac.jp/R/Steel.html

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

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

==========#

using StatsBase, Printf

function steel(data, group)
    index, ni = table(group)
    a = length(ni)
    control = data[group .== 1]
    n1 = length(control)
    t = zeros(a)
    rho = sum(n1 .== ni) == a ? 0.5 : getrho(ni)
    for i in 2:a
        r = tiedrank(vcat(control, data[group .== i]))
        R = sum(r[1:n1])
        N = n1 + ni[i]
        E = n1 * (N + 1) / 2
        V = n1 * ni[i] / N / (N - 1) * (sum(r .^ 2) - N * (N + 1) ^ 2 / 4)
        t[i] = abs(R - E) / sqrt(V)
    end
    @printf(" pair  %10s  %10s\n", "t", "rho")
    for i = 2:a
        @printf("%2d:%2d  %10.6f  %10.6f\n", 1, i, t[i], rho)
    end
    Dict(:pair => string.(repeat([1], a-1)) .* ":" .* string.(2:a),
         :t => t[2:a], :rho => rho)
end

function table(x) # indices が少ないとき
    indices = sort(unique(x))
    counts = zeros(Int, length(indices))
    for i in indexin(x, indices)
        counts[i] += 1
    end
    return indices, counts
end

func(x, y, ni1) = sqrt(x / (x + ni1) * y / (y + ni1))

function getrho(ni)
    k = length(ni)
    rho = func.(ni, transpose(ni), ni[1])
    [rho[i, i] = 0 for i =1:k]
    sum(rho[2:k, 2:k]) / (k - 2) / (k - 1)
end

data = [
  50, 55, 65, 63, 60, 68, 69, 60, 52, 49,
  80, 86, 74, 66, 79, 81, 70, 62, 60, 72,
  42, 48, 58, 63, 62, 55, 63, 60, 53, 45
];

group = vcat([repeat([i], 10) for i in 1:3]...);
steel(data, group)
# pair           t         rho
# 1: 2    2.952566    0.500000
# 1: 3    1.175674    0.500000