Julia の小ネタ--020 べき乗のべき乗，大きさの見積，長精度

「べき乗のべき乗」の大きさの見積 であるが，Julia の big 型を使えば以下のようになる

注意点は big(9.8)  などとせず，big"9.8" のように文字型を使って定義するところ。

julia> x = big"4.5"
4.5

julia> y = big"6.7"
6.699999999999999999999999999999999999999999999999999999999999999999999999999972

julia> z = big"9.8"
9.800000000000000000000000000000000000000000000000000000000000000000000000000028

julia> x^y^z
8.141852291891192175294360390819933153170222320673963475352901403217459199966931e+81393094

julia> x = big"2.1"
2.099999999999999999999999999999999999999999999999999999999999999999999999999986

julia> y = big"7.2"
7.199999999999999999999999999999999999999999999999999999999999999999999999999972

julia> z = big"3.6"
3.599999999999999999999999999999999999999999999999999999999999999999999999999986

julia> x^y^z
1.384119872013753626430173206328153257960203418740675402484696742668262448719525e+393

 

Julia の小ネタ--019 空文字列に文字列をアペンド

空ベクトル a = [] に 文字列をアペンドしたい。

以下のようにしたのでは，期待した結果にならない。

a = []
append!(a, "aaaaa")
a
# 5-element Vector{Any}:
#  'a': ASCII/Unicode U+0061 (category Ll: Letter, lowercase)
#  'a': ASCII/Unicode U+0061 (category Ll: Letter, lowercase)
#  'a': ASCII/Unicode U+0061 (category Ll: Letter, lowercase)
#  'a': ASCII/Unicode U+0061 (category Ll: Letter, lowercase)
#  'a': ASCII/Unicode U+0061 (category Ll: Letter, lowercase)

こういう風にしないとだめ！

a = []
append!(a, ["aaaaa"])
a
# 1-element Vector{Any}:
#  "aaaaa"

Julia に翻訳--202 主座標分析，古典的多次元尺度構成法

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

主座標分析
http://aoki2.si.gunma-u.ac.jp/R/princo.html

ファイル名: princo.jl　　関数名: princo, printprinco, plotprinco, similaritymatrix

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

各種類似度行列を作る関数を整備すると面白いかも。

==========#

using LinearAlgebra, NamedArrays, Plots

