Julia に翻訳--051 級内相関係数

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

級内相関係数
http://aoki2.si.gunma-u.ac.jp/R/ic.html

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

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

ドット演算子の使いすぎに注意。

==========#

using Statistics

function intraclasscorrelation(dat)
    nr, nc = size(dat)
    nt = length(dat)
    Sw = (nc - 1) * sum(var(dat, dims=2))
    Mw = Sw / (nt - nr)
    m = mean(dat, dims=2)
    gm = mean(dat)
    Sb = nc * sum((m .- gm) .^ 2)
    Mb = Sb / (nr - 1)
    (Mb - Mw) / (Mb + (nc - 1) * Mw)
end

dat = [
  71  71;
  79  82;
 105  99;
 115 114;
  76  70;
  83  82;
 114 113;
  57  44;
 114 113;
  94  91;
  75  83;
  76  72
  ];

intraclasscorrelation(dat) # 0.9656287577888716

Julia に翻訳--050 属性相関係数

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

属性相関係数
http://aoki2.si.gunma-u.ac.jp/R/zokusei-soukan.html

ファイル名: phi.jl
関数名:    phi, expectation, chisq, df, chisqpvalue, contingency, cramer,
          likelihoodratio, likelihoodratiopvalue, ln, sumln

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

（お遊びで，）必要になる機能を一行で関数指定して，組み合わせて目的を達成する。

==========#

using Statistics, Distributions

# 分割表の期待値
expectation(mat) = (sum(mat, dims=2) * sum(mat, dims=1)) / sum(mat)

# χ2検定統計量
chisq(mat) = sum((mat .- expectation(mat)) .^ 2 ./ expectation(mat))

# χ2検定統計量が従うχ2分布の自由度
df(mat) = (size(mat, 1) - 1) * (size(mat, 2) - 1)

# χ2検定の P 値
chisqpvalue(mat) = ccdf(Chisq(df(mat)), chisq(mat))

# φ係数
phi(mat) = sqrt(chisq(mat)/sum(mat))

# コンティンジェンシー係数
contingency(mat) = sqrt(chisq(mat) / (sum(mat) + chisq(mat)))

# クラメール係数
cramer(mat) = phi(mat) / sqrt(min(size(mat, 1), size(mat, 2)) - 1)

# 対数尤度
ln(n) = n == 0 ? 0 : n .* log.(n)

# 対数尤度の和
sumln(n) = sum(filter(x -> !isnan(x), ln(n)))

# 対数尤度比
likelihoodratio(mat) = 2 * (sumln(mat) - sumln(sum(mat, dims=2))
                        - sumln(sum(mat, dims=1)) + ln(sum(mat)))

# 対数尤度比検定の P 値
likelihoodratiopvalue(mat) =  ccdf(Chisq(df(mat)), likelihoodratio(mat))

mat = [4 5 2 0;
       0 7 6 1;
       1 0 3 1]

expectation(mat)
# 3×4 Array{Float64,2}:
#  1.83333   4.4  4.03333  0.733333
#  2.33333   5.6  5.13333  0.933333
#  0.833333  2.0  1.83333  0.333333

chisq(mat)       # 11.344332939787485
df(mat)          # 6
chisqpvalue(mat) # 0.07829981227251993
phi(mat)         # 0.6149344935245131
contingency(mat) # 0.5238192896958154
cramer(mat)      # 0.4348243503566983
likelihoodratio(mat)       # 15.364591286599591
likelihoodratiopvalue(mat) # 0.01760288650305146

Julia に翻訳--049 ロバストな回帰直線

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

ロバストな回帰直線
http://aoki2.si.gunma-u.ac.jp/R/least.sum.abs.html

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

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

テストデータでは R と Julia で，数値的にはかなり異なる結果になっている。

==========#

using Statistics, Optim, Printf, Plots

