Julia に翻訳--122 平均値の多重比較，ライアン法，テューキー法

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

平均値の多重比較（ライアンの方法，テューキーの方法）
http://aoki2.si.gunma-u.ac.jp/R/m_multi_comp.html

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

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

Rmath の qtukey() にバグがある。

qtukey(q::Number, nmeans::Number, df::Number, nranges::Number=1.0,
    lower_tail::Bool=true, log_p::Bool=false) =
    ccall((:qtukey ,libRmath), Float64,
    (Float64, Float64, Float64, Float64, Int32, Int32),
    p, nranges, nmeans, df, lower_tail, log_p)
end

# 最初の引数 'q::Number' は q ではなく，p である。すなわち，'p::Number'
# しかし，それを直しても lower_tail = false を指定すると NaN を返す。
# 回避策としては，
# qtukey(p, nmeans, df, false)
# を
# qtukey(1-p, nmeans, df)
# とする。

==========#

using Rmath, Roots, Printf

function mmulticomp(n, me, s,; alpha = 0.05, method = "ryan") # or method = "tukey"
    k = length(n)
    all(k == length(me) && k == length(s) && all(n .> 0) && all(floor.(n) .== n) && all(s .> 0)) || error("parameter error")
    o = sortperm(me)
    sn = n[o]
    sm = me[o]
    ss = s[o]
    nt = sum(sn)
    mt = sum(sn .* sm) / nt
    dfw = nt - k
    vw = sum((sn .- 1) .* ss .^ 2) / dfw
    numsignificant = nsn = 0
    nss = zeros(Int, k * (k - 1) ÷ 2)
    nsb = zeros(Int, k * (k - 1) ÷ 2)
    for m = k:-1:2
        for small = 1:k - m + 1
            big = small + m - 1
            if check(small, big, nsn, nss, nsb)
                t0 = (sm[big] - sm[small]) / sqrt(vw * (1 / sn[big] + 1 / sn[small]))
                if method == "ryan"
                    P = pt(t0, dfw, false) * 2
                    nominalalpha = 2 * alpha / (k * (m - 1))
                    result = P <= nominalalpha
                else
                    t0 *= sqrt(2)
                    q = (qtukey(0.95, k, dfw) + qtukey(0.95, m, dfw)) / 2
                    WSD = q * sqrt(vw * (1 / sn[big] + 1 / sn[small]) / 2)
                    P = find_zero(p -> ((qtukey(1 - p, k, dfw) +
                                 qtukey(1 - p, m, dfw)) / 2 - t0), [0, 1], Bisection())
                    result = sm[big] - sm[small] >= WSD
                end
                if result
                    numsignificant = 1
                    @printf("mean[ % 2i] = %7.5f vs mean[ % 2i] = %7.5f : diff = % 7.5f, ",
                        o[small], me[small], o[big], me[big], me[big] - me[small])
                    if method == "ryan"
                        @printf("t = %7.5f : P = %7.5f, alpha' = %7.5f\n", t0, P, nominalalpha)
                    else
                        @printf("WSD = %7.5f : t = %7.5f : P = %7.5f\n", WSD, t0, P)
                    end
                else
                    nsn += 1
                    nss[nsn] = small
                    nsb[nsn] = big
                end
            end
        end
    end
    if numsignificant == 0
        print("Not significant at all")
    end
end

function check(s, b, nsn, nss, nsb)
    if nsn > 1
        for i in 1:nsn
            if (nss[i] <= s <= nsb[i]) && (nss[i] <= b <= nsb[i])
                return false
            end
        end
    end
    return true
end

mmulticomp([8, 11, 22, 6], [135.83, 160.49, 178.35, 188.06], [19.59, 12.28, 15.01, 9.81], method = "ryan")
# mean[  1] = 135.83000 vs mean[  4] = 188.06000 : diff =  52.23000, t = 6.53866 : P = 0.00000, alpha' = 0.00833
# mean[  1] = 135.83000 vs mean[  3] = 178.35000 : diff =  42.52000, t = 6.96308 : P = 0.00000, alpha' = 0.01250
# mean[  2] = 160.49000 vs mean[  4] = 188.06000 : diff =  27.57000, t = 3.67279 : P = 0.00066, alpha' = 0.01250
# mean[  1] = 135.83000 vs mean[  2] = 160.49000 : diff =  24.66000, t = 3.58814 : P = 0.00085, alpha' = 0.02500
# mean[  2] = 160.49000 vs mean[  3] = 178.35000 : diff =  17.86000, t = 3.26998 : P = 0.00212, alpha' = 0.02500

