Julia に翻訳--173 Bhapkar 検定，一般化マクネマー検定

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

Python による統計処理　
http://aoki2.si.gunma-u.ac.jp/Python/PythonDir.tar.gz

Bhapkar 検定と一般化マクネマー検定
http://aoki2.si.gunma-u.ac.jp/Python/Bhapkar_test.pdf

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

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

==========#

using Rmath

function bhapkartest(tbl)
    r, c = size(tbl)
    r == c || error("'tbl' != a square matix.'")
    r1 = r - 1
    N = sum(tbl)
    pi = tbl ./ N
    picolSums = vec(sum(pi, dims=1))  # π(+, t)
    pirowSums = sum(pi, dims=2)  # π(s, +)
    d = (picolSums .- pirowSums)[1:r1]
    V0 = zeros(r1, r1)
    V1 = zeros(r1, r1)
    for s = 1:r1
        for t = 1:r1
            if s == t
                V0[s, s] = pirowSums[s] + picolSums[s] - 2 * pi[s, s]
                V1[s, s] = V0[s, s] - (pirowSums[s] - picolSums[s])^2
            else
                V0[s, t] = -pi[s, t] - pi[t, s]
                V1[s, t] = V0[s, t] - (picolSums[s] - pirowSums[s]) * (picolSums[t] - pirowSums[t])
            end
        end
    end
    Z0 = N * d' * inv(V0) * d
    p0 = pchisq(Z0, r1, false)
    Z1 = N * d' * inv(V1) * d
    p1 = pchisq(Z1, r1, false)
    println("Generalized McNemar test (Stuart-Maxwell test)")
    println("chisq. = $Z0,  df = $r1,  p value = $p0")
    println("Bhapkar's test")
    println("chisq. = $Z1,  df = $r1,  p value = $p1")
    Dict(:Z0 => Z0, :Z1 => Z1, :df => r1, :p0 => p0, :p1 => p1)
end

tbl = [1520  266  124   66
        234 1512  432   78
        117  362 1772  205
         36   82  179  492]
a = bhapkartest(tbl)
# Generalized McNemar test (Stuart-Maxwell test)
# chisq. = 11.956569622982544,  df = 3,  p value = 0.007533425054800585
# Bhapkar's test
# chisq. = 11.97572015552566,  df = 3,  p value = 0.007466797469724006

Julia に翻訳--172 Woolf 検定

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

Python による統計処理　
http://aoki2.si.gunma-u.ac.jp/Python/PythonDir.tar.gz

Woolf 検定
http://aoki2.si.gunma-u.ac.jp/Python/Woolf_test.pdf

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

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

Python の　map() と　Julia の map() は違う（あたりまえ）

==========#

using Rmath

function woolftest(x)
    all(x .!= 0) || (x = x + 1 / 2)
    k = size(x, 3)
    oddsratio = [x[1, 1, i] * x[2, 2, i] / (x[1, 2, i] * x[2, 1, i]) for i = 1:k]
    weight = [1 / sum(1 ./ x[:, :, i]) for i = 1:k]
    observed = log.(oddsratio)
    expected = sum(observed .* weight) / sum(weight)
    chisq = sum(weight .* (observed .- expected) .^ 2)
    df = k - 1
    pvalue = pchisq(chisq, df, false)
    println("chisq = $chisq,  df = $df,  p value = $pvalue")
    Dict(:chisq => chisq, :df => df, :pvalue => pvalue)
end

x0 = [9, 95, 7, 1841, 23, 105, 9, 1654, 54, 177, 19, 1863,
      121, 257, 48, 2357, 169, 273, 54, 1778, 269, 324, 88, 1712, 404,
      245, 117, 1324, 406, 225, 152, 967, 372, 132, 106, 526]

x = reshape(x0, 2, 2, 9)

woolftest(x)
# chisq = 26.203356695998252,  df = 8,  p value = 0.0009693542211758595

Julia に翻訳--171 Mood 検定，二群の等分散性の検定

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