function leastsumabs(x, y; n=1, p=1, sig=0.95,
                        col1=:black, col2=:red, xlab="", ylab="", legend=:best,
                        hist=false, FUNC=NelderMead) #FUNC=NelderMead, LBFGS or BFGS =#)
    intercept, slope = leastsumabs0(x, y, p, FUNC)
    if n == 1
        @printf("Intercept: %.7g\n", intercept)
        @printf("    Slope: %.7g\n", slope)
        results = Dict(:intercept => intercept, :slope => slope)
    else
        interceptvector = zeros(n)
        slopevector = zeros(n)
        suffix0 = collect(1:length(x))
        for i = 1:n
            suffix = rand(suffix0, length(x))
            interceptvector[i], slopevector[i] = leastsumabs0(x[suffix], y[suffix], p, FUNC)
        end
        alpha = (1-sig)/2
        LCLintercept, UCLintercept = quantile!(interceptvector, [alpha, 1-alpha])
        LCLslope, UCLslope = quantile!(slopevector, [alpha, 1-alpha])
        @printf("            Estimate      %6.3f%%   %6.3f%%\n", alpha, 1-alpha)
        @printf("Intercept:  %8.7g  [%8.7g, %8.7g]\n", intercept, LCLintercept, UCLintercept)
        @printf("Slope:      %8.7g  [%8.7g, %8.7g]\n", slope, LCLslope, UCLslope)
        results = Dict(:intercept => intercept, :slope => slope,
             :LCLintercept => LCLintercept, :UCLintercept => UCLintercept,
             :LCLslope => LCLslope, :UCLslope => UCLslope)
        !hist || simurationresults(interceptvector, slopevector)
    end
    plotleastsumabs(x, y, intercept, slope, col1, col2, xlab, ylab, legend)
    return results
end

function leastsumabs0(x, y, p, FUNC)
    evaluateerror(par) = sum(abs.(y .- par[2] .* x .- par[1]) .^ p)
    if p < 1
        println("p は 1 以上の値でなければなりません。p=1 として実行します。")
        p = 1
    end
    result = optimize(evaluateerror, zeros(2), FUNC())
    par = Optim.minimizer(result)
    return par
end

function simurationresults(interceptvector, slopevector)
    plthist1 = histogram(interceptvector, grid=false, tick_direction=:out,
                            alpha=0.2, label="intercept")
    plthist2 = histogram(slopevector, grid=false, tick_direction=:out,
                            alpha=0.2, label="slope")
    plthist0 = plot(plthist1, plthist2)
    display(plthist0)
end

function plotleastsumabs(x, y, intercept, slope, col1, col2, xlab, ylab, legend)
    abline(intercept, slope, col, label) =
        plot!([minx, maxx],
              [slope * minx + intercept, slope * maxx + intercept],
              color=col, label=label, legend=legend)
    function lm(x, y)
        b = cov(x, y)/var(x)
        return mean(y) - b * mean(x), b
    end
    minx, maxx = extrema(x)
    margin = (maxx - minx) * 0.05
    minx -= margin
    maxx += margin
    pyplot()
    plt = scatter(x, y, grid=false, tick_direction=:out,
                    xlabel=xlab, ylabel=ylab, color=col1, label="")
    abline(intercept, slope, col1, "least sum abs")
    a, b = lm(x, y)
    abline(a, b, col2, "linear regression")
    display(plt)
end

x = 1961:1993
y = [139.130, 140.293, 143.137, 147.350, 150.686, 144.317, 151.207,
    152.882, 156.867, 155.749, 157.735, 162.962, 159.036, 158.589, 149.213,
    148.725, 161.331, 161.363, 158.899, 142.862, 139.085, 162.026, 162.117,
    163.621, 152.982, 170.722, 162.175, 144.809, 167.581, 185.984, 176.457,
    134.477, 134.477];

leastsumabs(x, y)
# Intercept: -1267.525
#     Slope: 0.7212667
# Dict{Symbol,Float64} with 2 entries:
#   :intercept => -1267.52
#   :slope     => 0.721267