function princo(s; objectnames=[])
  n = size(s, 1)
  length(objectnames) > 0 || (objectnames = ["Obj-$i" for i = 1:n])
  m = vec(sum(s, dims=1)) ./ n
  h = s .+ (sum(s) / n ^ 2) .- (m .+ m')
  values, vectors = eigen(h, sortby=x-> -x)
  values = values[values .> 1e-5]
  ax = length(values)
  vectors = sqrt.(values)' .* vectors[:, 1:ax]
  axisnames = ["Axis-$i" for i = 1:ax]
  colnames(vectors) = names(values) = paste("解", 1:ax, sep = "")
  rownames(vectors) = objectnames
  Dict(:ax => ax, :n => n, :values => values, :vectors => vectors,
       :objectnames => objectnames, :axisnames => axisnames)
end

function printprinco(res)
  val = res[:values]
  val2 = val / sum(val)
  println(NamedArray(hcat(val, val2, cumsum(val2))',
         (["eigen values", "contribution", "cumu. contribution"],
          res[:axisnames])))
  println("\nベクトル\n", NamedArray(res[:vectors],
         (res[:objectnames], res[:axisnames])))
end

function plotprinco(res; color=:black, annotate=false)
  pyplot()
  vectors = res[:vectors]
  size(vectors, 2) >= 2 || error("解が一次元なので，二次元配置図は描けません。")
  plt = scatter(vectors[:, 1], vectors[:, 2], grid=false, tick_direction=:out,
      color=color, markerstrokecolor=color, xlabel="Axis-1", ylabel="Axis-2", label="")
  annotate && annotate!(vectors[:, 1], vectors[:, 2],
                text.(" " .* res[:objectnames], 8, :blue, :left),
                label="")
  display(plt)
end

function similaritymatrix(dat)
  nr, nc = size(dat)
  s = zeros(nr, nr)
  for i = 1:nr-1
    for j = i+1:nr
      s[i, j] = s[j, i] = -sum((dat[i, :] .- dat[j, :]) .^ 2)
    end
  end
  s
end

s = [
  0 -1 -2 -3
  -1 0 -3 -4
  -2 -3 0 -1
  -3 -4 -1 0
];

res = princo(s)
printprinco(res)
#=====
3×3 Named Adjoint{Float64, Matrix{Float64}}
             A ╲ B │   Axis-1    Axis-2    Axis-3
───────────────────┼─────────────────────────────
eigen values       │  5.23607       1.0  0.763932
contribution       │  0.74801  0.142857  0.109133
cumu. contribution │  0.74801  0.890867       1.0

ベクトル
4×3 Named Matrix{Float64}
A ╲ B │    Axis-1     Axis-2     Axis-3
──────┼────────────────────────────────
Obj-1 │  0.850651       -0.5   0.525731
Obj-2 │   1.37638        0.5   -0.32492
Obj-3 │ -0.850651       -0.5  -0.525731
Obj-4 │  -1.37638        0.5    0.32492
=====#

plotprinco(res)

ユークリッド平方距離を類似度に変換して分析する

using RDatasets
iris = dataset("datasets", "iris");
s = similaritymatrix(Matrix(iris[:, 1:4]))
res = princo(s, objectnames=string.(collect(1:150)))
printprinco(res)
plotprinco(res, color=vcat(repeat([:black], 50), repeat([:blue], 50), repeat([:red], 50)...))

savefig("fig1.png")

Julia に翻訳--201 順位データの双対尺度法

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

順位データの双対尺度法
http://aoki2.si.gunma-u.ac.jp/R/ro.dual.html

ファイル名: rodual.jl　　関数名: rodual, printdual, plotdual

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

==========#

using LinearAlgebra, NamedArrays, Plots

function rodual(F; rownames=[], colnames=[])
    N, n = size(F)
    length(rownames) > 0 || (rownames = ["Row-$i" for i = 1:N])
    length(colnames) > 0 || (colnames = ["Col-$i" for i = 1:n])
    E = (n + 1) .- 2 .* F
    Hn = (transpose(E) * E) ./ (N * n * (n - 1) ^ 2)
    values, vectors = eigen(Hn, sortby=x-> -x)
    ne = size(Hn, 1) - 1
    eta2 = values[1:ne]
    eta = sqrt.(eta2)
    contribution = eta2 ./ sum(values[1:ne]) .* 100
    cumcont = cumsum(contribution)
    axisnames = ["Axis-$i" for i = 1:ne]
    W = vectors[:, 1:ne]
    colscore = W .* sqrt(n)
    colscore2 = eta' .* colscore
    rowscore2 = (E * W) ./ (sqrt(n) * (n - 1))
    rowscore = rowscore2 ./ eta'
    Dict(:eta2 => eta2, :eta => eta, :contribution => contribution,
         :cumcont => cumcont, :colscore => colscore,
         :colscore2 => colscore2, :rowscore => rowscore,
         :rowscore2 => rowscore2, :rownames => rownames,
         :colnames => colnames, :axisnames => axisnames)
end

function printdual(x; weighted = true)
    println(NamedArray(hcat(x[:eta2], x[:eta], x[:contribution], x[:cumcont])',
        (["eta square", "correlation", "contribution", "cumulative contribution"],
        x[:axisnames])))
    if weighted
        println("\nweighted row score\n",
            NamedArray(x[:rowscore2], (x[:rownames], x[:axisnames])))
        println("\nweighted col score\n",
            NamedArray(x[:colscore2], (x[:colnames], x[:axisnames])))
    else
        println("\nrow score\n",
            NamedArray(x[:rowscore], (x[:rownames], x[:axisnames])))
        println("\ncol score\n",
            NamedArray(x[:colscore], (x[:colnames], x[:axisnames])))
    end
end

function plotdual(x; weighted=true, ax1=1, ax2=2)
    pyplot()
    if weighted
        row, col = x[:rowscore], x[:colscore]
    else
        row, col = x[:rowscore2], x[:colscore2]
    end
    plt = scatter(row[:, ax1], row[:, ax2], grid=false, tick_direction=:out,
                  color=:red, xlabel="Axis-$ax1", ylabel="Axis-$ax2",
                  aspect_ratio=1, size=(600, 600), label="")
    annotate!.(row[:, ax1], row[:, ax2],
               text.(" " .* x[:rownames], 8, :red, :left));
    scatter!(col[:, ax1], col[:, ax2], grid=false, tick_direction=:out,
            color=:blue, label="")
    annotate!.(col[:, ax1], col[:, ax2],
               text.(" " .* x[:colnames], 8, :blue, :left))
    display(plt)
end

F = [6 1 5 3 2 8 4 7
        3 8 1 6 7 5 4 2
        5 7 1 6 8 2 4 3
        4 6 2 3 8 7 1 5
        2 4 6 3 7 5 1 8
        2 4 5 3 8 7 1 6
        1 7 6 3 8 5 2 4
        7 5 3 1 8 4 6 2
        4 2 7 3 8 6 5 1
        5 1 2 4 7 6 3 8
        6 4 3 2 8 7 5 1
        3 8 4 2 5 6 1 7
        3 2 1 6 4 7 5 8
        5 8 1 4 7 3 6 2]

res = rodual(F)
printdual(res)
#=====
4×7 Named Adjoint{Float64, Matrix{Float64}}
                  A ╲ B │     Axis-1      Axis-2  …      Axis-6      Axis-7
────────────────────────┼──────────────────────────────────────────────────
eta square              │   0.141472    0.122668  …   0.0160181  0.00491559
correlation             │   0.376128     0.35024       0.126562   0.0701113
contribution            │    33.0102     28.6226        3.73755     1.14697
cumulative contribution │    33.0102     61.6328  …      98.853       100.0

weighted row score
14×7 Named Matrix{Float64}
 A ╲ B │      Axis-1       Axis-2  …       Axis-6       Axis-7
───────┼──────────────────────────────────────────────────────
Row-1  │   -0.407444    -0.348416  …    0.0933572    0.0370619
Row-2  │      0.5271     0.153163        0.222556    0.0255489
Row-3  │    0.491104     0.240422       -0.157631    0.0962444
Row-4  │    0.492153    -0.370614       0.0611955     0.108397
Row-5  │    0.129472    -0.584181       -0.174652   -0.0338449
Row-6  │    0.301451    -0.562441        0.031781     0.014255
Row-7  │    0.460194    -0.277195       0.0445545   -0.0775488
Row-8  │      0.4095     0.206126       -0.177598    -0.086752
Row-9  │     0.20572   -0.0267503       0.0698013  -0.00513034
Row-10 │   0.0203491    -0.446696       -0.184627    0.0212436
Row-11 │    0.415701    0.0522849        0.119628    0.0172468
Row-12 │    0.269953    -0.411087       0.0304552   -0.0294498
Row-13 │   -0.174422    -0.314766        0.132588    -0.137881
Row-14 │     0.49742     0.359662  …   -0.0244221    -0.101179

weighted col score
8×7 Named Matrix{Float64}
A ╲ B │      Axis-1       Axis-2  …       Axis-6       Axis-7
──────┼──────────────────────────────────────────────────────
Col-1 │    0.131864    -0.293808  …     0.121021    -0.103093
Col-2 │   -0.472329    -0.232793      -0.0557051    0.0137367
Col-3 │    0.325684    0.0866956       0.0296036   -0.0148614
Col-4 │    0.212309    -0.209285      -0.0738725   -0.0830142
Col-5 │   -0.753761     0.178061        0.115497   0.00743009
Col-6 │  -0.0154598     0.415622       -0.262662  -0.00767548
Col-7 │    0.221272    -0.512397      -0.0199407     0.136794
Col-8 │    0.350421     0.567905  …     0.146059     0.050683
=====#

printdual(res, weighted = false)
#=====
4×7 Named Adjoint{Float64, Matrix{Float64}}
                  A ╲ B │     Axis-1      Axis-2  …      Axis-6      Axis-7
────────────────────────┼──────────────────────────────────────────────────
eta square              │   0.141472    0.122668  …   0.0160181  0.00491559
correlation             │   0.376128     0.35024       0.126562   0.0701113
contribution            │    33.0102     28.6226        3.73755     1.14697
cumulative contribution │    33.0102     61.6328  …      98.853       100.0

row score
14×7 Named Matrix{Float64}
 A ╲ B │     Axis-1      Axis-2  …      Axis-6      Axis-7
───────┼──────────────────────────────────────────────────
Row-1  │   -1.08326   -0.994792  …    0.737638    0.528616
Row-2  │    1.40138    0.437308        1.75847    0.364405
Row-3  │    1.30568    0.686448       -1.24548     1.37274
Row-4  │    1.30847    -1.05817        0.48352     1.54607
Row-5  │   0.344224    -1.66794       -1.37997   -0.482731
Row-6  │   0.801457    -1.60587        0.25111    0.203319
Row-7  │     1.2235   -0.791443       0.352036    -1.10608
Row-8  │    1.08873    0.588526       -1.40324    -1.23735
Row-9  │    0.54694  -0.0763771       0.551516  -0.0731742
Row-10 │  0.0541015     -1.2754       -1.45878    0.302999
Row-11 │    1.10521    0.149283       0.945207    0.245991
Row-12 │   0.717716    -1.17373       0.240634   -0.420043
Row-13 │  -0.463729   -0.898715        1.04761     -1.9666
Row-14 │    1.32248      1.0269  …   -0.192965    -1.44312

col score
8×7 Named Matrix{Float64}
A ╲ B │     Axis-1      Axis-2  …      Axis-6      Axis-7
──────┼──────────────────────────────────────────────────
Col-1 │   0.350584   -0.838875  …    0.956215    -1.47042
Col-2 │   -1.25577   -0.664667      -0.440139    0.195927
Col-3 │   0.865885    0.247532       0.233905   -0.211968
Col-4 │   0.564459   -0.597548      -0.583684    -1.18403
Col-5 │     -2.004    0.508396       0.912572    0.105976
Col-6 │ -0.0411024     1.18668       -2.07536   -0.109476
Col-7 │    0.58829    -1.46299      -0.157556      1.9511
Col-8 │   0.931652     1.62147  …     1.15404    0.722894
=====#

plotdual(res)
savefig("fig1.png")