Julia に翻訳--164 ニュートン・ラフソン法， 1 変数方程式の解

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

ニュートン・ラフソン法による 1 変数方程式の解
http://aoki2.si.gunma-u.ac.jp/R/newton.html

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

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

ラムダ関数の指定法（-> を使う）はちょっとなじめないが，
  x  -> x ^ 3 - 2
f(x) =  x ^ 3 - 2
なので，慣れればどうということもない。

==========#

function newton(fun, x1; delta = 1e-5, epsilon = 1e-14, maxrotation = 100)
    fun2(x) = (fun(x + delta) - fun(x)) / delta
    for i = 1:maxrotation
        x2 = x1 - fun(x1) / fun2(x1)
        if abs((x2 - x1) / x2) < epsilon
            return x2
        end
        x1 = x2
    end
    error("収束しませんでした")
end

関数定義する場合

fun(x) = sin(x) / x
newton(fun , 1) # 3.141592653589793

無名関数（ラムダ関数）で定義する場合

newton(x -> x ^ 3 - 2, 1) # 1.2599210498948732

Julia に翻訳--163 日付に関する関数

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

日付に関する関数
http://aoki2.si.gunma-u.ac.jp/R/date.html

ファイル名: date.jl　　関数名: dw, jday, date の 3 関数

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

÷ を使うのが，ちょっと気恥ずかしい

==========#

using Printf

dw(j) = ["Sun", "Mon", "Tue", "Wed", "Thu", "Fri", "Sat"][(j + 1) % 7 + 1]

function jday(iy, jm, kd)
    tmp = -(jm < 3)
    kd - 32075 + (1461 * (iy + 4800 + tmp)) ÷ 4 + (367 * (jm - 2 - tmp * 12)) ÷ 12 - (3 * ((iy + 4900 + tmp) ÷ 100)) ÷ 4
end

function date(jul)
    l = jul + 68569
    n = (4 * l) ÷ 146097
    l = l - (146097 * n + 3) ÷ 4
    iy = (4000 * (l + 1)) ÷ 1461001
    l = l - (1461 * iy) ÷ 4 + 31
    jm = (80 * l) ÷ 2447
    kd = l - (2447 * jm) ÷ 80
    l = jm÷ 11
    jm = jm + 2 - 12 * l
    iy = 100 * (n - 49) + iy + l
    @sprintf("%4d/%02d/%02d", iy, jm, kd)
end

jday(2004, 3, 7) # 2453072
jday(2000, 1, 1) - jday(1900, 1, 1) # 36524
dw(jday(2004, 3, 7)) # "Sun"
date(jday(2004, 2, 28) + 1) # "2004/02/29"

Julia に翻訳--162 類似度指数

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

類似度指数
http://aoki2.si.gunma-u.ac.jp/R/Morisita.html

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

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

==========#

function morisita(a, b)
    length(a) == length(b) || error("length differ")
    na = sum(a)
    nb = sum(b)
    2 * sum(a .* b) / ((sum(a .* (a .- 1) / na / (na - 1)) + sum(b .* (b .- 1) / nb / (nb - 1))) * na * nb)
end

morisita([3, 2, 3, 4, 3, 4], [2, 6, 4, 3, 2, 1]) # 1.0184331797235024

Julia に翻訳--161 多様度指数

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

多様度指数
http://aoki2.si.gunma-u.ac.jp/R/Heterogeneity.html

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

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

base のデフォルト指定が冗長だった

==========#

function heterogeneity(x; base = ℯ)
    x = x[x .> 0] ./ sum(x)
    -sum(x .* log.(x) ./ log(base))
end

heterogeneity([3, 2, 1, 4, 3, 4]) # 1.7115472862863528

heterogeneity([3, 2, 1, 4, 3, 4], base = 2) # 2.469240782172284

heterogeneity([3, 2, 1, 4, 3, 4], base = 10) # 0.7433155419506481

 

atom と Julia の共同作業について

Julia が 1.6.0 になったのでインストールした後，いつものように atom で Julia プログラムをインクリメンタルに実行しようとして，今までと動作が違う（というか，ウインドウ構成も全く違ってしまった）のに気づいた。