leastsumabs(x, y, n=1000, hist=true)
# Estimate       0.025%    0.975%
# Intercept:  -1267.525  [-1991.714, 655.0887]
# Slope:      0.7212667  [-0.2552111,  1.08665]
# Dict{Symbol,Float64} with 6 entries:
# :UCLintercept => 655.089
# :UCLslope     => 1.08665
# :intercept    => -1267.52
# :slope        => 0.721267
# :LCLslope     => -0.255211
# :LCLintercept => -1991.71

Julia に翻訳--048 抵抗直線

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

抵抗直線
http://aoki2.si.gunma-u.ac.jp/R/resistant-line.html

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

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

テンプレートがたまってくると，翻訳処理が楽になる。

==========#

using Statistics, Printf, Plots

function resistantline(x, y; noiteration=false, n=1, sig=0.95,
                        col1=:black, col2=:red, xlab="", ylab="", legend=:best,
                        hist=false)
    intercept, slope = resistantline0(x, y, noiteration)
    if n == 1
        @printf("Intercept: %.7g\n", intercept)
        @printf("    Slope: %.7g\n", slope)
        results = Dict(:intercept => intercept, :slope => slope)
    else
        interceptvector = zeros(n)
        slopevector = zeros(n)
        suffix0 = collect(1:length(x))
        for i = 1:n
            suffix = rand(suffix0, length(x))
            interceptvector[i], slopevector[i] = resistantline0(x[suffix], y[suffix], noiteration)
        end
        alpha = (1-sig)/2
        LCLintercept, UCLintercept = quantile!(interceptvector, [alpha, 1-alpha])
        LCLslope, UCLslope = quantile!(slopevector, [alpha, 1-alpha])
        @printf("            Estimate      %6.3f%%   %6.3f%%\n", alpha, 1-alpha)
        @printf("Intercept:  %8.7g  [%8.7g, %8.7g]\n", intercept, LCLintercept, UCLintercept)
        @printf("Slope:      %8.7g  [%8.7g, %8.7g]\n", slope, LCLslope, UCLslope)
        results = Dict(:intercept => intercept, :slope => slope,
             :LCLintercept => LCLintercept, :UCLintercept => UCLintercept,
             :LCLslope => LCLslope, :UCLslope => UCLslope)
        !hist || simurationresults(interceptvector, slopevector)
    end
    plotresistantline(x, y, intercept, slope, col1, col2, xlab, ylab, legend)
    return results
end

function resistantline0(x, y, noiteration)
    n = length(x)
    o = sortperm(x)
    x = x[o]
    y = y[o]
    m = n ÷ 3
    n1 = 1:m
    n2 = m+1:n-m
    n3 = n-m+1:n
    mx1 = median(x[n1])
    mx2 = median(x[n2])
    mx3 = median(x[n3])
    a = b = 0
    while true
        my1 = median(y[n1])
        my2 = median(y[n2])
        my3 = median(y[n3])
        da = (my3-my1)/(mx3-mx1)
        db = (my1+my2+my3-da*(mx1+mx2+mx3))/3
        a += da
        b += db
        (noiteration || abs(da/a) < 1e-7) && return b, a
        y .-= da .* x .+ db
    end
end

function simurationresults(interceptvector, slopevector)
    plthist1 = histogram(interceptvector, grid=false, tick_direction=:out,
                            alpha=0.2, label="intercept")
    plthist2 = histogram(slopevector, grid=false, tick_direction=:out,
                            alpha=0.2, label="slope")
    plthist0 = plot(plthist1, plthist2)
    display(plthist0)
end

function plotresistantline(x, y, intercept, slope, col1, col2, xlab, ylab, legend)
    abline(intercept, slope, col, label) =
        plot!([minx, maxx],
              [slope * minx + intercept, slope * maxx + intercept],
              color=col, label=label, legend=legend)
    function lm(x, y)
        b = cov(x, y)/var(x)
        return mean(y) - b * mean(x), b
    end
    minx, maxx = extrema(x)
    margin = (maxx - minx) * 0.05
    minx -= margin
    maxx += margin
    pyplot()
    plt = scatter(x, y, grid=false, tick_direction=:out,
                    xlabel=xlab, ylabel=ylab, color=col1, label="")
    abline(intercept, slope, col1, "resistant line")
    a, b = lm(x, y)
    abline(a, b, col2, "linear regression")
    display(plt)
