Julia に翻訳--155

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

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

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

一回しか使わないであろうプログラムも関数にしよう。
条件設定はパラメータとして与える。
啓蒙目的ならやむを得ないが，なるべく既存の関数を使おう。

対象読者もそのようなプログラムを書くことができるるようになることを目指そう。

==========#

using Statistics, Rmath

function samplesizedeterminationforcorrelationcoefficient(x, y; mu=0, alpha=0.05, side=2, power=0.8)
    mdx = x .- mean(x) # mean deviation
    mdy = y .- mean(y)
    r = sum(mdx .* mdy) / (sqrt(sum(mdx .^ 2)) * sqrt(sum(mdy .^ 2))) # cor(x, y)
    n = length(x)
    t = (abs(r) - mu) / sqrt((1 - r ^ 2) / (n - 2))
    p_t = pt(t, n - 2, false) * 2
    Za = qnorm(alpha / side, false)
    Zb = qnorm(power)
    z = atanh(r) # 1 / 2 *  log((1 + r) / (1 - r))
    n = ((Za + Zb) / z) ^ 2 + 3
    Dict(:n => n, :r => r, :t => t, :p_t => p_t, :alpha => alpha, :side => side, :power => power, :n => n)
end

chukan = [64, 40, 71, 33, 30, 71, 92, 23, 41, 55, 93, 74];
kimatu = [55, 52, 76, 24, 48, 87, 100, 30, 35, 67, 86, 81];

samplesizedeterminationforcorrelationcoefficient(chukan, kimatu)
#========================
Dict{Symbol,Real} with 7 entries:
  :alpha => 0.05
  :n     => 6.12334
  :p_t   => 2.33485e-5
  :power => 0.8
  :r     => 0.919416
  :t     => 7.39269
  :side  => 2
========================#

Julia に翻訳--154 Kaplan-Meier 法による生命表

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

Kaplan-Meier 法による生命表
http://aoki2.si.gunma-u.ac.jp/R/km-surv.html

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

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

==========#

using Rmath, Printf, Plots

function kmsurv(time, event)
    str(x) = isnan(x) ? "" : @sprintf("%9.7f", x)
    fraction = minimum(time[time .> 0]) / 1000
    n = length(time)
    truncate = 1 .- event;
    time += truncate .* fraction;
    o = sortperm(time);
    time = time[o];
    truncate = truncate[o]
    time -= truncate .* fraction
    p = [truncate[i] == 1 ? 1 : (n - i) / (n - i + 1) for i = 1:n];
    P = cumprod(p);
    se = [(1 - truncate[i]) / (n - i + 1) / (n - i) for i = 1:n];
    temp = P .* sqrt.(cumsum(se))
    SE = [truncate[i] == 0 ? temp[i] : NaN for i = 1:n]
    @printf("%4s %8s  %8s  %8s  %8s\n", "time", "truncate", "p", "P", "SE")
    for i = 1:length(time)
        @printf("%4d %8d %9.7f %9.7f %9s\n",
                time[i], truncate[i], p[i], P[i], str(SE[i]))
    end
    time = vcat(0, time)
    P = vcat(1, P)
    pyplot()
    plt = plot(time, P, seriestype=:steppost, grid=false, tick_direction=:out,
               ylims=(-0.01, 1.01), ylabel="P", xlabel="time", label="")
    scatter!(time, P, label="")
    display(plt)
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];
agroup = group .== 1;
kmsurv(time[agroup], event[agroup])

time truncate         p         P        SE
   1        1 1.0000000 1.0000000          
   2        0 0.9333333 0.9333333 0.0644061
   2        0 0.9285714 0.8666667 0.0877707
   2        1 1.0000000 0.8666667          
   3        1 1.0000000 0.8666667          
   3        1 1.0000000 0.8666667          
   3        1 1.0000000 0.8666667          
   5        1 1.0000000 0.8666667          
   6        1 1.0000000 0.8666667          
   7        0 0.8571429 0.7428571 0.1371088
  11        1 1.0000000 0.7428571          
  11        1 1.0000000 0.7428571          
  13        0 0.7500000 0.5571429 0.1908971
  15        1 1.0000000 0.5571429          
  17        1 1.0000000 0.5571429          
  20        1 1.0000000 0.5571429          