あれこれやったせいで，何が原因かわからなくって，再インストールも何回かやって，丸一日ああだこうだやっていた。

やっと，理由がわかった。今までは Juno はインストールしていなかった。Hydrogen で逐次実行するとその結果がインラインで表示されていた。Juno をインストールしたので，動作が全く違ってしまったのだ。

uber-juno, julia-client, ink, toolbar を disable して，やっと今までの使い勝手にもどった。やれやれ。

Julia に翻訳--160 ED50，LD50

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

ED50 や LD50 の計算
http://aoki2.si.gunma-u.ac.jp/R/ed50.html

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

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

関数定義のどれか一つだけ有効になるように書いても，関数が再定義されたというウォーニングメッセージが出る。

==========#

using Rmath, Printf

function ed50(x, n, r; transform = "none") #  or "ln", or "log"
    dose = "x"
    if transform != "none"
        all(x .> 0) || error("all x_i must be positive")
        if transform == "ln"
            x = log.(x)
            dose = "ln(x)"
        elseif transform == "log"
            x = log10.(x)
            dose = "log10(x)"
        end
    end
    ollf, alpha, beta = estimation1(x, n, r)
    alpha, beta, ed50value, se, cll, clu =
        estimation2(x, n, r, alpha, beta; epsilon = 1e-13)
    println("P_hat = $alpha + $beta * $dose")
    println("ED50 = $(inversetransform(ed50value, transform))")
    cll = inversetransform(cll, transform)
    clu = inversetransform(clu, transform)
    println("95 % confidence interval = [$cll, $clu]")
            e, p0, p9, cllf, chi, df, p = GoodnessOfFitness(x, n, r, alpha, beta)
    @printf("%7s %7s %7s %10s %10s %9s\n", "dose", "n", "r", "e", "r/n", "e/n")
    for i = 1:length(x)
        @printf("%7.5g %7.0f %7.0f %10.5g %10.5f% 10.5f\n",
                x[i], n[i], r[i], e[i], p0[i], p9[i] )
    end
    println("chisq. = $chi,  df = $df,  p value = $p")
    println("　　観察値の対数尤度 = $ollf")
    println("あてはめ後の対数尤度 = $cllf")
end

function inversetransform(x, transform)
    if transform == "none"
        return x
    elseif transform == "ln"
        return exp.(x)
    else
        return 10 .^ x
    end
end

function estimation1(x, n, r)
    select = 0 .< r .< n
    x = x[select]
    n = n[select]
    r = r[select]
    p9 = r ./ n
    y = qnorm.(p9) .+ 5
    w = 1
    wx2 = sum(w * x .^ 2)
    wy = sum(w * y)
    wx = sum(w * x)
    wyx = sum(w * y .* x)
    ww = length(x)
    temp = ww * wx2 - wx^2
    alpha = (wx2 * wy - wx * wyx) / temp
    beta = (ww * wyx - wx * wy) / temp
    ollf = sum(r .* log.(p9) + (n - r) .* log.(1 .- p9))
    ollf, alpha, beta
end

function estimation2(x, n, r, alpha, beta; epsilon = 1e-5)
    sqrtpi2 = 1 / sqrt(2π)
    for ii = 1:100
        y = alpha .+ beta .* x
        x9 = y .- 5
        p9 = pnorm.(x9)
        z = exp.(-x9 .^ 2 / 2) * sqrtpi2
        y += (r ./ n .- p9) ./ z
        w = n .* z .^ 2 ./ ((1 .- p9) .* p9)
        wx2 = sum(w .* x .^ 2)
        wy = sum(w .* y)
        wx = sum(w .* x)
        wyx = sum(w .* y .* x)
        ww = sum(w)
        ssnwx = wx2 - wx^2 / ww
        beta1 = (wyx - wx * wy / ww) / ssnwx
        xbar = wx / ww
        alpha1 = (wy / ww) - beta1 * xbar
        if abs(alpha - alpha1) / alpha < epsilon &&
           abs(beta - beta1) / beta < epsilon
            ed50value = (5 - alpha1) / beta1
            se = sqrt(1 / ww + (ed50value - xbar)^2 / ssnwx) / abs(beta1)
            g = (qnorm(0.975) / beta1)^2 / ssnwx
            cl = ed50value + g * (ed50value - xbar) / (1 - g)
            pm = qnorm(0.975) / beta1 / (1 - g) *
                 sqrt((1 - g) / ww + (ed50value - xbar)^2 / ssnwx)
            return alpha1, beta1, ed50value, se, cl - pm, cl + pm
        end
        alpha = alpha1
        beta = beta1
    end
    error("収束しませんでした")