end

x = 1961:1993
y = [139.130, 140.293, 143.137, 147.350, 150.686, 144.317, 151.207,
    152.882, 156.867, 155.749, 157.735, 162.962, 159.036, 158.589, 149.213,
    148.725, 161.331, 161.363, 158.899, 142.862, 139.085, 162.026, 162.117,
    163.621, 152.982, 170.722, 162.175, 144.809, 167.581, 185.984, 176.457,
    134.477, 134.477];

resistantline(x, y)
# Intercept: -1189.844
#     Slope: 0.681875
# Dict{Symbol,Float64} with 2 entries:
#   :intercept => -1189.84
#   :slope     => 0.681875

resistantline(x, y, n=1000, hist=true)
# Estimate       0.025%    0.975%
# Intercept:  -1189.844  [-1939.097, 1179.805]
# Slope:      0.681875  [ -0.5216, 1.062364]
# Dict{Symbol,Float64} with 6 entries:
# :UCLintercept => 1179.8
# :UCLslope     => 1.06236
# :intercept    => -1189.84
# :slope        => 0.681875
# :LCLslope     => -0.5216
# :LCLintercept => -1939.1

Julia に翻訳--047 Deming 法による回帰直線

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

Deming 法による回帰直線
http://aoki2.si.gunma-u.ac.jp/R/Deming.html

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

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

sxx, sxy, sxy, syy = cov(hcat(x, y)) * (length(x)-1)
が，スマートなのか疑問？

==========#

using Statistics, Printf, Plots

function deming(x, y; n=1, a=1, sig=0.95, col1=:black, col2=:red,
                xlab="", ylab="", legend=:best, hist=false)
    intercept, slope = deming0(x, y, a)
    if n == 1
        @printf("Intercept: %.7g\n", intercept)
        @printf("    Slope: %.7g\n", slope)
        results = Dict(:intercept => intercept, :slope => slope)
    else
        interceptvector = zeros(n)
        slopevector = zeros(n)
        suffix0 = collect(1:length(x))
        for i = 1:n
            suffix = rand(suffix0, length(x))
            interceptvector[i], slopevector[i] = deming0(x[suffix], y[suffix], a)
        end
        alpha = (1-sig)/2
        LCLintercept, UCLintercept = quantile!(interceptvector, [alpha, 1-alpha])
        LCLslope, UCLslope = quantile!(slopevector, [alpha, 1-alpha])
        @printf("            Estimate      %6.3f%%   %6.3f%%\n", alpha, 1-alpha)
        @printf("Intercept:  %8.7g  [%8.7g, %8.7g]\n", intercept, LCLintercept, UCLintercept)
        @printf("Slope:      %8.7g  [%8.7g, %8.7g]\n", slope, LCLslope, UCLslope)
        results = Dict(:intercept => intercept, :slope => slope,
             :LCLintercept => LCLintercept, :UCLintercept => UCLintercept,
             :LCLslope => LCLslope, :UCLslope => UCLslope)
        !hist || simurationresults(interceptvector, slopevector)
    end
    plotdeming(x, y, intercept, slope, col1, col2, xlab, ylab, legend)
    return results
end

function deming0(x, y, a)
    sxx, sxy, sxy, syy = cov(hcat(x, y)) * (length(x)-1)
    if sxy != 0
        Slope = (syy-a*sxx+sqrt((syy-a*sxx)^2+4*a*sxy^2))/(2*sxy)
        Intercept = mean(y)-Slope*mean(x)
    else
        Slope = Intercept = NaN
    end
    return Intercept, Slope
end

function simurationresults(interceptvector, slopevector)
    plthist1 = histogram(interceptvector, grid=false, tick_direction=:out,
                            alpha=0.2, label="intercept")
    plthist2 = histogram(slopevector, grid=false, tick_direction=:out,
                            alpha=0.2, label="slope")
    plthist0 = plot(plthist1, plthist2)
    display(plthist0)
