Julia に翻訳--184 冪曲線回帰 y = a * x^b

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

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

冪曲線回帰 y = a * x**b
http://aoki2.si.gunma-u.ac.jp/Python/power.pdf

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

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

==========#

using Statistics, Printf, Plots

function power(x, y; printout=true, graph=true)
    f(x, a, b) = a * x ^ b
    model = "y = a * x^b"
    n = length(x)
    n >= 2 || error("sample size must be greater than 1")
    (all(x .> 0) && all(y .> 0)) || error("all x[i] && y[i] must be positive")
    a, b = sreg(log.(x), log.(y))
    a = exp(a)
    predicted = f.(x, a, b)
    resid = y .- predicted
    println(model)
    println("a = $a, b = $b")
    if printout
        @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
    if graph
        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)
        plot!(x2, y2, linewidth=0.5, color=:red, label="")
        display(plt)
    end
    Dict(:a => a, :b => b, :predicted => predicted, :resid => resid,
        :x => x, :y => y, :model => model)
end

function sreg(x, y)
    a = cov(x, y) / var(x)
    mean(y) - a * mean(x), a
end

x = collect(1:10);
y = [2, 16, 54, 128, 250, 432, 686, 1024, 1458, 2000];
power(x, y)
# y = a * x^b
# a = 1.9999999999999984, b = 2.9999999999999996
#            x            y        pred.       resid.
#     1.000000     2.000000     2.000000     0.000000
#     2.000000    16.000000    16.000000     0.000000
#     3.000000    54.000000    54.000000     0.000000
#     4.000000   128.000000   128.000000     0.000000
#     5.000000   250.000000   250.000000     0.000000
#     6.000000   432.000000   432.000000     0.000000
#     7.000000   686.000000   686.000000     0.000000
#     8.000000  1024.000000  1024.000000     0.000000
#     9.000000  1458.000000  1458.000000     0.000000
#    10.000000  2000.000000  2000.000000     0.000000

Julia に翻訳--183 特殊な指数曲線回帰 b^y = a * x

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

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

特殊な指数曲線回帰 b**y = a * x
http://aoki2.si.gunma-u.ac.jp/Python/exp3.pdf

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

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

==========#

using Statistics, Printf, Plots

function exp3(x, y; printout=true, graph=true)
    f(x, a, b) = a + b * x
    model = "b^y = a * x"
    n = length(x)
    n >= 2 || error("sample size must be greater than 1")
    all(x .> 0) || error("all x[i] must be positive")
    logx = log.(x)
    intercept, slope = sreg(logx, y)
    b = exp(1 / slope)
    a = exp(intercept / slope)
    #predicted = intercept + slope*logx
    predicted = f.(logx, intercept, slope)
    resid = y - predicted
    println(model)
    println("a = $a, b = $b")
    if printout
        @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
    if graph
        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.(log.(x2), intercept, slope)
        plot!(x2, y2, linewidth=0.5, color=:red, label="")
        display(plt)
    end
    Dict(:a => a, :b => b, :intercept => intercept, :slope => slope, :predicted => predicted, :resid => resid, :x => x, :logx => logx, :y => y, :model => model)
end

function sreg(x, y)
    a = cov(x, y) / var(x)
    mean(y) - a * mean(x), a
end

x = collect(1:10);
y = [0.631, 1.262, 1.631, 1.893, 2.096, 2.262, 2.402, 2.524, 2.631, 2.727];
exp3(x, y)
# b^y = a * x
# a = 2.000197950638151, b = 2.999971611554844
#            x            y        pred.       resid.
#     1.000000     0.631000     0.631025    -0.000025
#     2.000000     1.262000     1.261960     0.000040
#     3.000000     1.631000     1.631034    -0.000034
#     4.000000     1.893000     1.892896     0.000104
#     5.000000     2.096000     2.096011    -0.000011
#     6.000000     2.262000     2.261969     0.000031
#     7.000000     2.402000     2.402284    -0.000284
#     8.000000     2.524000     2.523831     0.000169
#     9.000000     2.631000     2.631043    -0.000043
#    10.000000     2.727000     2.726947     0.000053

Julia に翻訳--182 指数曲線回帰 y = a * b^x

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

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

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

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

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

==========#

using Statistics, Printf, Plots

function exp1(x, y; printout=true, graph=true)
    f(x, a, b) = a * b ^ x
    model = "y = a * b^x"
    n = length(x)
    n >= 2 || error("sample size must be greater than 1")
    all(y .> 0) || error("all y[i] must be positive")
    a, b = exp.(sreg(x, log.(y)))
    predicted = f.(x, a, b)
    resid = y .- predicted
    println(model)
    println("a = $a, b = $b")
    if printout
        @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
    if graph
        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)
        plot!(x2, y2, linewidth=0.5, color=:red, label="")
        display(plt)
    end
    Dict(:a => a, :b => b, :predicted => predicted, :resid => resid,
        :x => x, :y => y, :model => model)
end

function sreg(x, y)
    a = cov(x, y) / var(x)
    mean(y) - a * mean(x), a
end