end

function GoodnessOfFitness(x, n, r, alpha, beta)
    p9 = pnorm.(alpha .+ beta .* x .- 5)
    p0 = r ./ n
    cllf = sum(r .* log.(p9) + (n - r) .* log.(1 .- p9))
    q = 1 .- p9
    chi = sum(n .* (p0 - p9) .^ 2 ./ (q - q .^ 2))
    df = length(x) - 2
    p = pchisq(chi, df, false)
    n .* p9, p0, p9, cllf, chi, df, p
end

x = [10, 20, 30];
n = [10, 10, 10];
r = [1, 5, 8];
ed50(x, n, r)
# ED50 = 21.181831322664955
# 95 % confidence interval = [15.207598835509728, 28.038741634031247]
#    dose       n       r          e        r/n       e/n
#      10      10       1     1.2081    0.10000   0.12081
#      20      10       5     4.5075    0.50000   0.45075
#      30      10       8     8.2211    0.80000   0.82211
# chisq. = 0.17215430107451019,  df = 1,  p value = 0.6782042663897273
# 　　観察値の対数尤度 = -15.186325774895813
# あてはめ後の対数尤度 = -15.272763997646091

ed50(x, n, r, transform = "ln")
# P_hat = -0.7457541595935755 + 1.9289230037245237 * ln(x)
# ED50 = 19.662964979570706
# 95 % confidence interval = [13.888400744959936, 27.976920343850338]
#    dose       n       r          e        r/n       e/n
#  2.3026      10       1    0.96075    0.10000   0.09608
#  2.9957      10       5     5.1308    0.50000   0.51308
#  3.4012      10       8     7.9243    0.80000   0.79243
# chisq. = 0.012098701206857332,  df = 1,  p value = 0.9124140570261099
# 　　観察値の対数尤度 = -15.186325774895813
# あてはめ後の対数尤度 = -15.1923793050997

ed50(x, n, r, transform = "log")
# P_hat = -0.7457541595936261 + 4.441509353909426 * log10(x)
# ED50 = 19.662964979570695
# 95 % confidence interval = [13.888400744959997, 27.976920343850143]
#    dose       n       r          e        r/n       e/n
#       1      10       1    0.96075    0.10000   0.09608
#   1.301      10       5     5.1308    0.50000   0.51308
#  1.4771      10       8     7.9243    0.80000   0.79243
# chisq. = 0.012098701206858081,  df = 1,  p value = 0.9124140570261072
# 　　観察値の対数尤度 = -15.186325774895813
# あてはめ後の対数尤度 = -15.192379305099697

Julia に翻訳--159 Log rank 検定

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

Log rank 検定
http://aoki2.si.gunma-u.ac.jp/R/logrank.html

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

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

NaN を除いて統計量を計算するには
1行関数で sumrmnan(x) = sum(x[.!isnan.(x)]) のように定義すると良い

==========#

using Rmath, Printf