Julia に翻訳--153 バートレットの球面性検定

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

バートレットの球面性検定
http://aoki2.si.gunma-u.ac.jp/R/Bartlett.sphericity.test.html

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

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

==========#

using Statistics, Rmath, LinearAlgebra

function bartlettsphericitytest(x)
   n, p = size(x)
   chisq = (1 - n + (2 * p + 5) / 6) * log(det(cor(x)))
   df = p * (p - 1) ÷ 2
   pvalue = pchisq(chisq, df, false)
   println("chisq. = $chisq,  df = $df,  p value = $pvalue")
   Dict(:chisq => chisq, :df => df, :pvalue => pvalue)
end

x = [ 1  5 6 4
      2 14 5 3
      3  3 4 2
      4  2 6 6
      3  4 3 5 ]

bartlettsphericitytest(x)

# chisq. = 1.361750057342831,  df = 6,  p value = 0.9681438167949121

Julia に翻訳--152 多重共線性のチェック，従属性

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

多重共線性のチェック（従属性）
http://aoki2.si.gunma-u.ac.jp/R/find_multico.html

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

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

R だと数行なのに Julia でちゃんと書くと長くなる。

==========#

using Statistics, LinearAlgebra, Printf

function findmultico(data, epsilon = 1e-10)
    values, vectors = eigen(cor(data))
    n = length(values)
    result = fill("", n, n)
    exists = false
    for i = 1:n
        if values[i] < epsilon
            for j = 1:n
                if abs(vectors[j, i]) > epsilon
                    result[j, i] = "***"
                    exists = true
                end
            end
        end
    end
    if exists
        for j = 1:n
            @printf("%7s", "Var"*string(j))
            for i = 1:n
                @printf("%5s", result[j, i])
            end
            println()
        end
    else
        println("no multico")
    end
end

x = [ 1 2 4 3
      3 2 5 5
      4 3 7 7
      2 1 3 3
      5 4 7 9 ]

findmultico(x)
# Var1  ***               
# Var2  ***               
# Var3                    
# Var4  ***               

 

Julia に翻訳--151 多重共線性のチェック，トレランス

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

多重共線性のチェック（トレランス）
http://aoki2.si.gunma-u.ac.jp/R/tolerance.html

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

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

==========#

using Statistics, LinearAlgebra, Printf

function tolerance(x)
    m = size(x, 2)
    VIF = [inv(cor(x))[i, i] for i = 1:m]
    tolerance = 1 ./ VIF
    @printf("       %12s  %12s\n", "tolerance", "VIF")
    for i = 1:m
        @printf("%5s  %12.7f  %12.7f\n", "Var"*string(i), tolerance[i], VIF[i])
    end
end

x = [ 1 2 4
      3 2 5
      4 3 7
      2 1 3
      5 4 7 ]

tolerance(x)
#          tolerance           VIF
# Var1     0.2181818     4.5833333
# Var2     0.1318681     7.5833333
# Var3     0.0937500    10.6666667

 

Julia に翻訳--150 サンプリング適切性基準

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

サンプリング適切性基準
http://aoki2.si.gunma-u.ac.jp/R/kmo.html

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

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

diag() 欲しいな。

==========#

using Statistics, LinearAlgebra, Rmath

function kmo(x)
    m = size(x, 2)
    r = cor(x)
    r2 = r .^ 2
    i = inv(r)
    d = [i[k, k] for k =1:m]
    p2 = (-i ./ sqrt.(d * d')) .^ 2
    [r2[k, k] = p2[k, k] = 0 for k = 1:m]
    KMO = sum(r2) / (sum(r2) + sum(p2))
    MSA = sum(r2, dims=1) ./ (sum(r2, dims=1) + sum(p2, dims=1))
    println("KMO = $KMO")
    println("MSA = $MSA")
    Dict(:KMO => KMO, :MSA => MSA)
end

x = [
    1  5 6 4
    2 14 5 3
    3  3 4 2
    4  2 6 6
    3  4 3 5]

kmo(x)
# KMO = 0.5396372472781039
# MSA = [0.5499697817773649 0.6870648139657886 0.30803672198983056 0.5320125096821899]