Julia に翻訳--170 2×2 分割表のフィッシャーの正確確率検定

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

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

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

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

==========#

using SpecialFunctions, Printf

function fisher2x2(a, b, c, d; alternative="twosided")
    method = "Fisher's Exact Test for Count Data"
    a0 = a
    e, f, g, h, n = a + b, c + d, a + c, b + d, a + b + c + d
    x0, y0, z0 = stats(a, e, f, g, h, n)
    mx = min(e, g)
    mi = max(0, e + g - n)
    len = mx - mi + 1
    a1 = mi:mx
    b1 = e .- a1
    c1 = g .- a1
    d1 = f .- c1
    vchisq     = zeros(len);
    vprob      = similar(vchisq);
    voddsratio = similar(vchisq);
    for i =1:len
        vchisq[i], vprob[i], voddsratio[i] = stats(mi + i - 1, e, f, g, h, n)
    end
    if alternative == "less"
        mpearson = a1 .<= a0
        mfisher  = a1 .<= a0
        mor      = a1 .<= a0
    elseif alternative == "greater"
        # x = y = z = sum(vprob[a1 .>= a0])
        mpearson = a1 .>= a0
        mfisher  = a1 .>= a0
        mor      = a1 .>= a0
    else
        alternative = "twosided"
        mpearson = vchisq     .>= x0
        mfisher  = vprob      .<= y0
        mor      = voddsratio .>= z0
    end
    ppearson = sum(vprob[mpearson])
    pfisher  = sum(vprob[mfisher])
    por      = sum(vprob[mor])
    cumprob1 = cumsum(vprob)
    cumprob2 = reverse(cumsum(reverse(vprob)))
    println(method)
    @printf("%3s  %3s  %3s  %3s  %12s %s  %12s %s  %12s %s  %12s  %12s\n",
            "a", "b", "c", "d", "Chi square", " ", "Probability", " ", "Odds Ratio", " ",
            "Cum. Prob1", "Cum. Prob2")
    for i = 1:len
        @printf("%3d  %3d  %3d  %3d  %12.6f %s  %12.6f %s  %12.6f %s  %12.6f  %12.6f\n",
                a1[i], b1[i], c1[i], d1[i], vchisq[i], mark(mpearson[i]),
                vprob[i], mark(mfisher[i]), voddsratio[i], mark(mor[i]),
                cumprob1[i], cumprob2[i])
    end
    println("\nalternative = $alternative")
    println("p(Pearson) = $ppearson")
    println("p(Fisher) = $pfisher")
    println("p(Odds Ratio) = $por")
    Dict(:alternative => alternative, :pPearson => ppearson, :pFisher => pfisher, :pOddsRatio => por,
        :vchisq => vchisq, :oddsratio => voddsratio, :probability => vprob)
end

function oddsratio(a, b, c, d)
    a * b * c * d != 0 || ((a, b, c, d) = (a, b, c, d) .+ 0.5)
    OR = a * d / (b * c)
    max(OR, 1 / OR)
end

mark(x) = x ? "*" : " "

lbinomial(n, k) = logfactorial(n) - logfactorial(k) - logfactorial(n - k)

function stats(i, e, f, g, h, n)
    a, b, c = i, e - i, g - i
    d = f - c
    vchisq = n * (a * d - b * c) ^ 2 / (e * f * g * h)
    or = oddsratio(a, b, c, d)
    vprob = exp(lbinomial(e, a) + lbinomial(f, c) - lbinomial(n, g))
    vchisq, vprob, or
end

a, b, c, d = 10, 13, 16, 61
fisher2x2(a, b, c, d)
#=====
Fisher's Exact Test for Count Data
  a    b    c    d    Chi square     Probability      Odds Ratio      Cum. Prob1    Cum. Prob2
  0   23   26   51     10.494910 *      0.000332 *     24.184466 *      0.000332      1.000000
  1   22   25   52      7.278386 *      0.003815 *     10.576923 *      0.004147      0.999668
  2   21   24   53      4.648818        0.019796 *      4.754717 *      0.023943      0.995853
  3   20   23   54      2.606207        0.061587        2.839506        0.085530      0.976057
  4   19   22   55      1.150553        0.128773        1.900000        0.214302      0.914470
  5   18   21   56      0.281857        0.192239        1.350000        0.406541      0.785698
  6   17   20   57      0.000117        0.212475        1.005882        0.619016      0.593459
  7   16   19   58      0.305335        0.177934        1.335526        0.796950      0.380984
  8   15   18   59      1.197510        0.114602        1.748148        0.911552      0.203050
  9   14   17   60      2.676642        0.057301        2.268908        0.968853      0.088448
 10   13   16   61      4.742731 *      0.022357 *      2.932692 *      0.991210      0.031147
 11   12   15   62      7.395777 *      0.006818 *      3.788889 *      0.998028      0.008790
 12   11   14   63     10.635780 *      0.001623 *      4.909091 *      0.999651      0.001972
 13   10   13   64     14.462741 *      0.000300 *      6.400000 *      0.999952      0.000349
 14    9   12   65     18.876658 *      0.000043 *      8.425926 *      0.999995      0.000048
 15    8   11   66     23.877533 *      0.000005 *     11.250000 *      1.000000      0.000005
 16    7   10   67     29.465364 *      0.000000 *     15.314286 *      1.000000      0.000000
 17    6    9   68     35.640153 *      0.000000 *     21.407407 *      1.000000      0.000000
 18    5    8   69     42.401899 *      0.000000 *     31.050000 *      1.000000      0.000000
 19    4    7   70     49.750602 *      0.000000 *     47.500000 *      1.000000      0.000000
 20    3    6   71     57.686262 *      0.000000 *     78.888889 *      1.000000      0.000000
 21    2    5   72     66.208879 *      0.000000 *    151.200000 *      1.000000      0.000000
 22    1    4   73     75.318454 *      0.000000 *    401.500000 *      1.000000      0.000000
 23    0    3   74     85.014985 *      0.000000 *   1000.428571 *      1.000000      0.000000