x = collect(1:10);
y = [3, 4.5, 6.75, 10.125, 15.188, 22.781, 34.172, 51.258, 76.887, 115.330];
exp1(x, y)

# y = a * b^x
# a = 2.0000061235873563, b = 1.5000000443425876
#            x            y        pred.       resid.
#     1.000000     3.000000     3.000009    -0.000009
#     2.000000     4.500000     4.500014    -0.000014
#     3.000000     6.750000     6.750021    -0.000021
#     4.000000    10.125000    10.125032    -0.000032
#     5.000000    15.188000    15.187549     0.000451
#     6.000000    22.781000    22.781324    -0.000324
#     7.000000    34.172000    34.171987     0.000013
#     8.000000    51.258000    51.257982     0.000018
#     9.000000    76.887000    76.886975     0.000025
#    10.000000   115.330000   115.330465    -0.000465

Julia に翻訳--181 ハーディー・ワインベルグ平衡

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

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

ハーディー・ワインベルグ平衡
http://aoki2.si.gunma-u.ac.jp/Python/Hardy_Weinberg_test.pdf

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

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

当然かも知れないが，logfactorial() の引数は整数でなくてはならない。

==========#

using SpecialFunctions, Rmath

function hardyweinbergtest(n11, n12, n22)
    exactp(n11, n12, n22) = exp(log(2) * n12 - lgamma(n11+1) - lgamma(n12+1) - lgamma(n22+1) + base)
    MACHINEDOUBLEEPS = eps(Float64)
    chisq, asymptoticpvalue = chi2contingency(n11, n12 / 2, n12 / 2, n22)
    n1 = 2 * n11 + n12
    n2 = 2 * n22 + n12
    n = n11 + n12 + n22
    base = lgamma(n+1) + lgamma(n1+1) + lgamma(n2+1) - lgamma(2 * n+1)
    p0 = exactp(n11, n12, n22)
    p = [exactp((n1 - x12) / 2, x12, (n2 - x12) / 2)
         for x12 = n1 % 2:2: min(n1, n2)]
    exactpvalue = sum([p2 for p2 = p if p2 <= p0 + MACHINEDOUBLEEPS])
    println("Exact test for Hardy-Weinberg equilibrium")
    println("chisq. = $chisq,  df = 1,  asymptotic p value = $asymptoticpvalue")
    println("exact p value = $exactpvalue")
    Dict(:N11 => n11, :N12 => n12, :N22 => n22, :N1 => n1, :N2 => n2,
        :exactpvalue => exactpvalue, :chisq => chisq, :df => 1,
        :asymptoticpvalue => asymptoticpvalue)
end

function chi2contingency(a, b, c, d)
    n = a + b + c + d
    dif = a * d - b * c
    dif = dif > n / 2 ? dif - n / 2 : 0
    chisq = n * dif^2 / (a + b) / (c + d) / (a + c) / (b + d)
    chisq, pchisq(chisq, 1, false)
end

hardyweinbergtest(6, 1, 2)
# Exact test for Hardy-Weinberg equilibrium
# chisq. = 1.7914792899408287,  df = 1,  asymptotic p value = 0.1807460300456809
# exact p value = 0.05882352941176444

hardyweinbergtest(45, 25, 12)
# Exact test for Hardy-Weinberg equilibrium
# chisq. = 4.854012844364627,  df = 1,  asymptotic p value = 0.02758189277568519
# exact p value = 0.01685913919443042

Julia に翻訳--180 高次式方程式の解

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

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

高次式方程式の解
http://aoki2.si.gunma-u.ac.jp/Python/Bairstow.pdf

ファイル名: bairstow.jl　　関数名: bairstow, quadraticroot

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

julia の Complex 型はちょっと扱いづらい。

sqrt(−2) は自動的に計算されないので，sqrt(-2+0im) とか sqrt(Complex(-2)), sqrt(complex(-2, 0)) と書かなければならない。
Vector{Float64} に append！() しようとすると，inexact error になるので，Vector{Complex} に append！() しなくてはならない。
sort() は　by=abs を指定しないといけない。

==========#

function quadraticroot(p, q)
    d = p^2 - 4 * q
    if d >= 0
        x1 = (-p - sign(p) * sqrt(d)) / 2
        x2 = q / x1
    else
        x1 = (-p + sqrt(complex(d, 0))) / 2
        x2 = conj(x1)
    end
    x1, x2
end