function logrank(group, event, time; method = "SAS") # or method = "Tominaga"
    sumrmnan(x) = sum(x[.!isnan.(x)])

    len = length(group)
    length(event) == length(time) == len || error("length differ")
    indexx, indexy, tg = table(vcat(time, NaN, NaN, NaN, NaN),
          vcat(group, 1, 1, 2, 2) .* 10 .+ vcat(event, 1, 0, 1, 0));
    tg = tg[1:end-1, :]
    k = size(tg, 1)
    index, (nia, nib) = table(group)
    na = vcat(nia, (repeat([nia], k) - cumsum(tg[:, 1] + tg[:, 2]))[1:k - 1]);
    nb = vcat(nib, (repeat([nib], k) - cumsum(tg[:, 3] + tg[:, 4]))[1:k - 1]);
    da = tg[:, 2];
    db = tg[:, 4];
    dt = da + db;
    nt = na + nb;
    d = dt ./ nt;
    O = [sum(da), sum(db)];
    ea = na .* d;
    eb = nb .* d;
    E = [sumrmnan(ea), sumrmnan(eb)]
    @printf(" %4s  %3s  %3s  %3s  %3s  %3s  %3s  %10s  %10s  %10s\n",
            "time", "da", "db", "dt", "na", "nb", "nt", "d", "ea", "eb")
    for i = 1:k
        @printf("%5d %4d %4d %4d %4d %4d %4d %11.8f %11.8f %11.8f\n",
                indexx[i], da[i], db[i], dt[i], na[i], nb[i], nt[i],
                d[i], ea[i], eb[i])
    end
    if method == "Tominaga"
        method = "ログランク検定（富永）"
        chi = sum((O - E) .^ 2 ./ E)
    else
        method = "ログランク検定（一般的）"
        v = sumrmnan( dt .* (nt - dt) ./ (nt .- 1) .* na ./ nt .* (1 .- na ./ nt))
        chi = (sumrmnan(da) - sumrmnan(na .* d)) ^ 2 / v
    end
    P = pchisq(chi, 1, false)
    println(method)
    println("chisq. = $chi,  df = 1,  p value = $P")
    Dict(:chi => chi, :df => df, :pvalue => P, :method => method)
end

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

group = [1, 1, 2, 2, 1, 1, 2, 2, 1, 2, 2, 1, 1, 1, 2, 2, 2, 1, 1, 2, 1, 2, 1, 2, 1, 2, 2, 1, 1, 1];
event = [1, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0];
time = [2, 20, 5, 1, 3, 17, 2, 3, 15, 14, 12, 13, 11, 11, 10, 8, 8, 3, 7, 3, 6, 2, 5, 4, 2, 3, 1, 3, 2, 1];

logrank(group, event, time)
# time   da   db   dt   na   nb   nt           d          ea          eb
#    1    0    2    2   16   14   30  0.06666667  1.06666667  0.93333333
#    2    2    2    4   15   12   27  0.14814815  2.22222222  1.77777778
#    3    0    1    1   12   10   22  0.04545455  0.54545455  0.45454545
#    4    0    0    0    9    7   16  0.00000000  0.00000000  0.00000000
#    5    0    1    1    9    6   15  0.06666667  0.60000000  0.40000000
#    6    0    0    0    8    5   13  0.00000000  0.00000000  0.00000000
#    7    1    0    1    7    5   12  0.08333333  0.58333333  0.41666667
#    8    0    1    1    6    5   11  0.09090909  0.54545455  0.45454545
#   10    0    0    0    6    3    9  0.00000000  0.00000000  0.00000000
#   11    0    0    0    6    2    8  0.00000000  0.00000000  0.00000000
#   12    0    1    1    4    2    6  0.16666667  0.66666667  0.33333333
#   13    1    0    1    4    1    5  0.20000000  0.80000000  0.20000000
#   14    0    0    0    3    1    4  0.00000000  0.00000000  0.00000000
#   15    0    0    0    3    0    3  0.00000000  0.00000000  0.00000000
#   17    0    0    0    2    0    2  0.00000000  0.00000000  0.00000000
#   20    0    0    0    1    0    1  0.00000000  0.00000000  0.00000000
# ログランク検定（一般的）
# chisq. = 3.3805322873790673,  df = 1,  p value = 0.06597074594879886

logrank(group, event, time, method = "Tominaga")
# time   da   db   dt   na   nb   nt           d          ea          eb
#    1    0    2    2   16   14   30  0.06666667  1.06666667  0.93333333
#    2    2    2    4   15   12   27  0.14814815  2.22222222  1.77777778
#      略
#   20    0    0    0    1    0    1  0.00000000  0.00000000  0.00000000
# ログランク検定（富永）
# chisq. = 3.1527657452320845,  df = 1,  p value = 0.07579838765953799