end

function plotdeming(x, y, intercept, slope, col1, col2, xlab, ylab, legend)
    abline(intercept, slope, col, label) =
        plot!([minx, maxx],
              [slope * minx + intercept, slope * maxx + intercept],
              color=col, label=label, legend=legend)
    function lm(x, y)
        b = cov(x, y)/var(x)
        return mean(y) - b * mean(x), b
    end
    minx, maxx = extrema(x)
    margin = (maxx - minx) * 0.05
    minx -= margin
    maxx += margin
    pyplot()
    plt = scatter(x, y, grid=false, tick_direction=:out,
                    xlabel=xlab, ylabel=ylab, color=col1, label="")
    abline(intercept, slope, col1, "Deming")
    a, b = lm(x, y)
    abline(a, b, col2, "linear regression")
    display(plt)
end

x = [594, 422, 516, 432, 300, 452, 498, 281, 445, 717, 507, 582,
    656, 512, 330, 477, 449, 634, 430, 637, 577, 391, 550, 347, 544,
    426, 518, 493, 508, 442, 418, 519, 443, 479, 521, 450, 624, 616,
    338, 544, 383, 526, 706, 488, 589, 425, 655, 527, 519, 562];
y = [616, 401, 424, 396, 416, 335, 447, 322, 507, 680, 588, 614,
    667, 445, 403, 431, 634, 440, 395, 676, 616, 517, 521, 457, 561,
    371, 556, 437, 491, 379, 434, 473, 443, 524, 479, 406, 643, 533,
    480, 519, 422, 423, 651, 406, 531, 617, 720, 468, 478, 608];

deming(x, y)
# Intercept: 0.8898059
#     Slope: 0.9983003
# Dict{Symbol,Float64} with 2 entries:
#   :intercept => 0.889806
#   :slope     => 0.9983

deming(x, y, n=1000, hist=true)
# Estimate       0.025%    0.975%
# Intercept:  0.8898059  [-154.9576, 113.4234]
# Slope:      0.9983003  [0.7733284, 1.287288]
# Dict{Symbol,Float64} with 6 entries:
# :UCLintercept => 113.423
# :UCLslope     => 1.28729
# :intercept    => 0.889806
# :slope        => 0.9983
# :LCLslope     => 0.773328
# :LCLintercept => -154.958

Julia に翻訳--046 Passing & Bablok 法による回帰直線

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

Passing & Bablok 法による回帰直線
http://aoki2.si.gunma-u.ac.jp/R/PassingBablok.html

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

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

こりた。関数内の関数定義はなるべくやらないようにしよう。

==========#

using Statistics, StatsBase, Combinatorics, Random, Printf, Plots

function passingbablok(x, y; n=1, sig=0.95, col1=:black, col2=:red,
                       xlab="", ylab="", legend=:best, hist=false)
    corsign = sign(corkendall(x, y))
    intercept, slope = passingbablok0(x, y, corsign)
    if n == 1
        @printf("Intercept: %.7g\n", intercept)
        @printf("    Slope: %.7g\n", slope)
        results = Dict(:intercept => intercept, :slope => slope)
    else
        interceptvector = zeros(n);
        slopevector = zeros(n);
        suffix0 = collect(1:length(x))
        for i = 1:n
            suffix = rand(suffix0, length(x))
            interceptvector[i], slopevector[i] = passingbablok0(x[suffix], y[suffix], corsign)
        end
        alpha = (1-sig)/2
        LCLintercept, UCLintercept = quantile!(interceptvector, [alpha, 1-alpha])
        LCLslope, UCLslope = quantile!(slopevector, [alpha, 1-alpha])
        @printf("            Estimate      %6.3f%%   %6.3f%%\n", alpha, 1-alpha)
        @printf("Intercept:  %8.7g  [%8.7g, %8.7g]\n", intercept, LCLintercept, UCLintercept)
        @printf("Slope:      %8.7g  [%8.7g, %8.7g]\n", slope, LCLslope, UCLslope)
        results = Dict(:intercept => intercept, :slope => slope,
             :LCLintercept => LCLintercept, :UCLintercept => UCLintercept,
             :LCLslope => LCLslope, :UCLslope => UCLslope)
        !hist || simurationresults(interceptvector, slopevector)
    end
    plotpassingbablok(x, y, intercept, slope, col1, col2, xlab, ylab, legend)
    return results