function bairstow(a; EPSILON = 1e-14, MAXLOOP = 10000)
    b = copy(a)
    while b[end] == 0
        pop!(b)
    end
    ans0 = zeros(Complex, length(a) - length(b))
    a = vec(b)
    a /= a[1]
    while true
        n = length(a) - 1
        if n == 2
            append!(ans0, quadraticroot(a[2], a[3]))
            break
        elseif n == 1
            append!(ans0, -a[2] / a[1])
            break
        end
        p = q = 1
        b = ones(length(a))
        c = ones(length(a))
        ill = true
        d = 999
        for i = 1:MAXLOOP
            b[2] = a[2] - p * b[1]
            c[2] = b[2] - p * c[1]
            for j = 3:n+1
                b[j] = a[j] - p * b[j - 1] - q * b[j - 2]
                c[j] = b[j] - p * c[j - 1] - q * c[j - 2]
            end
            d = c[n - 1]^2 - c[n - 2] * (c[n] - b[n])
            deltap = (b[n] * c[n - 1] - b[n+1] * c[n - 2]) / d
            deltaq = (b[n+1] * c[n - 1] - b[n] * (c[n] - b[n])) / d
            if abs(deltap) < EPSILON && abs(deltaq) < EPSILON
                ill = false
                break
            end
            p += deltap
            q += deltaq
        end
        !ill || error("ill condition. p = $p, q = $q, d = $d")
        append!(ans0, quadraticroot(p, q))
        a = copy(b[1:n-1])
    end
    sort!(ans0, by=abs)
    for x in ans0
        println(imag(x) == 0 ? real(x) : x)
    end
    ans
end

a = [1, -15, 85, -225, 274, -120];
bairstow(a)
# 0.9999999999999998
# 2.000000000000007
# 2.999999999999973
# 4.0000000000000435
# 4.999999999999977

# 5-element Vector{Complex}:
#  0.9999999999999998 + 0.0im
#   2.000000000000007 + 0.0im
#   2.999999999999973 + 0.0im
#  4.0000000000000435 + 0.0im
#   4.999999999999977 + 0.0im

x = bairstow([1,2,3,4,5,6])
# 0.5516854634589815 + 1.2533488602772063im
# 0.5516854634589815 - 1.2533488602772063im
# -0.8057864693890313 + 1.2229047133744098im
# -0.8057864693890313 - 1.2229047133744098im
# -1.4917979881399004

# 5-element Vector{Complex}:
#   0.5516854634589815 + 1.2533488602772063im
#   0.5516854634589815 - 1.2533488602772063im
#  -0.8057864693890313 + 1.2229047133744098im
#  -0.8057864693890313 - 1.2229047133744098im
#  -1.4917979881399004 + 0.0im

Julia に翻訳--178 生存期間の差の検定に必要なサンプルサイズ

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

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

生存期間の差の検定に必要なサンプルサイズ
http://aoki2.si.gunma-u.ac.jp/Python/surv_sample.pdf

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

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

==========#

using Rmath

function survsample(;duration, more, altime, d=1.5, α=0.05, power=0.8)
    φ2(al, du) = du / al^3 / (du / al - 1 + exp(-du / al))
    φ3(al, du, mo) = 1 / al^2 / (1 - (exp(-(mo - du) / al) - exp(-mo / al)) * al / du)
    method = "生存期間の差の検定に必要なサンプルサイズ"
    if α <= 0 || α >= 1 || power <= 0 || power >= 1 || d == 0 || duration < 0 || altime < 0 || more < 0
        error("input value error")
    end
    Zα = qnorm.([α, α/2], false)
    Zβ = qnorm(power)
    if duration == 0 && more == 0
        n1, n2 = 2(Zα .+ Zβ).^2 / log(d)^2
        message = "一括登録の"
    elseif duration != 0 && more == 0
        altime > 0 || error("対照群の生存期間を入力のこと")
        n1, n2 = (Zα .+ Zβ).^2 * (φ2(altime, duration) + φ2(altime * d, duration)) / (1 / altime / d - 1 / altime)^2
        message = "継続登録の"
    elseif duration != 0 && more != 0
        altime > 0 || error("対照群の生存期間を入力のこと")
        more += duration
        n1, n2 = (Zα .+ Zβ).^2 * (φ3(altime, duration, more) + φ3(altime * d, duration, more)) / (1 / altime / d - 1 / altime)^2
        message = "登録終了後に観察期間を設ける"
    else
        error("入力値が変です")
    end
    message *= "場合について計算した，各群あたりのサンプルサイズです"
    println(method)
    println("片側検定の場合: n = $n1,  両側検定の場合: n = $n2,  有意水準 = $α,  検出力 = $power")
    println(message)
    Dict(:n1 => n1, :n2 => n2, :α => α, :power => power, :method => method, :message => message)
end

survsample(duration=0, more=0, altime=0)
# 生存期間の差の検定に必要なサンプルサイズ
# 片側検定の場合: n = 75.21269772788015,  両側検定の場合: n = 95.4840200919946,  有意水準 = 0.05,  検出力 = 0.8
# 一括登録の場合について計算した，各群あたりのサンプルサイズです

survsample(duration=5, more=0, altime=3)
# 生存期間の差の検定に必要なサンプルサイズ
# 片側検定の場合: n = 170.80410875138924,  両側検定の場合: n = 216.83922322290763,  有意水準 = 0.05,  検出力 = 0.8
# 継続登録の場合について計算した，各群あたりのサンプルサイズです

survsample(duration=4, more=1, altime=3)
# 生存期間の差の検定に必要なサンプルサイズ
# 片側検定の場合: n = 144.75696749266226,  両側検定の場合: n = 183.77185781227448,  有意水準 = 0.05,  検出力 = 0.8
# 登録終了後に観察期間を設ける場合について計算した，各群あたりのサンプルサイズです