Julia の小ネタ--013 下三角行列に要素を代入

下三角行列に要素を代入するとき，R では

n = 4
mat = matrix(0, n, n)
# A matrix: 4 × 4 of type dbl
# 0    0    0    0
# 0    0    0    0
# 0    0    0    0
# 0    0    0    0
x = c(0.1, 0.3, 0.6, 1.2, 3.1, 5.7)
mat[lower.tri(mat, diag=FALSE)] <- x
mat
# A matrix: 4 × 4 of type dbl
# 0.0    0.0    0.0    0
# 0.1    0.0    0.0    0
# 0.3    1.2    0.0    0
# 0.6    3.1    5.7    0

とするのが定道。

Julia では lower.tri のような関数はない（と思う）。

以下のようにするしかないか。

n = 4
mat = zeros(n, n)
# 4×4 Array{Float64,2}:
#  0.0  0.0  0.0  0.0
#  0.0  0.0  0.0  0.0
#  0.0  0.0  0.0  0.0
#  0.0  0.0  0.0  0.0
x = [0.1, 0.3, 0.6, 1.2, 3.1, 5.7] # length(x) == (n * (n - 1) ÷ 2)
print(collect(x)) # [0.1, 0.3, 0.6, 1.2, 3.1, 5.7]
lowertri = [i > j for i = 1:n, j = 1:n]
# 4×4 Array{Bool,2}:
#  0  0  0  0
#  1  0  0  0
#  1  1  0  0
#  1  1  1  0
mat[lowertri] = x;
mat
# 4×4 Array{Float64,2}:
#  0.0  0.0  0.0  0.0
#  0.1  0.0  0.0  0.0
#  0.3  1.2  0.0  0.0
#  0.6  3.1  5.7  0.0

lowertri = [i > j for i = 1:n, j = 1:n]

のように中間変数に代入しないようにすると

mat[[i > j for i = 1:n, j = 1:n]] = x

のようにも書けるから，[i > j for i = 1:n, j = 1:n] が lower.tri() に相当している。

 

Julia に翻訳--117 シェッフェの一対比較法

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

シェッフェの一対比較法
http://aoki2.si.gunma-u.ac.jp/R/ScheffePairedComparison.html

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

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

==========#

using Statistics, Combinatorics, Plots

function scheffepairedcomparison(A, B; digits=5)
    nr = size(A, 1)
    n = floor(Int, (1 + sqrt(1 + 8nr)) / 2)
    AB = A * B
    x = combinations(1:n, 2)
    nc = binomial(n, 2)
    indep = zeros(nc, n)
    for (i, (j,k)) in enumerate(x)
        indep[i, j] = 1
        indep[i, k] = -1
    end
    score = mreg(hcat(indep[:, 2:n], AB))
    sortedscore = sort(score)
    printscheffepairedcomparison(score, sortedscore, digits)
    plotscheffepairedcomparison(score, n)
    Dict(:score => score, :sortedscore => sortedscore)
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

function printscheffepairedcomparison(score, sortedscore, digits)
    println("スコア")
    println(round.(score, digits=digits))
    println("\nソートされたスコア")
    println(round.(sortedscore, digits=digits))
end

function plotscheffepairedcomparison(score, nc)
    pyplot(size=(600, 60))
    plt = scatter(score, repeat([25], nc), grid=false, tick_direction=:out,
            yshowaxis=false, ylims=(0, 60), yticks=false, label="")
    for i = 1:nc
        annotate!(score[i], 50, text(string(i), 10, :red), label="")
    end
    display(plt)
end

A = [
   10   13   41   33   10
    3   12   47   26   19
    2    9   32   12   52
   23   32   30   12   10
   27   11   31   13   25
   21    7   10   33   36
   ];

B = [4, 2, 0, - 2, - 4];

scheffepairedcomparison(A, B)
# スコア
# [73.0, 155.0, 129.5, 272.5]

# ソートされたスコア
# [73.0, 129.5, 155.0, 272.5]