alternative = twosided
p(Pearson) = 0.035294133134962415
p(Fisher) = 0.05508990606358157
p(Odds Ratio) = 0.05508990606358157

Julia に翻訳--169 Box-Pierce 検定，Ljung-Box 検定

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

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

Box-Pierce 検定または Ljung-Box 検定
http://aoki2.si.gunma-u.ac.jp/Python/Box_test.pdf

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

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

Python --> Julia は，R --> Julia より面倒。
理由はいろいろあるが，配列のインデックスが 0 始まりか 1 始まりかの違い。
「C や Python のように 0 始まりが普通だろう」という人がいるが，FORTRAN 師匠は 1 始まりが当たり前だと仰っておられる。

==========#

using Statistics, Rmath

function boxtest(x; lag=1, method="Box-Pierce", fitdf=0)
    x = vec(x)
    ndims(x) == 1 || error("x != a vector")
    cor = acf(x, lag)
    n = length(x)
    df = lag - fitdf
    if method == "Box-Pierce"
        method = "Box-Pierce test"
        chisq = n * sum(cor .^ 2)
    else
        method = "Ljung-Box test"
        chisq = n * (n + 2) * sum(cor .^ 2 ./ (n - 1:-1:n - lag))
    end
    pvalue = pchisq(chisq, lag - fitdf, false)
    println("$method\nchisq. = $chisq,  df = $df,  p value = $pvalue")
    Dict(:chisq => chisq, :df => df, :pvalue => pvalue, :method => method)
end

acf(x, lag) = [acf1(x, i) for i = 1:lag]

function acf1(x, lag)
    n = length(x)
    (n >= 3 &&  n - lag >= 2 && lag >= 0) || error("invalid argument")
    y = x .- mean(x)
    sum(y[1:n - lag] .* y[lag+1: n]) / (var(x) * (n - 1))
end

x = [-0.56, -0.23, 1.559, 0.071, 0.129, 1.715, 0.461, -1.265, -0.687,
     -0.446, 1.224, 0.36, 0.401, 0.111, -0.556, 1.787, 0.498, -1.967,
     0.701, -0.473, -1.068, -0.218, -1.026, -0.729, -0.625, -1.687,
     0.838, 0.153, -1.138, 1.254, 0.426, -0.295, 0.895, 0.878, 0.822,
     0.689, 0.554, -0.062, -0.306, -0.38, -0.695, -0.208, -1.265,
     2.169, 1.208, -1.123, -0.403, -0.467, 0.78, -0.083, 0.253, -0.029,
     -0.043, 1.369, -0.226, 1.516, -1.549, 0.585, 0.124, 0.216, 0.38,
     -0.502, -0.333, -1.019, -1.072, 0.304, 0.448, 0.053, 0.922, 2.05,
     -0.491, -2.309, 1.006, -0.709, -0.688, 1.026, -0.285, -1.221,
     0.181, -0.139, 0.006, 0.385, -0.371, 0.644, -0.22, 0.332, 1.097,
     0.435, -0.326, 1.149, 0.994, 0.548, 0.239, -0.628, 1.361, -0.6,
     2.187, 1.533, -0.236, -1.026];

boxtest(x, lag=1)
# Box-Pierce test
# chisq. = 0.06539836770776197,  df = 1,  p value = 0.7981585211157581

boxtest(x, lag=1, method="Ljung-Box")
# Ljung-Box test
# chisq. = 0.067380136426179,  df = 1,  p value = 0.7951902022577569

boxtest(x, lag=3, fitdf=1)
# Box-Pierce test
# chisq. = 3.5535597565015773,  df = 2,  p value = 0.16918205789358604

boxtest(x, lag=3, method="Ljung-Box", fitdf=1)
# Ljung-Box test
# chisq. = 3.721733299128315,  df = 2,  p value = 0.15553777519526657