Julia に翻訳--188 ポリシリアル相関係数

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

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

ポリシリアル相関係数
http://aoki2.si.gunma-u.ac.jp/Python/polyserial.pdf

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

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

Julia には（まだ？） pmvnorm（） がないので，R の pmvnorm() を呼ぶ。

optimize() も，hessian を明示的には得られないので，トリッキーなことをする（方法が，インターネットから得られた）。

==========#

using Rmath, Statistics, Optim, Printf # LinearAlgebra

function polyserial(x, y; twostep=false, bins=4)
    function f(pars)
        rho = pars[1]
        if length(pars) == 1
            cts = vcat(-Inf, rowCuts, Inf)
        else
            temp = pars[2:end]
            if any(diff(temp) .< 0)
                return Inf
            end
            cts = vcat(-Inf, temp, Inf)
        end
        tau = -(rho .* z .- cts') ./ sqrt(1 - rho^2)
        choose = pnorm.(tau[i, y[i]+1] for i = 1:n)
        choose2 = pnorm.(tau[i, y[i]] for i = 1:n)
        -sum(log.(dnorm.(z) .* (choose .- choose2)))
    end

    method = "Polyserial Correlation:  " * (twostep ? "two step " : "") * "ML estimator"
    z = zscore(x)
    val, tab = table(y)
    n = sum(tab)
    s = length(tab)
    rowCuts = qnorm.(cumsum(tab) ./ n)[1:s-1]
    initval = std(y, corrected=false) * cor(x, y) / sum(dnorm.(rowCuts))
    if twostep
        result = optimize(f, [initval], BFGS())
        rho = Optim.minimizer(result)
        SE = sqrt(invH)
    else # pars = vcat(initval, rowCuts)
        result = optimize(f, vcat(initval, rowCuts), BFGS())
        x = Optim.minimizer(result)
        rho, rowCuts = x[1], x[2:end]
        hes = sqrt.(invH[i, i] for i = 1:s)
        SE, ser = hes[1], hes[2:end]
    end
    chisq = chisqfunc(y, z, rho, rowCuts, bins)
    df = s * bins - s - 1
    pvalue = pchisq(chisq, df, false)
    println("$method\n  ML estimator = $rho,  Std.Err. = $SE")
    println("Test of bivariate normality:\n  chisq = $chisq,  df = $df,  p value = $pvalue")
    if !twostep
        @printf("%13s  %9s\n", "Treshold", "Std. Err.")
        for i = 1:length(rowCuts)
            @printf("%2d  %9.6f  %9.6f\n", i, rowCuts[i], ser[i])
        end
    end
    Dict(:method => method, :rho => rho, :SE => SE, :chisq => chisq, :df => df, :pvalue => pvalue, :rowCuts => rowCuts)
end

Optim.after_while!(d, state::Optim.BFGSState, method::BFGS, options) = global invH = state.invH

function binBvn(rho, rowCuts, colCuts)
    rowCuts = vcat(-Inf, rowCuts, Inf)
    colCuts = vcat(-Inf, colCuts, Inf)
    r = length(rowCuts) - 1
    c = length(colCuts) - 1
    P = zeros(r, c)
    S = [1 rho; rho 1]
    mu = zeros(2)
    for i = 1:r
        for j = 1:c
            P[i, j] = pmvnorm([rowCuts[i], colCuts[j]],
                              [rowCuts[i + 1], colCuts[j + 1]], mu, S)
        end
    end
    P
end

using RCall
pmvnorm(lower, upper, mu, S) = rcopy(R"library(mvtnorm); pmvnorm($lower, $upper, $mu, $S)")

zscore(x) = (x .- mean(x)) / std(x, corrected=false)

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 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

function chisqfunc(y, z, rho, rowCuts, bins)
    colCuts = qnorm.((1:bins-1) ./ bins)
    P = binBvn(rho, rowCuts, colCuts)
    catz = [sum(w .> vcat(-Inf, colCuts, Inf)) for w in zscore(z)]
    indexx, indexy, tab = table(y, catz)
    nonzero = tab .!= 0
    n = length(y)
    2sum(tab[nonzero] .* log.(tab[nonzero] ./ (P[nonzero] .* n)))
end

x = [0.7492, -0.22295, 0.11472, 0.53715, -0.51242, 0.35807, 0.49264,
    -0.51353, -0.94197, 1.15984, 1.12984, -1.02436, -1.07608, -0.30467,
    0.54926, 1.35279, 2.40187, 0.37274, -0.74321, 1.71418, 0.47398,
    -0.78159, 1.2034, -1.21809, 0.22509, -0.01787, 0.14278, -0.57617,
    0.88083, 1.46454, -0.20298, 0.94596, -1.04114, 1.2703, 0.1639,
    -0.15033, -0.41042, 1.59887, 0.58806, 0.72434, 0.89338, 0.34713,
    1.42137, 0.56944, 0.73173, -1.15237, 1.72124, -0.76932, 0.07062,
    -0.75912, -0.08296, 0.06515, -0.00246, 0.10176, -0.34965, -1.53095,
    1.55971, -0.20873, 0.11238, 2.08786, 1.83424, 0.70741, -0.70681,
    -1.48295, -2.36075, 1.81494, 0.20038, 2.29407, -0.80354, 0.0154,
    -1.43625, -0.96287, 0.15512, -1.27466, 0.90138, -1.42222, -0.02011,
    -0.16987, 1.55233, -0.10292, -0.39256, -0.35841, -1.03065, 0.29056,
    1.69807, 0.54179, 0.24826, -1.52018, 2.50281, -0.2621, -0.89337,
    -0.09901, 0.80746, 1.50104, -0.29506, 0.39425, 0.83788, 1.47223,
    0.70706, -0.42811]

y = [3, 2, 1, 2, 2, 4, 3, 3, 2, 3, 4, 1, 3, 2, 2, 4, 4, 3, 1, 2,
    1, 3, 4, 2, 1, 4, 2, 2, 3, 1, 1, 4, 1, 4, 1, 2, 4, 3, 2, 3, 3,
    4, 4, 2, 3, 1, 2, 2, 4, 2, 1, 1, 2, 3, 2, 1, 4, 2, 3, 2, 4, 4,
    1, 3, 1, 3, 1, 2, 2, 2, 1, 1, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
    1, 2, 3, 3, 2, 3, 1, 2, 3, 2, 3, 2, 2, 2, 3, 2, 2]

polyserial(x, y, twostep=true)
# Polyserial Correlation:  two step ML estimator
#   ML estimator = [0.46699921329233146],  Std.Err. = [0.07931271523690714]
# Test of bivariate normality:
#   chisq = 5.488940449811883,  df = 11,  p value = 0.9051999325887348

polyserial(x, y)
# Polyserial Correlation:  ML estimator
#   ML estimator = 0.46619488420452654,  Std.Err. = 0.08060954580245325
# Test of bivariate normality:
#   chisq = 5.484065294038735,  df = 11,  p value = 0.9054810695989302
#      Treshold  Std. Err.
#  1  -0.845491   0.133867
#  2   0.305882   0.118075
#  3   1.042709   0.144813

polyserial(x, y, bins=6)
# Polyserial Correlation:  ML estimator
# Polyserial Correlation:
#   ML estimator = 0.46619488420452654,  Std.Err. = 0.08060954580245325
# Test of bivariate normality:
#   chisq = 12.914187782275436,  df = 19,  p value = 0.8429348531100296
#      Treshold  Std. Err.
#  1  -0.845491   0.133867
#  2   0.305882   0.118075
#  3   1.042709   0.144813

Julia に翻訳--187 ゴンペルツ曲線回帰 y = a * b^exp(-c*x)

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

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

ゴンペルツ曲線回帰 y = a * b**exp(-c*x)
http://aoki2.si.gunma-u.ac.jp/Python/gomperts.pdf

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

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

for ループが変数スコープを形成するので，for ループ内から結果を表示する関数を呼ぶようにした。

==========#

using Statistics, Printf, Plots

function gomperts(y, maxit=1000, tol=1e-6, printout=true, graph=true)
    model = "y = a * b^exp(-c*x)"
    n = length(y)
    n > 3 || error("sample size must be greater than 3")
    all(y .> 0) || error("all y[i] must be positive")
    d = zeros(n, 3)
    u1 = u0 = init(y, n)
    d[:, 3] = log.(y)
    x = 1:n
    for i = 1:maxit
        d[:, 1] = exp.(-u1 .* x)
        d[:, 2] = x .* d[:, 1]
        a, b, c = mreg(d)
        u0 = u1
        u1 -= c / b
        if abs((u0 - u1) / u0) < tol
            a = exp(a)
            b = exp(b)
            c = u1
            predicted = f.(x, a, b, c)
            resid = y - predicted
            printout && printout2(model, a, b, c, n, x, y, predicted, resid)
            graph    && graph2(model, a, b, c, x, y)
            return Dict(:a => a, :b => b, :c => c, :predicted => predicted,
                        :resid => resid,
                        :x => x, :y => y, :model => model)
        end
    end
    error("not converged")
end

f(x, a, b, c) = a * b ^ exp(-c * x)

function printout2(model, a, b, c, n, x, y, predicted, resid)
    println(model)
    println("a = $a, b = $b, c = $c")
    @printf("%12s %12s %12s %12s\n", "x", "y", "pred.", "resid.")
    for i =1:n
        @printf("%12.6f %12.6f %12.6f %12.6f\n", x[i], y[i], predicted[i], resid[i])
    end
end

function graph2(model, a, b, c, x, y)
    model = "y = a * b^{\\exp(-c*x)}"
    plt = scatter(x, y, grid=false, tick_direction=:out, title="\$$model\$",
                  color=:blue, markerstrokecolor=:blue, label="")
    min, max = extrema(x)
    w = 0.05(max - min)
    x2 = range(min - w, max + w, length=1000)
    y2 = f.(x2, a, b, c)
    plot!(x2, y2, linewidth=0.5, color=:red, label="")
    display(plt)
end

function init(y, n)
    m = length(y) ÷ 3
    s1 = sum(log.(y[1:m]))
    s2 = sum(log.(y[m+1:2m]))
    s3 = sum(log.(y[2m+1:end]))
    log((s1 - s2) / (s2 - s3)) / m
end

function mreg(dat) # 回帰分析を行い，偏回帰係数を返す関数
    r = cor(dat)
    betas = r[1:end-1, 1:end-1] \ r[1:end-1, end]
    SS = cov(dat, corrected=false) .* size(dat, 1)
    prop = [SS[i, i] for i in 1:size(SS, 1)]
    b = betas .* sqrt.(prop[end] ./ prop[1:end-1])
    means = mean(dat, dims=1)
    b0 = means[end] - sum(b .* means[1:end-1])
    vcat(b0, b)
end

y = [2.9, 5.2, 9.1, 15.5, 25.0, 37.8, 52.6, 66.9, 78.6, 87.0, 92.4, 95.7, 97.6, 98.6, 99.2];
gomperts(y)
# y = a * b^exp(-c*x)
# a = 134.13510703273636, b = 0.006807921520544414, c = 0.22478593362143764
#            x            y        pred.       resid.
#     1.000000     2.900000     2.493440     0.406560
#     2.000000     5.200000     5.561857    -0.361857
#     3.000000     9.100000    10.555993    -1.455993
#     4.000000    15.500000    17.609913    -2.109913
#     5.000000    25.000000    26.501587    -1.501587
#     6.000000    37.800000    36.732508     1.067492
#     7.000000    52.600000    47.674639     4.925361

Julia に翻訳--186 ロジスティック曲線回帰 y = a / (1 + b*exp(-c*x))

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

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

ロジスティック曲線回帰 y = a / (1 + b*exp(-c*x))
http://aoki2.si.gunma-u.ac.jp/Python/logistic.pdf

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

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

for ループが変数スコープを形成するので，for ループ内から結果を表示する関数を呼ぶようにした。

==========#

using Statistics, Printf, Plots

function logistic(y; maxit=1000, tol=1e-6, printout=true, graph=true)
    model = "y = a / (1 + b * exp(-c * x))"
    n = length(y)
    n > 3 || error("sample size must be greater than 3")
    all(y .> 0) || error("all y[i] must be positive")
    d = zeros(n, 3)
    u1 = u0 = init(y)
    d[:, 3] = 1 ./ y
    x = 1:n
    for i = 1:maxit
        d[:, 1] = exp.(-u1 .* x)
        d[:, 2] = x .* d[:, 1]
        a, b, c = mreg(d)
        u0 = u1
        u1 -= c / b
        if abs((u0 - u1) / u0) < tol
            a = 1 / a
            b *= a
            c = u1
            predicted = f.(x, a, b, c)
            resid = y - predicted
            printout && printout2(model, a, b, c, n, x, y, predicted, resid)
            graph    && graph2(model, a, b, c, x, y)
            return Dict(:a => a, :b => b, :c => c, :predicted => predicted,
                        :resid => resid,
                        :x => x, :y => y, :model => model)
        end
    end
    error("not converged")
end

f(x, a, b, c) = a / (1 + b * exp(-c * x))

function printout2(model, a, b, c, n, x, y, predicted, resid)
    println(model)
    println("a = $a, b = $b, c = $c")
    @printf("%12s %12s %12s %12s\n", "x", "y", "pred.", "resid.")
    for i =1:n
        @printf("%12.6f %12.6f %12.6f %12.6f\n", x[i], y[i], predicted[i], resid[i])
    end
end

function graph2(model, a, b, c, x, y)
    model = replace(model, "exp" => "\\exp")
    plt = scatter(x, y, grid=false, tick_direction=:out, title="\$$model\$",
                  color=:blue, markerstrokecolor=:blue, label="")
    min, max = extrema(x)
    w = 0.05(max - min)
    x2 = range(min - w, max + w, length=1000)
    y2 = f.(x2, a, b, c)
    plot!(x2, y2, linewidth=0.5, color=:red, label="")
    display(plt)
end

function init(y)
    m = length(y) ÷ 3
    s1 = sum(log.(y[1:m]))
    s2 = sum(log.(y[m+1:2m]))
    s3 = sum(log.(y[2m+1:end]))
    log((s1 - s2) / (s2 - s3)) / m
end

function mreg(dat) # 回帰分析を行い，偏回帰係数を返す関数
    r = cor(dat)
    betas = r[1:end-1, 1:end-1] \ r[1:end-1, end]
    SS = cov(dat, corrected=false) .* size(dat, 1)
    prop = [SS[i, i] for i in 1:size(SS, 1)]
    b = betas .* sqrt.(prop[end] ./ prop[1:end-1])
    means = mean(dat, dims=1)
    b0 = means[end] - sum(b .* means[1:end-1])
    vcat(b0, b)
end

y = [2.9, 5.2, 9.1, 15.5, 25.0, 37.8, 52.6, 66.9, 78.6, 87.0, 92.4,95.7, 97.6, 98.6, 99.2];
logistic(y)
# y = a / (1 + b * exp(-c * x))
# a = 99.10621619280843, b = 60.7563899945766, c = 0.6055235054694381
#            x            y        pred.       resid.
#     1.000000     2.900000     2.901223    -0.001223
#     2.000000     5.200000     5.189233     0.010767
#     3.000000     9.100000     9.110771    -0.010771
#     4.000000    15.500000    15.506534    -0.006534
#     5.000000    25.000000    25.138002    -0.138002
#     6.000000    37.800000    38.030374    -0.230374
#     7.000000    52.600000    52.813748    -0.213748
#     8.000000    66.900000    67.036298    -0.136298
#     9.000000    78.600000    78.586916     0.013084
#    10.000000    87.000000    86.744501     0.255499
#    11.000000    92.400000    91.954133     0.445867
#    12.000000    95.700000    95.070403     0.629597
#    13.000000    97.600000    96.862005     0.737995
#    14.000000    98.600000    97.868622     0.731378
#    15.000000    99.200000    98.426898     0.773102

Julia に翻訳--185 漸近指数曲線回帰 y = a * b^x + c

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

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

漸近指数曲線回帰 y = a * b**x + c
http://aoki2.si.gunma-u.ac.jp/Python/asyr.pdf

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

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

for ループが変数スコープを形成するので，for ループ内から結果を表示する関数を呼ぶようにした。

==========#

using Statistics, Printf, Plots

function asyr(x, y; maxiter=1000, tol=1e-6, printout=true, graph=true)
    model = "y = a * b^x + c"
    n = length(x)
    n > 3 || error("sample size must be greater than 3")
    d = zeros(n, 3)
    d[:, 3] = y
    u1 = init(y, n)
    for loop = 1:maxiter
        d[:, 1] = u1 .^ x
        d[:, 2] = x .* u1 .^ (x .- 1)
        c, a, b = mreg(d)
        u0 = u1
        u1 += b / a
        if abs((u0 - u1) / u0) < tol
            b = u1
            predicted = f.(x, a, b, c)
            resid = y - predicted
            printout && printout2(model, a, b, c, n, x, y, predicted, resid)
            graph    && graph2(model, a, b, c, x, y)
            return Dict(:a => a, :b => b, :c => c, :predicted => predicted,
                        :resid => resid,
                        :x => x, :y => y, :model => model)
        end
    end
    error("not converged")
end

f(x, a, b, c) = a * b^x + c

function printout2(model, a, b, c, n, x, y, predicted, resid)
    println(model)
    println("a = $a, b = $b, c = $c")
    @printf("%12s %12s %12s %12s\n", "x", "y", "pred.", "resid.")
    for i =1:n
        @printf("%12.6f %12.6f %12.6f %12.6f\n", x[i], y[i], predicted[i], resid[i])
    end
end

function graph2(model, a, b, c, x, y)
    plt = scatter(x, y, grid=false, tick_direction=:out, title="\$$model\$",
                  color=:blue, markerstrokecolor=:blue, label="")
    min, max = extrema(x)
    w = 0.05(max - min)
    x2 = range(min - w, max + w, length=1000)
    y2 = f.(x2, a, b, c)
    plot!(x2, y2, linewidth=0.5, color=:red, label="")
    display(plt)
end

function init(d, n)
    if n == 4
        (4d[4]+d[3]-5d[2])/(4d[3]+d[2]-5d[1])
    elseif n == 5
        (4d[5]+3d[4]-d[3]-6d[2])/(4d[4]+3d[3]-d[2]-6d[1])
    elseif n == 6
        (4d[6]+4d[5]+2d[4]-3d[3]-7d[2])/(4d[5]+4d[4]+2d[3]-3d[2]-7d[1])
    elseif n == 7
        (d[7]+d[6]+d[5]-d[3]-2d[2])/(d[6]+d[5]+d[4]-d[2]-2d[1])
    else
        k = n/7
        (d[floor(Int, 6*k)+1]+d[floor(Int, 5*k)+1]+d[floor(Int, 4*k)+1]-d[floor(Int, 2*k)+1]-2d[floor(Int, k)+1]) /
            (d[floor(Int, 5*k)+1]+d[floor(Int, 4*k)+1]+d[floor(Int, 3*k)+1]-d[floor(Int, k)+1]-2d[1])
    end
end

function mreg(dat) # 回帰分析を行い，偏回帰係数を返す関数
    r = cor(dat)
    betas = r[1:end-1, 1:end-1] \ r[1:end-1, end]
    SS = cov(dat, corrected=false) .* size(dat, 1)
    prop = [SS[i, i] for i in 1:size(SS, 1)]
    b = betas .* sqrt.(prop[end] ./ prop[1:end-1])
    means = mean(dat, dims=1)
    b0 = means[end] - sum(b .* means[1:end-1])
    vcat(b0, b)
end

x = 0:0.5:5
y = [52, 53, 56, 60, 68, 81, 104, 144, 212, 331, 536];
asyr(x, y)
# y = a * b^x + c
# a = 2.0282343471965527, b = 2.991960122729875, c = 49.73508957252342
#            x            y        pred.       resid.
#     0.000000    52.000000    51.763324     0.236676
#     0.500000    53.000000    53.243384    -0.243384
#     1.000000    56.000000    55.803486     0.196514
#     1.500000    60.000000    60.231767    -0.231767
#     2.000000    68.000000    67.891489     0.108511
#     2.500000    81.000000    81.140729    -0.140729
#     3.000000   104.000000   104.058313    -0.058313
#     3.500000   144.000000   143.699509     0.300491
#     4.000000   212.000000   212.268009    -0.268009
#     4.500000   331.000000   330.872886     0.127114
#     5.000000   536.000000   536.027103    -0.027103