end

function passingbablok0(x, y, corsign)
    corsign2 = sign(corkendall(x, y))
    suffix = combinations(1:length(x), 2);
    slope = [(y[i] - y[j]) / (x[i] - x[j]) for (i, j) in suffix if x[i] != x[j]];
    length(slope)
    sum(slope .== -1)
    if corsign2 == corsign
        slope = slope[slope .!= -1]# length(slope)
        k = sum(slope .< -1)
        sort!(slope)
    else
        slope = slope[slope .!= 1]
        k = sum(slope .> 1)
        sort!(slope, rev=true)
    end
    n = length(slope)
    Slope = n % 2 == 0 ? (slope[(n+1)÷2+k]+slope[(n+1)÷2+k+1])/2 :
                          slope[(n+1)÷2+k]
    return median(y .- Slope .* x), Slope
end

function simurationresults(interceptvector, slopevector)
    plthist1 = histogram(interceptvector, grid=false, tick_direction=:out,
                            alpha=0.2, label="intercept")
    plthist2 = histogram(slopevector, grid=false, tick_direction=:out,
                            alpha=0.2, label="slope")
    plthist0 = plot(plthist1, plthist2)
    display(plthist0)
end

function plotpassingbablok(x, y, intercept, slope, col1, col2, xlab, ylab, legend)
    abline(intercept, slope, col, label) =
        plot!([minx, maxx],
              [slope * minx + intercept, slope * maxx + intercept],
              color=col, label=label, legend=legend)
    function lm(x, y)
        b = cov(x, y)/var(x)
        return mean(y) - b * mean(x), b
    end
    minx, maxx = extrema(x)
    margin = (maxx - minx) * 0.05
    minx -= margin
    maxx += margin
    pyplot()
    plt = scatter(x, y, xlabel=xlab, ylabel=ylab, color=col1, label="")
    abline(intercept, slope, col1, "Passing-Bablok")
    a, b = lm(x, y)
    abline(a, b, col2, "linear regression")
    display(plt)
end

x = [594, 422, 516, 432, 300, 452, 498, 281, 445, 717, 507, 582,
    656, 512, 330, 477, 449, 634, 430, 637, 577, 391, 550, 347, 544,
    426, 518, 493, 508, 442, 418, 519, 443, 479, 521, 450, 624, 616,
    338, 544, 383, 526, 706, 488, 589, 425, 655, 527, 519, 562];
y = [616, 401, 424, 396, 416, 335, 447, 322, 507, 680, 588, 614,
    667, 445, 403, 431, 634, 440, 395, 676, 616, 517, 521, 457, 561,
    371, 556, 437, 491, 379, 434, 473, 443, 524, 479, 406, 643, 533,
    480, 519, 422, 423, 651, 406, 531, 617, 720, 468, 478, 608];

passingbablok(x, y)
# Intercept: -31.19101
#     Slope: 1.026217
# Dict{Symbol,Float64} with 2 entries:
#   :intercept => -31.191
#   :slope     => 1.02622

passingbablok(x, y, n=1000, hist=true)
# Estimate       0.025%    0.975%
# Intercept:  -31.19101  [-207.6835, 92.82652]
# Slope:      1.026217  [0.8147902, 1.373737]
# Dict{Symbol,Float64} with 6 entries:
# :UCLintercept => 92.8265
# :UCLslope     => 1.37374
# :intercept    => -31.191
# :slope        => 1.02622
# :LCLslope     => 0.81479
# :LCLintercept => -207.683