Python による統計処理　
http://aoki2.si.gunma-u.ac.jp/Python/PythonDir.tar.gz

Mood 検定（二群の等分散性の検定）
http://aoki2.si.gunma-u.ac.jp/Python/Mood_test.pdf

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

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

R は，特に指定のない場合には，変数は実数型である。Julia はそうではないので，演算中のオーバーフローに注意しないといけない。
わかってはいるけど，うっかりが怖い。テストデータで x, y に長さ10000のデータを指定して結果が合わないことに気づくという，大失態。
特に問題だったのは， variance を計算する式。他の所も羮に懲りて膾を吹いたところがあるし，もっと大きなデータが与えられたらアウトの所も残っているかも。

==========#

using StatsBase, Rmath

function moodtest(x, y; alternative = "twosided") # or "less" or "greater"
  nx = length(x)
  ny = length(y)
  n = nx + ny
  n >= 3 || error("not enough observations")
  expectation = float(nx) * (n ^ 2 - 1) / 12
  variance = float(nx) * ny * (n + 1) * (n + 2) * (n - 2) / 180
  z = vcat(x, y);
  u = sort!(unique(z))
  a = tabulate2(x, u)
  t = tabulate2(z, u)
  p = cumsum(((1:length(z)) .- (n + 1) / 2) .^ 2)
  ft = float(t)
  variance -= (nx * ny) / (180 * n * (n - 1)) * 
               sum(ft .* (ft .^ 2 .- 1) .* (ft .^ 2 .- 4 .+ 15 .* (n .- ft) .^ 2))
  statistic = sum(a .* diff(vcat(0, p[cumsum(t)])) ./ ft)
  z = (statistic - expectation) / sqrt(variance)
  p = pnorm(z)
  if alternative == "less"
    pvalue = p
  elseif alternative == "greater"
    pvalue = 1 - p
  else # twosided
    pvalue = 2 * min(p, 1 - p)
  end
  println("S = $statistic,  E(S) = $expectation,  V(S) = $variance")
  println("z = $z,  p value = $pvalue")
  Dict(:S => statistic,  :E => expectation,  :V => variance,
       :z => z, :pvalue => pvalue)
end

function tabulate2(x, u)
  y = vcat([indexin([x0], u) for x0 in x]...)
  ret = zeros(Int, length(u))
  [ret[i] += 1 for i in y]
  ret
end

x = [111, 107, 100, 99, 102, 106, 109, 108, 104, 99, 101, 96, 97, 102, 107, 113, 116, 113, 110, 98];
y = [107, 108, 106, 98, 105, 103, 110, 105, 104, 100, 96, 108, 103, 104, 114, 114, 113, 108, 106, 99];

moodtest(x, y)
# S = 3051.3333333333335,  E(S) = 2665.0,  V(S) = 138758.6324786325
# z = 1.0371275614960966,  p value = 0.2996764118570592

moodtest(x, y, alternative="less")
# S = 3051.3333333333335,  E(S) = 2665.0,  V(S) = 138758.6324786325
# z = 1.0371275614960966,  p value = 0.8501617940714704

moodtest(x, y, alternative="greater")
# S = 3051.3333333333335,  E(S) = 2665.0,  V(S) = 138758.6324786325
# z = 1.0371275614960966,  p value = 0.1498382059285296

x = [111, 107, 100, 99, 102, 106, 109, 108, 104];
y = [101, 96, 97, 113, 116, 110, 98, 105, 103, 114];

moodtest(x, y)
# S = 137.0,  E(S) = 270.0,  V(S) = 3570.0
# z = -2.225960907290223,  p value = 0.026016800044908195

moodtest(x, y, alternative="less")
# S = 137.0,  E(S) = 270.0,  V(S) = 3570.0
# z = -2.225960907290223,  p value = 0.013008400022454097

moodtest(x, y, alternative="greater")
# S = 137.0,  E(S) = 270.0,  V(S) = 3570.0
# z = -2.225960907290223,  p value = 0.9869915999775459