mmulticomp([8, 11, 22, 6], [135.83, 160.49, 178.35, 188.06], [19.59, 12.28, 15.01, 9.81], method = "tukey")
# mean[  1] = 135.83000 vs mean[  4] = 188.06000 : diff =  52.23000, WSD = 21.34697 : t = 9.24706 : P = 0.00000
# mean[  1] = 135.83000 vs mean[  3] = 178.35000 : diff =  42.52000, WSD = 15.57114 : t = 9.84728 : P = 0.00000
# mean[  2] = 160.49000 vs mean[  4] = 188.06000 : diff =  27.57000, WSD = 19.14118 : t = 5.19411 : P = 0.00258
# mean[  1] = 135.83000 vs mean[  2] = 160.49000 : diff =  24.66000, WSD = 16.11328 : t = 5.07440 : P = 0.00195
# mean[  2] = 160.49000 vs mean[  3] = 178.35000 : diff =  17.86000, WSD = 12.80554 : t = 4.62444 : P = 0.00482

Julia に翻訳--121 13日の金曜日

次の 13 日の金曜日はいつか？
https://blog.goo.ne.jp/r-de-r/e/f806d8152abc1740ad261335aa4618cc

python だと

def Zeller(y, m, d):
    w = ['Sat', 'Sun', 'Mon', 'Tue', 'Wed', 'Thu', 'Fri']
    if m == 1 or m ==2:
        m += 12
        y -= 1
    C = y // 100
    Y = y % 100
    return w[(d + 26*(m+1)//10 + Y + Y//4 -2*C + C // 4) % 7]

for y in range(2021, 2030):
    for m in range(1, 13):
        if Zeller(y, m, 13) == 'Fri':
            print(y, m)

であった。

Julia でも，移植は簡単だろうと，安易に以下のようにすると，間違った答えが出てしまう。
Python の // は Julia では ÷ なのは間違いないのだが。

function zeller(y, m, d)
    w = ["Sat", "Sun", "Mon", "Tue", "Wed", "Thu", "Fri"]
    if m == 1 | m ==2
        m += 12
        y -= 1
    end
    C = floor(Int, y / 100)
    Y = y % 100
    return w[((d + 26*(m+1) ÷ 10 + Y + Y ÷ 4 -2*C + C ÷ 4) % 7) + 1]
end

for y = 2021:2029
    for m = 1:12
        if zeller(y, m, 13) == "Fri"
            println("$y, $m")
        end
    end
end

出力結果

2021, 8
2022, 5
2023, 2
2023, 10
2024, 1
2024, 9
2024, 12
2025, 6
2026, 3
2026, 11
2027, 8
2028, 2
2028, 10
2029, 4
2029, 7

デバッグ宜しく（答えを見るなら，スクロールしてね）

 

 

 

 

 

 

 

 

 

 

 

バグは
if m == 1 | m ==2
のところ。

ここは，ビット演算子ではなく論理演算子でなくてはならない。
つまり，
if m == 1 || m ==2

これを修正して，以下のようにすると，正しい答えが得られる。

function zeller(y, m, d)
    w = ["Sat", "Sun", "Mon", "Tue", "Wed", "Thu", "Fri"]
    if m == 1 || m ==2
        m += 12
        y -= 1
    end
    C = y ÷ 100
    Y = y % 100
    return w[((d + 26*(m+1) ÷ 10 + Y + Y ÷ 4 -2*C + C ÷ 4) % 7) + 1]
end

for y = 2021:2029
    for m = 1:12
        if zeller(y, m, 13) == "Fri"
            println("$y, $m")
        end
    end
end

出力結果

2021, 8
2022, 5
2023, 1
2023, 10
2024, 9
2024, 12
2025, 6
2026, 2
2026, 3
2026, 11
2027, 8
2028, 10
2029, 4
2029, 7

 

Julia に翻訳--120 比率の多重比較，ライアン法，テューキー法

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

比率の多重比較（ライアンの方法，テューキーの方法）
http://aoki2.si.gunma-u.ac.jp/R/p_multi_comp.html

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

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

Rmath の qtukey() にバグがある。

qtukey(q::Number, nmeans::Number, df::Number, nranges::Number=1.0,
    lower_tail::Bool=true, log_p::Bool=false) =
    ccall((:qtukey ,libRmath), Float64,
    (Float64, Float64, Float64, Float64, Int32, Int32),
    p, nranges, nmeans, df, lower_tail, log_p)
end

# 最初の引数 'q::Number' は q ではなく，p である。すなわち，'p::Number'
# しかし，それを直しても lower_tail = false を指定すると NaN を返す。
# 回避策としては，
# qtukey(p, nmeans, df, false)
# を
# qtukey(1-p, nmeans, df)
# とする。

==========#

using Rmath, Printf

function pmulticomp(n, r; alpha = 0.05, method = "ryan") # or method = "tukey"
    k = length(n)
    all(k == length(r) && all(n .> 0) && all(r .>= 0) && all(r .<= n) &&
        all(floor.(n) .== n) && all(floor.(r) .== r)) || error("parameter error")
    o = sortperm(r ./ n)
    sr = r[o]
    sn = n[o]
    srsn = sr ./ sn
    numsignificant = nsn = 0
    nss = zeros(Int, k * (k - 1) ÷ 2)
    nsb = zeros(Int, k * (k - 1) ÷ 2)
    # m = k; small = 1
    q = NaN # dummy これがないと動かない
    for m = k:-1:2
        for small = 1:k - m + 1
            big = small + m - 1
            if check(small, big, nsn, nss, nsb)
                prop = sum(sr[small:big]) / sum(sn[small:big])
                se = sqrt(prop * (1 - prop) * (1 / sn[small] + 1 / sn[big]))
                diff = srsn[big] - srsn[small]
                if method == "ryan"
                    nominalalpha = 2 * alpha / (k * (m - 1))
                    rd = se * qnorm(nominalalpha / 2, false)
                    z = se == 0 ? 0 : diff / se
                    p = pnorm(z, false) * 2
                    result = p <= nominalalpha
                else
                    if m == k
                        q = qtukey(1-alpha, k, Inf)
                        qq = q
                    else
                        qq = (q + qtukey(1-alpha, m, Inf)) / 2
                    end
                    WSD = se * qq / sqrt(2)
                    result = diff != 0 && diff >= WSD
                end
                if result
                    numsignificant = 1
                    @printf("p[%2i] = %7.5f vs p[%2i] = %7.5f : diff = %7.5f, ",
                        o[small], srsn[small], o[big], srsn[big], diff)
                    if method == "ryan"
                        @printf("RD = %7.5f : P = %7.5f, alpha' = %7.5f\n", rd, p, nominalalpha)
                    else
                        @printf("WSD = %7.5f\n", WSD)
                    end
                else
                    nsn += 1
                    nss[nsn] = small
                    nsb[nsn] = big
                end
            end
        end
    end
    if numsignificant == 0
        println("Not significant at all")
    end
end

function check(s, b, nsn, nss, nsb)
    if nsn > 1
        for i in 1:nsn
            if (nss[i] <= s <= nsb[i]) && (nss[i] <= b <= nsb[i])
                return false
            end
        end
    end
    return true
end

n = [30, 35, 47, 21, 45];
r = [2, 4, 14, 13, 39];

pmulticomp(n, r)
# p[ 1] = 0.06667 vs p[ 5] = 0.86667 : diff = 0.80000, RD = 0.32472 : P = 0.00000, alpha' = 0.00500
# p[ 1] = 0.06667 vs p[ 4] = 0.61905 : diff = 0.55238, RD = 0.33341 : P = 0.00001, alpha' = 0.00667
# p[ 2] = 0.11429 vs p[ 5] = 0.86667 : diff = 0.75238, RD = 0.30528 : P = 0.00000, alpha' = 0.00667
# p[ 1] = 0.06667 vs p[ 3] = 0.29787 : diff = 0.23121, RD = 0.23054 : P = 0.00979, alpha' = 0.01000
# p[ 2] = 0.11429 vs p[ 4] = 0.61905 : diff = 0.50476, RD = 0.32612 : P = 0.00007, alpha' = 0.01000
# p[ 3] = 0.29787 vs p[ 5] = 0.86667 : diff = 0.56879, RD = 0.26479 : P = 0.00000, alpha' = 0.01000
# p[ 3] = 0.29787 vs p[ 4] = 0.61905 : diff = 0.32118, RD = 0.29877 : P = 0.01239, alpha' = 0.02000

pmulticomp(n, r, method="tukey")
# p[ 1] = 0.06667 vs p[ 5] = 0.86667 : diff = 0.80000, WSD = 0.31555
# p[ 1] = 0.06667 vs p[ 4] = 0.61905 : diff = 0.55238, WSD = 0.32546
# p[ 2] = 0.11429 vs p[ 5] = 0.86667 : diff = 0.75238, WSD = 0.29800
# p[ 1] = 0.06667 vs p[ 3] = 0.29787 : diff = 0.23121, WSD = 0.22695
# p[ 2] = 0.11429 vs p[ 4] = 0.61905 : diff = 0.50476, WSD = 0.32104
# p[ 3] = 0.29787 vs p[ 5] = 0.86667 : diff = 0.56879, WSD = 0.26067
# p[ 3] = 0.29787 vs p[ 4] = 0.61905 : diff = 0.32118, WSD = 0.30102

Julia に翻訳--119 一対比較データの双対尺度法，西里

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

一対比較データの双対尺度法（西里）
http://aoki2.si.gunma-u.ac.jp/R/pc.dual.html

ファイル名: pcdual.jl　　関数名: pcdual, plotpcdual

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

考えて見れば当たり前なんだが，annotate! にもドットを付けることができるのだ。
その際，必要ならば text にもドットを付けることを忘れずに。

==========#

using Rmath, Combinatorics, LinearAlgebra, Printf, Plots

function pcdual(F, onetwo = true)
    N, C = size(F)
    if onetwo
        F[F .== 2] .= - 1
    end
    reshape(F, 14, :)
    n = floor(Int, 1 + sqrt(1 + 8 * C) / 2)
    x = combinations(1:n, 2)
    nc = binomial(8, 2)
    A = zeros(nc, n)
    for (i, (j, k)) in  enumerate(x)
        A[i, j] = 1
        A[i, k] = -1
    end
    E = F * A
    Hn = transpose(E) * E ./ (N * n * (n - 1) ^ 2)
    values, vectors = eigen(Hn, sortby=x-> -x)　
    ne = size(Hn, 1) - 1
    eta2 = values[1:ne]
    eta = sqrt.(eta2)
    contribution = eta2 ./ sum(values[1:ne]) * 100
    cumcont = cumsum(contribution)
    W = vectors[:, 1:ne]
    colscore = W .* sqrt(n)
    colscore2 = colscore .* transpose(eta)
    rowscore2 = E * (W ./ (sqrt(n) * (n - 1)))
    rowscore = rowscore2 ./ transpose(eta)
    printpcdual(N, n, ne, eta2, eta, contribution, cumcont,
                rowscore, colscore, rowscore2, colscore2)
    Dict(:eta2 => eta2, :eta => eta, :contribution => contribution,
         :cumcont => cumcont, :rowscore => rowscore,
         :colscore => colscore, :rowscoreweighted => rowscore2,
         :colscoreweighted => colscore2)
end

function printrow(name, x)
    @printf("%9s", name)
    for y in x
        @printf("%8.3f", y)
    end
    println()
end

function printheader(ne)
    @printf("%9s", "")
    for i = 1:ne
        @printf("%8s", "Axis-" * string(i))
    end
    println()
end

function printpcdual(N, n, ne, eta2, eta, contribution, cumcont,
                     rowscore, colscore, rowscore2, colscore2)
    printheader(ne)
    printrow("eta sq.",   eta2)
    printrow("cor.",  eta)
    printrow("cont", contribution)
    printrow("cum.cont.",   cumcont)
    println("\nrow score")
    printheader(ne)
    for i = 1:N
        printrow("Row-" * string(i), rowscore[i, :])
    end
    println("\ncolumn score")
    printheader(ne)
    for i = 1:n
        printrow("Col-" * string(i), colscore[i, :])
    end
    println("\nweighted row score")
    printheader(ne)
    for i = 1:N
        printrow("Row-" * string(i), rowscore2[i, :])
    end
    println("\nweighted column score")
    printheader(ne)
    for i = 1:n
        printrow("Col-" * string(i), colscore2[i, :])
    end
end

function plotpcdual(results; weighted=true, ax1=1, ax2=2)
    pyplot()
    if weighted
        row, col = results[:rowscore], results[:colscore]
    else
        row, col = results[:rowscoreweighted], results[:colscoreweighted]
    end
    rownames = [" R" * string(i) for i = 1:size(row, 1)]
    colnames = [" C" * string(i) for i = 1:size(col, 1)]
    plt = scatter(row[:, ax1], row[:, ax2], grid=false, tick_direction=:out,
        color=:red, xlabel="Axis-"*string(ax1), ylabel="Axis-"*string(ax2),
        aspect_ratio=1, size=(600, 600), label="")
    annotate!.(row[:, ax1], row[:, ax2], text.(rownames, 8, :red, :left));
    scatter!(col[:, ax1], col[:, ax2], grid=false, tick_direction=:out,
            color=:blue, label="")
    annotate!.(col[:, ax1], col[:, ax2], text.(colnames, 8, :blue, :left))
    display(plt)
end

F = [   1 1 2 1 1 2 1 2 2 2 2 2 2 2 1 1 2 1 1 1 2 1 1 2 1 2 1 2
        2 2 2 1 2 1 2 1 2 1 2 1 2 2 1 1 1 2 1 1 1 2 2 2 2 1 2 2
        1 1 1 1 1 2 1 1 1 1 1 2 1 1 1 1 2 1 1 1 2 1 2 2 2 2 1 1
        2 1 2 1 1 1 2 1 1 1 1 1 2 2 1 2 2 2 1 1 1 2 2 2 2 2 2 2
        2 2 2 1 2 1 2 2 2 1 2 2 2 2 1 2 1 2 1 1 1 1 2 2 2 1 2 2
        1 1 1 1 1 1 1 2 2 1 2 2 2 2 1 2 2 2 1 1 1 1 2 2 2 2 2 1
        1 1 1 1 1 2 1 1 2 1 1 2 1 2 1 1 2 1 1 1 2 1 2 2 2 2 2 1
        1 1 1 1 1 2 1 1 2 1 2 2 1 2 1 2 2 1 1 2 2 1 2 2 1 2 1 1
        1 2 2 1 1 2 1 2 2 1 1 2 2 1 1 1 2 1 1 1 2 1 2 2 2 2 2 1
        1 2 1 1 2 2 2 2 2 1 2 2 2 1 1 1 1 1 1 2 2 2 2 2 2 1 1 2
        1 2 1 1 1 1 1 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2
        2 2 2 2 1 2 2 1 2 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 2 2 1
        1 2 1 1 2 1 2 2 2 2 2 2 2 1 1 1 1 1 1 2 1 2 2 2 2 1 1 2
        2 2 2 2 1 2 1 2 1 1 1 1 1 1 1 1 1 1 2 1 2 1 1 2 1 2 1 1 ];

results = pcdual(F)
#==
           Axis-1  Axis-2  Axis-3  Axis-4  Axis-5  Axis-6  Axis-7
eta sq.     0.141   0.110   0.065   0.055   0.028   0.013   0.006
cor.        0.376   0.331   0.255   0.235   0.169   0.114   0.075
cont       33.719  26.241  15.591  13.176   6.800   3.122   1.352
cum.cont.  33.719  59.960  75.550  88.726  95.526  98.648 100.000

row score
        Axis-1  Axis-2  Axis-3  Axis-4  Axis-5  Axis-6  Axis-7
 Row-1   1.101  -0.417  -0.436  -0.156   1.989   0.483   1.119
 Row-2   0.334   1.473   1.417  -0.449  -0.520   0.458  -0.510
 Row-3   1.211  -1.160   0.105   0.158  -1.527  -0.377  -0.818
 Row-4   0.418   0.507   2.145   0.783  -0.898  -1.179   0.386
 Row-5   0.669   1.481   0.872  -0.332   0.255   1.894   0.030
 Row-6   1.498   0.195   0.345   1.129   0.984  -0.226  -0.539
 Row-7   1.422  -0.968   0.693   0.294  -0.299  -0.354   0.098
 Row-8   1.086  -0.927  -0.518   1.193  -0.501   2.117  -1.057
 Row-9   1.528  -0.366   0.095  -1.011  -0.122  -0.241   2.161
Row-10   0.808   1.149  -1.173  -0.618  -1.593   0.316   1.124
Row-11   1.296   0.613  -0.546  -0.939   1.037  -1.691  -1.601
Row-12  -0.120  -0.651   1.513  -1.906   0.797   0.856  -0.689
Row-13   0.678   1.453  -1.257  -0.585  -0.387  -0.071  -0.987
Row-14  -0.024  -1.289  -0.098  -1.995  -0.922   0.160  -0.555

column score
        Axis-1  Axis-2  Axis-3  Axis-4  Axis-5  Axis-6  Axis-7
 Col-1   1.042  -0.412  -0.488   1.244   0.373  -1.364  -1.400
 Col-2  -0.923  -0.949   1.504  -0.140  -1.563  -0.037  -0.722
 Col-3   0.506   0.361  -1.017  -2.246  -0.468  -0.411  -0.382
 Col-4   0.796   0.248   1.388  -0.443   1.567   1.240  -0.435
 Col-5  -1.950  -0.720  -0.721  -0.077   1.403  -0.293   0.317
 Col-6  -0.389   1.078  -1.113   0.972  -0.775   1.644  -0.449
 Col-7   1.091  -1.380  -0.361   0.326  -0.396   0.429   1.824
 Col-8  -0.173   1.776   0.807   0.365  -0.140  -1.209   1.246

weighted row score
        Axis-1  Axis-2  Axis-3  Axis-4  Axis-5  Axis-6  Axis-7
 Row-1   0.414  -0.138  -0.111  -0.037   0.335   0.055   0.084
 Row-2   0.126   0.488   0.362  -0.105  -0.088   0.052  -0.038
 Row-3   0.455  -0.384   0.027   0.037  -0.258  -0.043  -0.062
 Row-4   0.157   0.168   0.548   0.184  -0.151  -0.135   0.029
 Row-5   0.251   0.491   0.223  -0.078   0.043   0.216   0.002
 Row-6   0.563   0.065   0.088   0.265   0.166  -0.026  -0.041
 Row-7   0.534  -0.321   0.177   0.069  -0.050  -0.040   0.007
 Row-8   0.408  -0.307  -0.132   0.280  -0.084   0.242  -0.079
 Row-9   0.574  -0.121   0.024  -0.237  -0.021  -0.028   0.162
Row-10   0.304   0.381  -0.300  -0.145  -0.269   0.036   0.085
Row-11   0.487   0.203  -0.139  -0.220   0.175  -0.193  -0.120
Row-12  -0.045  -0.216   0.386  -0.448   0.134   0.098  -0.052
Row-13   0.255   0.482  -0.321  -0.137  -0.065  -0.008  -0.074
Row-14  -0.009  -0.427  -0.025  -0.468  -0.155   0.018  -0.042

weighted column score
        Axis-1  Axis-2  Axis-3  Axis-4  Axis-5  Axis-6  Axis-7
 Col-1   0.391  -0.137  -0.125   0.292   0.063  -0.156  -0.105
 Col-2  -0.347  -0.315   0.384  -0.033  -0.264  -0.004  -0.054
 Col-3   0.190   0.119  -0.260  -0.527  -0.079  -0.047  -0.029
 Col-4   0.299   0.082   0.355  -0.104   0.264   0.142  -0.033
 Col-5  -0.732  -0.239  -0.184  -0.018   0.237  -0.033   0.024
 Col-6  -0.146   0.357  -0.284   0.228  -0.131   0.188  -0.034
 Col-7   0.410  -0.457  -0.092   0.077  -0.067   0.049   0.137
 Col-8  -0.065   0.588   0.206   0.086  -0.024  -0.138   0.094
==#

plotpcdual(results)

Julia に翻訳--118 AHP，Analytic Hierachy Process

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

AHP (Analytic Hierachy Process)
http://aoki2.si.gunma-u.ac.jp/R/AHP.html

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

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

内包表記と for ループ どちらがよいだろうか？

==========#

using LinearAlgebra, Rmath, Plots

function ahp(x, y)
    n = items(length(x))
    weightx = makematrix(x)
    nitemsy = items(size(y, 1))
    weighty = zeros(nitemsy, n)
    for i = 1:n
        weighty[:, i] = makematrix(y[:, i])
    end
    score = weighty * weightx
    sortedscore = sort(score)
    printahp(weightx, weighty, score, sortedscore)
    plotahp(score, nitemsy)
    Dict(:weightx => weightx, :weighty => weighty,
         :score => score, :sortedscore => sortedscore)
end

function items(n)
    retval = (1 + sqrt(1 + 8 * n)) / 2
    retval != floor(retval) ? NaN : Int(retval)
end

function makematrix(x)
    n = items(length(x))
    !isnan(n) || error("size error")
    mat = Array{Float64}(I, n, n)
    lowertri = [i > j for i = 1:n, j = 1:n]
    mat[lowertri] = x
    mat = transpose(mat) + mat
    uppertri = [i < j for i = 1:n, j = 1:n]
    mat[uppertri] = 1 ./ mat[uppertri]
    [mat[i, i] = 1 for i in 1:n]
    values, vectors = eigen(mat)
    val = real(values[n])
    vec = real(vectors[:, n])
    weight = vec / sum(vec)
    ci = (val - n) / (n - 1)
    cr = ci / [0, 0, 0.58, 0.9, 1.12, 1.24, 1.32, 1.41, 1.45, 1.49, 1.51, 1.53][n]
    if ci > 0.1 || cr > 0.1
        println("\nCI = $ci, CR = $cr")
        println(mat)
        W = weight ./ weight'
        println(W)
        println(mat - W)
    end
    return weight
end

function printahp(weightx, weighty, score, sortedscore, digits = 5)
    println("\n評価基準の重み")
    println(round.(weightx, digits = digits))
    println("\n代替案の評価結果")
    println(round.(weighty, digits = digits))
    println("\nスコア")
    println(round.(score, digits = digits))
    println("\nソートされたスコア")
    println(round.(sortedscore, digits = digits))
end

function plotahp(score, nc)
    pyplot(size=(600, 60))
    plt = scatter(score, repeat([25], nc), grid=false, tick_direction=:out,
            yshowaxis=false, ylims=(0, 60), yticks=false, label="")
    for i = 1:nc
        annotate!(score[i], 50, text(string(i), 10, :red), label="")
    end
    display(plt)
end

x = [1 / 3, 1 / 5, 1 / 7, 1 / 5, 1 / 7, 1 / 3];
y = hcat([1 / 2, 1 / 3, 1 / 2], [5, 2, 1 / 7], [1 / 3, 1 / 2, 2], [2, 2, 1]);

ahp(x, y)
# CI = 0.05947572619908503, CR = 0.10254435551566386
# [1.0 0.2 0.5; 5.0 1.0 7.0; 2.0 0.14285714285714285 1.0]
# [1.0 0.14189834119703848 0.7047298732064893; 7.0472987320648866 1.0 4.96644194189634; 1.4189834119703837 0.20135139234471133 1.0]
# [0.0 0.05810165880296153 -0.20472987320648928; -2.0472987320648866 0.0 2.03355805810366; 0.5810165880296163 -0.058494249487568484 0.0]

# 評価基準の重み
# [0.54374, 0.31093, 0.09745, 0.04789]

# 代替案の評価結果
# [0.53961 0.10564 0.53961 0.2; 0.29696 0.74446 0.16342 0.4; 0.16342 0.1499 0.29696 0.4]

# スコア
# [0.38842, 0.42802, 0.18356]

# ソートされたスコア
# [0.18356, 0.38842, 0.42802]