空ベクトル 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 の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを 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 の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを 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")
#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。
クロス集計表・分布表の双対尺度法
http://aoki2.si.gunma-u.ac.jp/R/dual.html
ファイル名: dual.jl 関数名: dual, printdual, plotdual
翻訳するときに書いたメモ
==========#
using Rmath, LinearAlgebra, Printf, NamedArrays, Plots
function dual(tbl; rownames=[], colnames=[])
nr, nc = size(tbl)
length(rownames) > 0 || (rownames = ["Row-$i" for i = 1:nr])
length(colnames) > 0 || (colnames = ["Col-$i" for i = 1:nc])
ft = sum(tbl)
fr = sum(tbl, dims=2)
fc = vec(sum(tbl, dims=1))
temp = sqrt.(fc * fc')
w = transpose(tbl ./ fr) * tbl ./ temp .- temp ./ ft
values, vectors = eigen(w, sortby=x-> -x)
ne = sum(values .> 1e-5)
vectors = vectors[:, 1:ne]
values = values[1:ne]
colscore = vectors .* sqrt.(ft ./ fc)
rowscore = tbl * colscore ./ (fr .* sqrt.(values)')
colscore2 = transpose(sqrt.(values)) .* colscore
rowscore2 = transpose(sqrt.(values)) .* rowscore
cont = values ./ tr(w) * 100
chisq = -(ft - 1 - (nr + nc - 1) / 2) .* log.(1 .- values)
df = (nr + nc - 1) .- 2 .* (1:ne)
P = pchisq.(chisq, df, false)
axisnames = ["Axis-$i" for i = 1:ne]
Dict(:values => values, :correlation => sqrt.(values), :cont => cont,
:cumulativecontribution => cumsum(cont), :chisq => chisq, :df => df,
:P => P, :rowscore => rowscore, :colscore => colscore,
:rowscore2 => rowscore2, :colscore2 => colscore2,
:rownames => rownames, :colnames => colnames, :axisnames => axisnames)
end
function printout(name, vector; str=false)
@printf("%-25s", name)
for i = 1:length(vector)
str ? @printf("%9s", vector[i]) : @printf("%9.3f", vector[i])
end
println()
end
function printdual(x; weighted = true)
printout("", x[:axisnames], str=true)
printout("eta square", x[:values])
printout("correlation", sqrt.(x[:values]))
printout("contribution", x[:cont])
printout("cumulative contribution", cumsum(x[:cont]))
printout("Chi square value", x[:chisq])
printout("degrees of freedom", x[:df], str=true)
printout("P value", x[:P])
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(results; weighted=true, ax1=1, ax2=2)
pyplot()
if weighted
row, col = results[:rowscore], results[:colscore]
else
row, col = results[:rowscore2], results[: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.(" " .* results[:rownames], 8, :red, :left));
scatter!(col[:, ax1], col[:, ax2], grid=false, tick_direction=:out,
color=:blue, label="")
annotate!.(col[:, ax1], col[:, ax2],
text.(" " .* results[:colnames], 8, :blue, :left))
display(plt)
end
tbl = [
2 3 5 6
5 1 7 5
5 3 4 3
]
res = dual(tbl)
printdual(res)
#=====
Axis-1 Axis-2
eta square 0.049 0.038
correlation 0.221 0.194
contribution 56.572 43.428
cumulative contribution 56.572 100.000
Chi square value 2.263 1.727
degrees of freedom 4 2
P value 0.687 0.422
weighted row score
3×2 Named Matrix{Float64}
A ╲ B │ Axis-1 Axis-2
──────┼─────────────────────────
Row-1 │ -0.317955 -0.00775345
Row-2 │ 0.147306 0.219535
Row-3 │ 0.162385 -0.255172
weighted col score
4×2 Named Matrix{Float64}
A ╲ B │ Axis-1 Axis-2
──────┼───────────────────────
Col-1 │ 0.343349 -0.0831775
Col-2 │ -0.206018 -0.419061
Col-3 │ 0.0256527 0.153724
Col-4 │ -0.220608 0.10514
=====#
plotdual(res) # savefig("fig1.png")
plotdual(res, weighted=false)
これは，々とおなじで「おなじ」といれて変換する。
他に「おなじ」は，ゝ，ゞ，ヽ，ヾ，仝 がある。
これは簡単。
「これ」っていれて，変換する。名前によく使われるので「ゆき」のほうがよく使われているかも。「し」でもよいけど，これは場合によっては変換キーを何回も押さなくてはならなくて，行きすぎたりする。
石碑などの刻印に，「○○建之」とあるのは「○○がこれを建てた」ということ。
文字列の連結 3 通り
greet = "Hello"
whom = "world"
から，"Hello, world.\n" を作る
第1の方法
string(greet, ", ", whom, ".\n")
第2の方法
greet * ", " * whom * ".\n"
第3の方法
"$greet, $whom.\n"
#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。
因子分析の適合度検定
http://aoki2.si.gunma-u.ac.jp/R/fa.fit.test.html
ファイル名: fafittest.jl 関数名: fafittest
翻訳するときに書いたメモ
計算の途中で結果が Matrix{Any} になることがある。今回の場合は
sc = [i == j ? 1 / sqrt(uniquenesses[i]) : 0 for i = 1:p, j = 1:p]
では，sc が Matrix{Real} になるのが原因。それを回避するには
sc = [i == j ? 1.0 / sqrt(uniquenesses[i]) : 0.0 for i = 1:p, j = 1:p]
とすればよいようだが，なんだかなあ！
==========#
using DataFrames, Statistics, StatsBase, Rmath, LinearAlgebra
function fafittest(data, factors, uniquenesses)
n, p = size(data)
r = cor(data)
sc = [i == j ? 1.0 / sqrt(uniquenesses[i]) : 0.0 for i = 1:p, j = 1:p]
r = sc * r * sc
values = eigvals(r) # Julia ではデフォルトで昇順
e = values[1:p - factors]
s = -sum(log.(e) .- e) + factors - p
chisq = (n - 1 - (2 * p + 5) / 6 - (2 * factors) / 3) * s
df = ((p - factors) ^ 2 - p - factors) / 2
pvalue = pchisq(chisq, df, false)
println("chisq. = $chisq, df = $df, p value = $pvalue")
Dict(:chisq => chisq, :df => df, :pvalue => pvalue)
end
data = [
935 955 926 585 1010 925 1028 807 769 767
817 905 901 632 1004 950 957 844 781 738
768 825 859 662 893 900 981 759 868 732
869 915 856 448 867 874 884 802 804 857
787 878 880 592 871 874 884 781 782 807
738 848 850 569 814 950 957 700 870 764
763 862 839 658 887 900 1005 604 709 753
795 890 841 670 853 874 859 701 680 772
903 877 919 460 818 849 884 700 718 716
761 765 881 485 846 900 981 728 781 714
747 792 800 564 796 849 932 682 746 767
771 802 840 609 824 874 859 668 704 710
785 878 805 527 911 680 884 728 709 747
657 773 820 612 810 849 909 698 746 771
696 785 791 578 774 725 932 765 706 795
724 785 870 509 746 849 807 763 724 760
712 829 838 516 875 725 807 754 762 585
756 863 815 474 873 725 957 624 655 620
622 759 786 619 820 769 807 673 698 695
668 753 751 551 834 849 807 601 655 642];
uniquenesses = [0.005, 0.186849474704386, 0.235609633298087,
0.594548313805094, 0.005, 0.005, 0.644059344661751,
0.005, 0.466720965448361, 0.573892982555005]
fafittest(data, 4, uniquenesses)
# chisq. = 15.163572526523414, df = 11.0, p value = 0.1751303147464357
どこかでだれかが，「演算子 \ なんて使い道あるのか？」なんて書いていたが，ありますよ。
a*x = b のとき x を求める，つまり x = b /a = 1/a * b という記述を，\ を使って x = a \ b と記述できるということですよ。つまり，a の逆数と b の積を意味するんですね。
julia> a = 3
3
julia> b = 7
7
julia> x = b / a
2.3333333333333335
julia> x = a \ b
2.3333333333333335
これだけだとあまりありがたみがないんだけど，a や b が配列やベクトルであっても同じように使えますよ。
A * x = B のとき，x は A の逆行列と B の行列積ですね。上と同じように A の逆数と B の積で x = A^(-1) * B = inv(A) * B です。\ を使えば A \ B な訳です。
julia> A = [2 3; 5 7]
2×2 Matrix{Int64}:
2 3
5 7
julia> B = [11, 13]
2-element Vector{Int64}:
11
13
julia> inv(A) * B
2-element Vector{Float64}:
-38.00000000000006
29.00000000000004
julia> A \ B
2-element Vector{Float64}:
-38.00000000000006
29.00000000000004
で，同じになることがわかります。
じゃあ，わざわざ新しい演算子なんか持ち出すことはないだろう，ドッチでもよいじゃないか。
いやいや，そうじゃないんです。
inv(A) * B は A \ B より，3倍ぐらい遅いんです。
とは言っても，8000x8000 の行列 A と要素数 8000 のベクトル B の演算だと，inv(A) * B は 18 秒，A \ B は 5 秒ぐらいなので，五十歩百歩ではありますが。
R ではそれぞれ solve(A)*B ，solve(A, B) に対応するのだけど，同じ条件で，前者は 550 秒，後者は 125 秒ぐらいなので，Julia は相当に速いということではあります。
今日の「日本人の3割しか知らないこと」で，「おなじ」と入力して変換するというのが解であったが，
「どう」と，二文字の入力して変換するのでも入力できる。
まあ，初期状態で何回変換キーを押せばよいかの違いはあるだろうが。一度入力すれば学習されるので次回以降も「どう」を入力するのが吉。
だいぶ前に知っていたが，なつかしい。
ちなみに，だいぶ前，私の愛車のナビではどうやっても入力できなかった記憶がある。今はどうなっているのかな？
#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。
正準相関分析
http://aoki2.si.gunma-u.ac.jp/R/cancor.html
ファイル名: cancor.jl 関数名: cancor
翻訳するときに書いたメモ
==========#
using Statistics, StatsBase, LinearAlgebra, Rmath, NamedArrays
function cancor(x, gr1, gr2)
function geneig2(a, b, k, sd)
if size(a, 1) == 1
values = a / b
vectors = [1]
else
values, vectors = eigen(b, sortby=x-> -x)
g = [i == j ? 1 / sqrt(values[i]) : 0 for i = 1:k, j = 1:k]
values, vectors2 = eigen(float.(g * transpose(vectors) * a * vectors * g), sortby=x-> -x)
vectors = vectors * g * vectors2
end
stdvectors = vectors[:, 1:k]
unstdvectors = stdvectors ./ sd
values[1:k], stdvectors, unstdvectors
end
k = min(length(gr1), length(gr2))
r = cor(x)
S11 = r[gr1, gr1]
S22 = r[gr2, gr2]
S12 = r[gr1, gr2]
x1 = x[:, gr1]
x2 = x[:, gr2]
sd1 = std(x1, dims=1)
sd2 = std(x2, dims=1)
values1, stdvectors1, unstdvectors1 = geneig2(S12 * (S22 \ S12'), S11, k, sd1)
values2, stdvectors2, unstdvectors2 = geneig2(transpose(S12) * (S11 \ S12), S22, k, sd2)
score1 = scale(x1 * unstdvectors1)
score2 = scale(x2 * unstdvectors2)
Dict(:canonicalcorrelationcoefficients => sqrt.(values1),
:standardizedcoefficientsgroup1 => stdvectors1,
:standardizedcoefficientsgroup2 => stdvectors2,
:coefficientsgroup1 => unstdvectors1,
:coefficientsgroup2 => unstdvectors2,
:canonicalscoresgroup1 => score1,
:canonicalscoresgroup2 => score2)
end
function scale(x)
x = float.(x)
m = mean(x, dims=1)
s = std(x, dims=1)
for i = 1:length(m)
x[:, i] .= (x[:, i] .- m[i]) ./ s[i]
end
x
end
function printcancor(obj)
println("\n正準相関係数 = $(obj[:canonicalcorrelationcoefficients])")
println("\n第1群の変数に対する標準化された係数\n",
NamedArray(obj[:standardizedcoefficientsgroup1]))
println("\n第2群の変数に対する標準化された係数\n",
NamedArray(obj[:standardizedcoefficientsgroup2]))
println("\n第1群の変数に対する係数\n",
NamedArray(obj[:coefficientsgroup1]))
println("\n第2群の変数に対する係数\n",
NamedArray(obj[:coefficientsgroup2]))
println("\n第1群の変数に対する正準得点\n",
NamedArray(obj[:canonicalscoresgroup1]))
println("\n第2群の変数に対する正準得点\n",
NamedArray(obj[:canonicalscoresgroup2]))
end
x = [
2 1 2 2
1 2 -1 -1
0 0 0 0
-1 -2 -2 1
-2 -1 1 -2];
gr1 = 1:2
gr2 = 3:4
obj = cancor(x, 1:2, 3:4)
printcancor(obj)
#=====
正準相関係数 = [0.9999999999999992, 0.30151134457776346]
第1群の変数に対する標準化された係数
2×2 Named Matrix{Any}
A ╲ B │ 1 2
──────┼─────────────────────────
1 │ 1.66667 1.11022e-16
2 │ -1.33333 -1.0
第2群の変数に対する標準化された係数
2×2 Named Matrix{Any}
A ╲ B │ 1 2
──────┼───────────────────────────
1 │ -2.77556e-16 -1.00504
2 │ 1.0 0.100504
第1群の変数に対する係数
2×2 Named Matrix{Float64}
A ╲ B │ 1 2
──────┼─────────────────────────
1 │ 1.05409 7.02167e-17
2 │ -0.843274 -0.632456
第2群の変数に対する係数
2×2 Named Matrix{Float64}
A ╲ B │ 1 2
──────┼───────────────────────────
1 │ -1.75542e-16 -0.635642
2 │ 0.632456 0.0635642
第1群の変数に対する正準得点
5×2 Named Matrix{Float64}
A ╲ B │ 1 2
──────┼─────────────────────
1 │ 1.26491 -0.632456
2 │ -0.632456 -1.26491
3 │ 0.0 0.0
4 │ 0.632456 1.26491
5 │ -1.26491 0.632456
第2群の変数に対する正準得点
5×2 Named Matrix{Float64}
A ╲ B │ 1 2
──────┼─────────────────────
1 │ 1.26491 -1.14416
2 │ -0.632456 0.572078
3 │ 0.0 0.0
4 │ 0.632456 1.33485
5 │ -1.26491 -0.76277
=====#
#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。
正準判別分析
http://aoki2.si.gunma-u.ac.jp/R/candis.html
ファイル名: candis.jl 関数名: candis
翻訳するときに書いたメモ
後の利用のために，結果を Dict で返すようにしているが，Dict が長すぎるか？？
==========#
using Statistics, StatsBase, Rmath, LinearAlgebra, NamedArrays, Plots, StatsPlots
function candis(data, group)
vnames = names(data);
data = Matrix(data);
n, p = size(data)
gname, ni = table(group)
k = length(ni)
grandmeans = vec(mean(data, dims=1));
groupmeans = zeros(p, k);
groupvars = zeros(p, k);
w = zeros(p, p);
for (i, g) in enumerate(gname)
data2 = data[group .== g, :]
groupmeans[:, i] = mean(data2, dims=1)
groupvars[:, i] = var(data2, dims=1)
w .+= cov(data2, corrected=false) * ni[i]
end
means = hcat(grandmeans, groupmeans);
b = (cov(data, corrected=false) .* n) .- w
df1, df2 = k - 1, n - k
Fvalue = [onewaytest(ni, groupmeans[i, :], groupvars[i, :]) for i = 1:p]
pvalue = pf.(Fvalue, df1, df2, false)
wilksFromF = df2 ./ (Fvalue .* df1 .+ df2)
ss = cov(data, corrected=false) .* n .- w
withincov = w ./ (n - k)
sd = sqrt.(diag(withincov))
withinr = withincov ./ (sd * sd')
eigenvalues, eigenvectors = eigen(b, w, sortby=x-> -x)
nax = sum(eigenvalues .> 1e-10)
axnames = "axis " .* string.(1:nax)
eigenvalues = eigenvalues[1:nax]
eigenvectors = eigenvectors[:, 1:nax]
Λ = reverse(cumprod(1 ./ (1 .+ reverse(eigenvalues))))
chisq = ((p + k) / 2 - n + 1) .* log.(Λ)
lL = (1:nax) .- 1
df = (p .- lL) .* (k .- lL .- 1)
pwilks = pchisq.(chisq, df, false)
canonicalcorrcoef = sqrt.(eigenvalues ./ (1 .+ eigenvalues))
temp = diag(transpose(eigenvectors) * w * eigenvectors)
temp = 1 ./ sqrt.(temp ./ (n - k))
coeff = eigenvectors .* temp'
const0 = - grandmeans' * coeff
coefficient = cat(coeff, const0, dims=1)
stdcoefficient = coeff .* sqrt.(diag(withincov));
centroids = groupmeans' * coeff .+ const0;
canscore = data * coeff .+ const0;
structure = zeros(p, nax);
for i = 1:nax
data2 = hcat(canscore[:, i], data);
ss = zeros(p + 1, p + 1)
for j = 1:k
dss = data2[group .== gname[j], :]
ss .+= cov(dss, corrected=false) * ni[j]
structure[:, i] = (ss[:, 1] ./ sqrt.(ss[1, 1] .* diag(ss)))[2:end]
end
end
d = zeros(n, k);
for i = 1:n
for j = 1:k
d[i, j] = sum((canscore[i, :] .- centroids[j, :]) .^ 2)
end
end
pvalue2 = pchisq.(d, p, false)
pBayes = exp.(-d ./ 2) .* ni';
pBayes = pBayes ./ sum(pBayes, dims=2)
classification = [gname[argmax(pBayes[i, :])] for i = 1:n]
real, pred, correcttable = table(group, classification)
correctratio = 100sum(diag(correcttable)) / n
obj = Dict(:means => means, :b => b, :w => w, :withincov => withincov,
:withinr => withinr, :eigenvalues => eigenvalues,
:wilksFromF => wilksFromF, :Fvalue => Fvalue, :df1 => df1,
:df2 => df2, :pvalue => pvalue, :Λ => Λ, :chisq => chisq,
:df => df, :pwilks => pwilks,
:canonicalcorrcoef => canonicalcorrcoef,
:stdcoefficient => stdcoefficient, :structure => structure,
:coefficient => coefficient, :centroids => centroids,
:canscore => canscore, :pBayes => pBayes, :pvalue2 => pvalue2,
:classification => :classification, :correcttable => correcttable,
:real => real, :pred => pred, :correctratio => correctratio,
:group => group, :ngroup => k, :nax => nax, :vnames => vnames,
:gname => gname, :axnames => axnames)
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
function onewaytest(ni, mi, vi)
n = sum(ni)
k = length(ni)
gm = sum(ni .* mi) / n
df1, df2 = k - 1, n - k
return (sum(ni .* (mi .- gm).^2) / df1) / (sum((ni .- 1) .* vi) / df2)
end
function printcandis(obj)
vnames = obj[:vnames]
gname = obj[:gname]
axnames = obj[:axnames]
p = length(vnames)
n = size(obj[:canscore], 1)
println("\n全体および各群の平均値\n", NamedArray(obj[:means],
(vnames, vcat("grand mean", gname))))
println("\n群間平方和\n", NamedArray(obj[:b], (vnames, vnames)))
println("\n群内平方和\n", NamedArray(obj[:w], (vnames, vnames)))
println("\n平均値の差の検定（一変量;一元配置分散分析）\n",
NamedArray(hcat(obj[:wilksFromF], obj[:Fvalue], repeat([obj[:df1]], p), repeat([obj[:df2]], p), obj[:pvalue]),
(vnames, ["Wilks", "F value", "df1", "df2", "p value"])))
println("\nプールした分散・共分散行列\n", NamedArray(obj[:withincov], (vnames, vnames)))
println("\nプールした相関係数行列\n", NamedArray(obj[:withinr], (vnames, vnames)))
println("\n固有値 = $(obj[:eigenvalues])")
println("\nWilks のΛ = $(obj[:Λ])")
println("\nカイ二乗値（近似値） = $(obj[:chisq])")
println("\n自由度 = $(obj[:df])")
println("\nP 値 = $(obj[:pwilks])")
println("\n正準相関係数 = $(obj[:canonicalcorrcoef])")
println("\n標準化判別係数\n", NamedArray(obj[:stdcoefficient],
(vnames, axnames)))
println("\n構造行列\n", NamedArray(obj[:structure], (vnames, axnames)))
println("\n判別係数\n", NamedArray(obj[:coefficient],
(vcat(vnames, "constant"), axnames)))
println("\n各群の重心\n", NamedArray(obj[:centroids], (gname, axnames)))
println("\n各ケースの判別スコア\n", NamedArray(obj[:canscore], (1:n, axnames)))
println("\n各ケースが各群に所属するベイズ確率\n", NamedArray(obj[:pBayes], (1:n, gname)))
println("\n各ケースがそれぞれの群に属するとしたとき，その判別値をとる確率\n", NamedArray(obj[:pvalue2], (1:n, gname)))
println("\n判別結果\n", NamedArray(obj[:correcttable],
(obj[:real], obj[:pred])))
println("\n正判別率 = $(obj[:correctratio])\n")
end
function plotcandis(obj; axis = 1:2, which = "boxplot", # or "barplot"
nclass = 20)
score = obj[:canscore];
group = obj[:group];
gname = obj[:gname];
k = length(gname)
pyplot()
if size(score, 2) >= 2
plt = plot(grid=false, tick_direction=:out, label="")
for g in gname # g = "virginica"
suf = group .== g
scatter!(score[suf, axis[1]], score[suf, axis[2]], label=g)
end
else
if which == "boxplot"
plt = boxplot(string.(group), obj[:canscore], xlabel = "群", ylabel = "判別値", label="")
else
canscore = obj[:canscore];
minx, maxx = extrema(canscore)
w = (maxx - minx) / (nclass - 1)
canscore2 = floor.(Int, (canscore .- minx) ./ w)
index1, index2, res = table(canscore2, group)
plt = groupedbar(res, xlabel = "判別値($nclass 階級に基準化)", label="")
end
end
display(plt)
end
using RDatasets
iris = dataset("datasets", "iris");
data = iris[:, 1:4];
group = iris[:, 5];
obj = candis(data, group)
printcandis(obj)
#=====
全体および各群の平均値
4×4 Named Matrix{Float64}
A ╲ B │ grand mean setosa versicolor virginica
────────────┼───────────────────────────────────────────────
SepalLength │ 5.84333 5.006 5.936 6.588
SepalWidth │ 3.05733 3.428 2.77 2.974
PetalLength │ 3.758 1.462 4.26 5.552
PetalWidth │ 1.19933 0.246 1.326 2.026
群間平方和
4×4 Named Matrix{Float64}
A ╲ B │ SepalLength SepalWidth PetalLength PetalWidth
────────────┼───────────────────────────────────────────────────
SepalLength │ 63.2121 -19.9527 165.248 71.2793
SepalWidth │ -19.9527 11.3449 -57.2396 -22.9327
PetalLength │ 165.248 -57.2396 437.103 186.774
PetalWidth │ 71.2793 -22.9327 186.774 80.4133
群内平方和
4×4 Named Matrix{Float64}
A ╲ B │ SepalLength SepalWidth PetalLength PetalWidth
────────────┼───────────────────────────────────────────────────
SepalLength │ 38.9562 13.63 24.6246 5.645
SepalWidth │ 13.63 16.962 8.1208 4.8084
PetalLength │ 24.6246 8.1208 27.2226 6.2718
PetalWidth │ 5.645 4.8084 6.2718 6.1566
平均値の差の検定（一変量;一元配置分散分析）
4×5 Named Matrix{Float64}
A ╲ B │ Wilks F value df1 df2 p value
────────────┼────────────────────────────────────────────────────────────────
SepalLength │ 0.381294 119.265 2.0 147.0 1.66967e-31
SepalWidth │ 0.599217 49.16 2.0 147.0 4.49202e-17
PetalLength │ 0.0586283 1180.16 2.0 147.0 2.85678e-91
PetalWidth │ 0.0711171 960.007 2.0 147.0 4.16945e-85
プールした分散・共分散行列
4×4 Named Matrix{Float64}
A ╲ B │ SepalLength SepalWidth PetalLength PetalWidth
────────────┼───────────────────────────────────────────────────
SepalLength │ 0.265008 0.0927211 0.167514 0.0384014
SepalWidth │ 0.0927211 0.115388 0.0552435 0.0327102
PetalLength │ 0.167514 0.0552435 0.185188 0.0426653
PetalWidth │ 0.0384014 0.0327102 0.0426653 0.0418816
プールした相関係数行列
4×4 Named Matrix{Float64}
A ╲ B │ SepalLength SepalWidth PetalLength PetalWidth
────────────┼───────────────────────────────────────────────────
SepalLength │ 1.0 0.530236 0.756164 0.364506
SepalWidth │ 0.530236 1.0 0.377916 0.470535
PetalLength │ 0.756164 0.377916 1.0 0.484459
PetalWidth │ 0.364506 0.470535 0.484459 1.0
固有値 = [32.19192919827802, 0.28539104262307013]
Wilks のΛ = [0.023438630650878412, 0.7779733690685453]
カイ二乗値（近似値） = [546.1152964877398, 36.52966437258939]
自由度 = [8, 3]
P 値 = [8.870784815903418e-113, 5.786050138401111e-8]
正準相関係数 = [0.9848208944320842, 0.471197019230231]
標準化判別係数
4×2 Named Matrix{Float64}
A ╲ B │ axis 1 axis 2
────────────┼─────────────────────
SepalLength │ -0.426955 0.0124075
SepalWidth │ -0.521242 0.735261
PetalLength │ 0.947257 -0.401038
PetalWidth │ 0.575161 0.58104
構造行列
4×2 Named Matrix{Float64}
A ╲ B │ axis 1 axis 2
────────────┼─────────────────────
SepalLength │ 0.222596 0.310812
SepalWidth │ -0.119012 0.863681
PetalLength │ 0.706065 0.167701
PetalWidth │ 0.633178 0.737242
判別係数
5×2 Named Matrix{Float64}
A ╲ B │ axis 1 axis 2
────────────┼─────────────────────
SepalLength │ -0.829378 0.0241021
SepalWidth │ -1.53447 2.16452
PetalLength │ 2.20121 -0.931921
PetalWidth │ 2.81046 2.83919
constant │ -2.10511 -6.66147
各群の重心
3×2 Named Matrix{Float64}
A ╲ B │ axis 1 axis 2
───────────┼───────────────────
setosa │ -7.6076 0.215133
versicolor │ 1.82505 -0.7279
virginica │ 5.78255 0.512767
各ケースの判別スコア
150×2 Named Matrix{Float64}
A ╲ B │ axis 1 axis 2
──────┼───────────────────────
1 │ -8.0618 0.300421
2 │ -7.12869 -0.78666
3 │ -7.48983 -0.265384
4 │ -6.8132 -0.670631
5 │ -8.13231 0.514463
6 │ -7.70195 1.46172
7 │ -7.21262 0.355836
8 │ -7.60529 -0.0116338
9 │ -6.56055 -1.01516
10 │ -7.34306 -0.947319
11 │ -8.39739 0.647363
12 │ -7.2193 -0.109646
⋮ ⋮ ⋮
139 │ 3.93985 0.61402
140 │ 5.20383 1.14477
141 │ 6.65309 1.80532
142 │ 5.10556 1.99218
143 │ 5.50748 -0.035814
144 │ 6.79602 1.46069
145 │ 6.84736 2.42895
146 │ 5.645 1.67772
147 │ 5.17956 -0.363475
148 │ 4.96774 0.821141
149 │ 5.88615 2.34509
150 │ 4.68315 0.332034
各ケースが各群に所属するベイズ確率
150×3 Named Matrix{Float64}
A ╲ B │ setosa versicolor virginica
──────┼──────────────────────────────────────
1 │ 1.0 3.89636e-22 2.61117e-42
2 │ 1.0 7.21797e-18 5.04214e-37
3 │ 1.0 1.46385e-19 4.67593e-39
4 │ 1.0 1.26854e-16 3.56661e-35
5 │ 1.0 1.63739e-22 1.08261e-42
6 │ 1.0 3.88328e-21 4.56654e-40
7 │ 1.0 1.11347e-18 2.30261e-37
8 │ 1.0 3.87759e-20 1.0745e-39
9 │ 1.0 1.90281e-15 9.48294e-34
10 │ 1.0 1.1118e-18 2.72406e-38
11 │ 1.0 1.18528e-23 3.23708e-44
12 │ 1.0 1.62165e-18 1.8332e-37
⋮ ⋮ ⋮ ⋮
139 │ 4.53863e-29 0.192526 0.807474
140 │ 2.14023e-36 0.00082909 0.999171
141 │ 6.5709e-45 1.18081e-6 0.999999
142 │ 6.20259e-36 0.00042764 0.999572
143 │ 5.2138e-38 0.00107825 0.998922
144 │ 1.07395e-45 1.02852e-6 0.999999
145 │ 4.04825e-46 2.52498e-7 1.0
146 │ 4.97007e-39 7.47336e-5 0.999925
147 │ 4.61661e-36 0.00589878 0.994101
148 │ 5.54896e-35 0.00314587 0.996854
149 │ 1.61369e-40 1.25747e-5 0.999987
150 │ 2.85801e-33 0.0175423 0.982458
各ケースがそれぞれの群に属するとしたとき，その判別値をとる確率
150×3 Named Matrix{Float64}
A ╲ B │ setosa versicolor virginica
──────┼──────────────────────────────────────
1 │ 0.994689 1.765e-20 2.27291e-40
2 │ 0.872645 1.60097e-16 2.31897e-35
3 │ 0.993095 5.76252e-18 3.69796e-37
4 │ 0.841471 2.39456e-15 1.4239e-33
5 │ 0.985247 7.00518e-21 8.82307e-41
6 │ 0.815447 8.67075e-20 1.93074e-38
7 │ 0.996356 4.32657e-17 1.80199e-35
8 │ 0.999675 1.72791e-18 9.504e-38
9 │ 0.625063 1.86796e-14 2.01467e-32
10 │ 0.840489 2.35172e-17 1.18057e-36
11 │ 0.937023 4.28313e-22 2.19176e-42
12 │ 0.99246 6.00485e-17 1.3824e-35
⋮ ⋮ ⋮ ⋮
139 │ 6.93678e-28 0.179657 0.492348
140 │ 1.23891e-34 0.00486351 0.947027
141 │ 2.02812e-43 5.56203e-6 0.657478
142 │ 1.3775e-34 0.00114802 0.618523
143 │ 3.763e-36 0.00717075 0.984347
144 │ 4.32626e-44 6.18508e-6 0.749428
145 │ 3.95243e-45 4.24745e-7 0.307834
146 │ 2.24534e-37 0.0004203 0.84835
147 │ 2.18721e-34 0.0225559 0.88926
148 │ 3.05651e-33 0.0154124 0.943857
149 │ 2.82485e-39 3.26024e-5 0.498209
150 │ 1.19695e-31 0.0541962 0.871247
判別結果
3×3 Named Matrix{Int64}
A ╲ B │ setosa versicolor virginica
───────────┼───────────────────────────────────
setosa │ 50 0 0
versicolor │ 0 48 2
virginica │ 0 1 49
正判別率 = 98.0
=====#
plotcandis(obj) # savefig("fig1.png")
obj2 = candis(iris[51:150, 1:4], iris[51:150, 5])
printcandis(obj2)
plotcandis(obj2, which="barplot") # savefig("fig2.png")
plotcandis(obj2, which="boxplot") # savefig("fig3.png")
#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。
判別分析（二次の判別関数）
http://aoki2.si.gunma-u.ac.jp/R/quad_disc.html
ファイル名: quaddisc.jl 関数名: quaddisc
翻訳するときに書いたメモ
==========#
using Statistics, Rmath, LinearAlgebra, NamedArrays, Plots
function quaddisc(data, group)
vname = names(data)
data = Matrix(data)
n, p = size(data)
gname, ni = table(group)
k = length(ni)
groupmeans = zeros(p, k)
vars = zeros(p, p, k)
invvars = similar(vars)
for (i, g) in enumerate(gname)
data2 = data[group .== g, :]
groupmeans[:, i] = mean(data2, dims=1)
cov2 = cov(data2)
vars[:, :, i] = cov2
invvars[:, :, i] = inv(cov2)
end
scores = zeros(n, k)
for i = 1:n
g = group[i]
temp = data[i, :]
for j = 1:k
temp2 = temp .- groupmeans[:, j]
scores[i, j] = temp2' * invvars[:, :, j] * temp2
end
end
println("\nユークリッドの二乗距離\n", NamedArray(scores, (1:n, gname)))
pvalues = pchisq.(scores, p, false)
prediction = [gname[argmax(pvalues[i, :])] for i = 1:n]
println("\n各群への所属確率\n", NamedArray(cat(pvalues, prediction, dims=2),
(1:n, vcat(gname, "prediction"))))
real, pred, correcttable = table(group, prediction)
correctrate = sum(diag(correcttable)) / n * 100
println("\n予測結果\n", NamedArray(correcttable, (real, pred)))
println("\n正判別率 = $correctrate %")
Dict(:groupmeans => groupmeans, :vars => vars, :invvars => invvars,
:scores => scores, :pvalues => pvalues, :prediction => prediction,
:correcttable => correcttable, :correctrate => correctrate)
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
using RDatasets
iris = dataset("datasets", "iris")
quaddisc(iris[:, 1:4], iris[:, 5])
#=====
ユークリッドの二乗距離
150×3 Named Matrix{Float64}
A ╲ B │ setosa versicolor virginica
──────┼───────────────────────────────────
1 │ 0.449114 114.804 182.936
2 │ 2.08109 83.3154 153.975
3 │ 1.28434 94.9204 160.494
4 │ 1.70621 82.7788 140.641
5 │ 0.761685 120.481 184.037
6 │ 3.71265 120.48 183.298
7 │ 3.4242 96.0035 154.019
8 │ 0.343439 103.815 167.444
9 │ 2.99648 73.947 132.766
10 │ 3.20009 93.4959 156.739
11 │ 1.89095 129.675 198.805
12 │ 2.01488 99.6173 154.294
⋮ ⋮ ⋮ ⋮
139 │ 489.209 8.6334 3.06747
140 │ 688.018 21.4058 2.31141
141 │ 808.103 45.4817 2.17213
142 │ 678.615 48.093 8.31315
143 │ 582.651 19.2259 1.89443
144 │ 848.939 31.0498 1.33526
145 │ 851.482 49.9147 3.1402
146 │ 698.916 45.633 4.54549
147 │ 578.018 23.3641 4.0077
148 │ 618.186 16.741 1.11181
149 │ 720.692 33.2029 3.94189
150 │ 550.788 10.1125 2.69108
各群への所属確率
150×4 Named Matrix{Any}
A ╲ B │ setosa versicolor virginica prediction
──────┼───────────────────────────────────────────────────────
1 │ 0.978262 6.86992e-24 1.74568e-38 "setosa"
2 │ 0.720846 3.45379e-17 2.86279e-32 "setosa"
3 │ 0.864028 1.18489e-19 1.14537e-33 "setosa"
4 │ 0.78959 4.48817e-17 2.05743e-29 "setosa"
5 │ 0.94351 4.21616e-25 1.01264e-38 "setosa"
6 │ 0.446289 4.21814e-25 1.45938e-38 "setosa"
7 │ 0.489498 6.97144e-20 2.80168e-32 "setosa"
8 │ 0.98684 1.51497e-21 3.698e-35 "setosa"
9 │ 0.558415 3.32733e-15 9.97411e-28 "setosa"
10 │ 0.524917 2.38001e-19 7.31636e-33 "setosa"
11 │ 0.755807 4.56942e-27 6.78869e-42 "setosa"
12 │ 0.733022 1.18664e-20 2.44609e-32 "setosa"
⋮ ⋮ ⋮ ⋮ ⋮
139 │ 1.44521e-104 0.0709452 0.546599 "virginica"
140 │ 1.36991e-147 0.000263072 0.678692 "virginica"
141 │ 1.34966e-173 3.15701e-9 0.704135 "virginica"
142 │ 1.48777e-145 9.02582e-10 0.0807578 "virginica"
143 │ 8.80541e-125 0.000709559 0.755167 "virginica"
144 │ 1.92311e-182 2.99056e-6 0.855365 "virginica"
145 │ 5.41049e-183 3.76203e-10 0.534644 "virginica"
146 │ 5.9837e-150 2.93632e-9 0.337188 "virginica"
147 │ 8.85759e-124 0.000107086 0.404965 "virginica"
148 │ 1.79562e-132 0.00217025 0.892394 "virginica"
149 │ 1.15234e-154 1.08551e-6 0.413927 "virginica"
150 │ 6.90934e-118 0.0385747 0.610776 "virginica"
予測結果
3×3 Named Matrix{Int64}
A ╲ B │ setosa versicolor virginica
───────────┼───────────────────────────────────
setosa │ 50 0 0
versicolor │ 0 47 3
virginica │ 0 0 50
正判別率 = 98.0 %
Dict{Symbol, Any} with 8 entries:
:pvalues => [0.978262 6.86992e-24 1.74568e-38; 0.720846 3.45379e-17 2.86…
:vars => [0.124249 0.0992163 0.0163551 0.0103306; 0.0992163 0.14369 0…
:scores => [0.449114 114.804 182.936; 2.08109 83.3154 153.975; … ; 720.…
:prediction => ["setosa", "setosa", "setosa", "setosa", "setosa", "setosa",…
:correcttable => [50 0 0; 0 47 3; 0 0 50]
:correctrate => 98.0
:invvars => [18.9434 -12.4048 -4.50021 -4.77613; -12.4048 15.5705 1.1110…
:groupmeans => [5.006 5.936 6.588; 3.428 2.77 2.974; 1.462 4.26 5.552; 0.24…
=====#
#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。
線形判別関数における説明変数の有意性検定
http://aoki2.si.gunma-u.ac.jp/R/PartialF.html
ファイル名: partialf.jl 関数名: partialf
翻訳するときに書いたメモ
多変量解析では結果出力に変数名が必須なので，データはデータフレームで与えるのがよい。
ただ，その後の計算ではデータフレームのままでは不便なので Matrix に変換する。
==========#
using CSV, DataFrames, Statistics, StatsBase, LinearAlgebra, Rmath, Printf
function partialf(data::DataFrame, group)
variables = names(data)
data = Matrix(data)
ncase, p = size(data)
gname, num = table(group)
ng = length(gname)
t = cov(data, corrected=false) .* ncase
w = zeros(p, p)
for (n, g) in zip(num, gname)
w .+= cov(data[group .== g, :], corrected=false) .* n
end
f = diag(inv(t)) ./ diag(inv(w))
df1 = ng - 1
df2 = ncase - (ng - 1) - p
f = df2 ./ df1 .* (1 .- f) ./ f
pvalue = pf.(f, df1, df2, false)
@printf("%15s %12s %12s\n", "", "Partial F", "pvalue")
for i = 1:p
@printf("%15s %12.6f %12.6f\n", variables[i], f[i], pvalue[i])
end
Dict(:variables => variables, :f => f, :pvalue => pvalue,
:df1 => df1, :df2 => df2)
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
using RDatasets
iris = dataset("datasets", "iris");
group = iris[:, 5];
data = iris[:, 1:4];
partialf(data, group)
# Partial F pvalue
# SepalLength 4.721152 0.010329
# SepalWidth 21.935928 0.000000
# PetalLength 35.590175 0.000000
# PetalWidth 24.904333 0.000000
#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。
判別分析（線形判別関数；ステップワイズ変数選択）
http://aoki2.si.gunma-u.ac.jp/R/sdis.html
ファイル名: sdis.jl 関数名: sdis
翻訳するときに書いたメモ
==========#
using CSV, DataFrames, Rmath, Statistics, LinearAlgebra, NamedArrays, Printf, Plots
function sdis(data::Array{Float64,2}, group; name=[], stepwise=true, Pin=0.05, Pout=0.05, predict=false, verbose=false)
getitem(t, lxi) = [t[i, i] for i in lxi]
formatpval(p) = p >= 0.000001 ? @sprintf("%.6f", p) : "< 0.000001"
function stepout(isw)
step += 1
ncasek = ncase - ng
isw != 0 && verbose && println("\n ***** ステップ $step ***** \n",
" $(["編入", "除去"][isw])変数: $(name[ip])")
lxi = lx[1:ni]
a = zeros(ni, ng)
a0 = zeros(ng)
for g = 1:ng
a[:, g] = -(w[lxi, lxi] * means[lxi, g]) * 2ncasek
a0[g] = transpose(means[lxi, g]) * w[lxi, lxi] * means[lxi, g] * ncasek
end
idf1 = ng - 1
idf2 = ncase - (ng - 1) - ni
temp = idf2 / idf1
f = getitem(t, lxi) ./ getitem(w, lxi)
f = temp * (1 .- f) ./ f
P = pf.(f, idf1, idf2, false)
alp = ng - 1
b = ncase - 1 - 0.5 * (ni + ng)
qa = ni ^ 2 + alp ^ 2
c = 1
qa != 5 && (c = sqrt((ni ^ 2 * alp ^ 2 - 4) / (qa - 5)))
wldf1 = ni * alp
wldf2 = b * c + 1 - 0.5 * ni * alp
wl = detw / dett
cl = exp(log(wl) / c)
wlf = wldf2 * (1 - cl) / (wldf1 * cl)
wlp = pf(wlf, wldf1, wldf2, false)
results = merge(results, Dict(:rownames => name[lxi], :a => a, :a0 => a0,
:f => f, :idf1 => idf1, :idf2 => idf2, :P => P,
:wl => wl, :wlf => wlf, :wldf1 => wldf1,
:wldf2 => wldf2, :wlp => wlp))
end
function printclassfunc()
println("\n***** 分類関数 *****")
@printf("%12s", "")
for i = 1:ng
@printf(" %12s", gname[i])
end
@printf(" %12s %12s\n", "偏F値", "P値")
rownames = results[:rownames]
for j = 1:length(rownames)
@printf("%12s", rownames[j])
for i = 1:ng
@printf(" %12.6f", results[:a][j, i])
end
@printf(" %12.6f %12s\n", results[:f][j], formatpval(results[:P][j]))
end
@printf("%12s", "定数項")
for i = 1:ng
@printf(" %12.6f", results[:a0][i])
end
println("\n\nウィルクスのΛ: $(results[:wl])")
println("等価なF値: $(results[:wlf])")
println("自由度: ($(results[:wldf1]), $(results[:wldf2]))")
println("P値: $(results[:wlp])")
end
function fmax()
kouho = 1:p
if ni > 0
suf = trues(p)
[suf[i] = false for i in lx[1:ni]]
kouho = (1:p)[suf]
end
temp = getitem(w, kouho) ./ getitem(t, kouho)
temp = (1 .- temp) ./ temp
ip = argmax(temp)
return temp[ip], kouho[ip]
end
function fmin()
kouho = lx[1:ni]
temp = getitem(t, kouho) ./ getitem(w, kouho)
temp = (1 .- temp) ./ temp
ip = argmin(temp)
return temp[ip], lx[ip]
end
function sweepsdis!(r, det, ip, p)
ap = r[ip, ip]
abs(ap) > EPSINV || error("正規方程式の係数行列が特異行列です")
det *= ap
for i = 1:p
if i != ip
temp = r[ip, i] / ap
for j = 1:p
if j != ip
r[j, i] -= r[j, ip] * temp
end
end
end
end
r[:, ip] /= ap
r[ip, :] /= -ap
r[ip, ip] = 1 / ap
det
end
function discriminantfunction()
lxi = lx[1:ni]
side = name[lxi]
ncasek = ncase - ng
p0 = length(lxi)
m = ng * (ng - 1) ÷ 2
dfunc = zeros(p0 + 1, m)
stddfunc = zeros(p0, m)
header = fill("", m)
dist = zeros(m)
errorp = zeros(m)
k = 0
for g1 in 1:ng - 1
for g2 in g1 + 1:ng
k += 1
header[k] = gname[g1] * ":" * gname[g2]
xx = means[lxi, g1] .- means[lxi, g2]
fn = w[lxi, lxi] * xx .* ncasek
stddfunc[1:p0, k] = sd[lxi] .* fn
fn0 = -sum(fn .* (means[lxi, g1] .+ means[lxi, g2]) .* 0.5)
dfunc[1:p0, k] = fn
dfunc[p0 + 1, k] = fn0
dist[k] = sqrt(sum(xx .* fn))
errorp[k] = pnorm(dist[k] * 0.5, false)
end
end
printdfunc(header, side, dfunc, stddfunc, dist, errorp)
results = merge(results, Dict(:header => header, :side => side, :dfunc => dfunc,
:dist => dist, :errorp => errorp))
end
function printdfunc(header, side, dfunc, stddfunc, dist, errorp)
simpleheader = "Func." .* string.(1:length(header))
array1 = NamedArray(dfunc,
(vcat(side, "Constant"), simpleheader))
println("\n判別関数\n", array1)
array2 = NamedArray(stddfunc, (side, simpleheader))
println("\n標準化判別関数\n", array2)
array3 = NamedArray(vcat(dist', errorp'),
(["マハラノビスの汎距離", "理論的誤判別率"], simpleheader))
println("\n", array3)
for i = 1:length(header)
println(" Func.$i: $(header[i])")
end
end
function procpredict()
nc0 = 0
ncasek = ncase - ng
lxi = lx[1:ni]
data = data[:, lxi]
tdata = transpose(data)
dis = zeros(ncase, ng)
for g = 1:ng
temp = tdata .- means[lxi, g]
for j = 1:ncase
dis[j, g] = temp[:, j]' * w[lxi, lxi] * temp[:, j]
end
end
dis .*= ncase - ng
P = pchisq.(dis, p, false)
prediction = [gname[argmax(P[j, :])] for j = 1:ncase]
index1, tbl = table(prediction)
correct = prediction .== group
factor1, factor2, correcttable = table(group, prediction)
correctrate = sum(diag(correcttable)) / ncase * 100
if ng == 2
fn = results[:dfunc][1:end-1]
fn0 = results[:dfunc][end]
dfv = data * fn .+ fn0
else
dfv = NaN
end
results = merge(results, Dict(:dis => dis, :P => P,
:prediction => prediction, :correct => correct,
:correcttable => correcttable,
:factor1 => factor1, :factor2 => factor2,
:correctrate => correctrate, :dfv => dfv))
end
EPSINV = 1e-6
Pout >= Pin || (Pout = Pin)
step = 0
ip = 0
lxi = []
ncase, p = size(data)
p > 1 || (stepwise = false)
length(name) != 0 || (name = ["x" * string(i) for i = 1:p])
gname, num = table(group)
ng = length(gname)
ng > 1 || error("1群しかありません")
any(num .>= 2) || error("ケース数が1以下の群があります")
gmeans = vec(mean(data, dims=1))
t = cov(data, corrected=false) .* ncase
means = zeros(p, ng)
vars = zeros(p, p, ng)
for i = 1:ng
gdata = data[group .== gname[i], :]
means[:, i] = vec(mean(gdata, dims=1))
vars[:, :, i] = cov(gdata, corrected=false) .* num[i]
end
if verbose
println("有効ケース数： $ncase")
println("判別するグループ： $gname")
println("***** 平均値 *****")
println(NamedArray(hcat(means, gmeans),
(name, vcat(gname, "全体"))))
end
w = sum(vars, dims=3)[:, :, 1]
detw = dett = 1
sd2 = sqrt.(diag(w) ./ ncase)
r = w ./ (sd2 * sd2') ./ ncase
if verbose
println("\n***** プールされた群内相関係数行列 *****")
println(NamedArray(r, (name, name)))
end
sd = sqrt.(diag(t) ./ ncase)
results = Dict(:ncase => ncase, :gname => gname,
:means => means, :gmeans => gmeans, :r => r)
if stepwise == false
for ip = 1:p
detw = sweepsdis!(w, detw, ip, p)
dett = sweepsdis!(t, dett, ip, p)
end
lx = 1:p
ni = p
stepout(0)
else
verbose && println("\n変数編入基準 Pin: $Pin\n",
"変数除去基準 Pout: $Pout")
lx = zeros(Int, p)
ni = 0
while ni != p
ansmax, ip = fmax()
P = (ncase - ng - ni) / (ng - 1) * ansmax
P = pf(P, ng - 1, ncase - ng - ni, false)
verbose && println("編入候補変数: $(name[ip]) P: $(formatpval(P))")
if P > Pin
verbose && println(" 編入されませんでした")
break
end
verbose && println(" 編入されました")
ni += 1
lx[ni] = ip
detw = sweepsdis!(w, detw, ip, p)
dett = sweepsdis!(t, dett, ip, p)
stepout(1)
verbose && printclassfunc()
while true
ansmin, ip = fmin()
P = (ncase - ng - ni + 1) / (ng - 1) * ansmin
P = pf(P, ng - 1, ncase - ng - ni + 1, false)
verbose && println("\n除去候補変数: $(name[ip]) P: $(formatpval(P))")
if P <= Pout
verbose && println(" 除去されませんでした")
break
else
verbose && println(" 除去されました")
lx = lx[lx .!= ip]
ni -= 1
w, detw = sweepsdis(w, detw)
t, dett = sweepsdis(t, dett)
stepout(2)
verbose && printclassfunc()
end
end
end
end
ni == 0 && error("条件（Pin < $Pin）を満たす独立変数がありません")
verbose && println("\n========== 結果 ==========")
printclassfunc()
discriminantfunction()
procpredict()
if predict
println("\n********** 各ケースの判別結果 **********")
println("\n各群への二乗距離")
if ng == 2
println(NamedArray(hcat(results[:dis], results[:dfv]), (1:ncase, vcat(gname, "判別値"))))
else
println(NamedArray(results[:dis], (1:ncase, gname)))
end
println("\nP 値")
println(NamedArray(results[:P], (1:ncase, gname)))
println(" メモ:「二乗距離」とは，各群の重心までのマハラノビスの汎距離の二乗です。")
println(" P値は各群に属する確率です。")
println("\n判別，判別の正否")
println(NamedArray(hcat(string.(group), results[:prediction], results[:correct]),
(1:ncase, ["実際の群", "判別された群", "正否"])))
end
println("\n判別結果")
println(NamedArray(results[:correcttable], (results[:factor1], results[:factor2])))
println("\n正判別率 = $(results[:correctrate])\n")
results
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
function sdis(data::DataFrame, group; stepwise=true, Pin=0.05, Pout=0.05, predict=false, verbose=false)
sdis(Matrix(data), group, name=names(data), stepwise, Pin, Pout, predict, verbose)
end
function sdis(data::Array{Int64,2}, group; name=[], stepwise=true, Pin=0.05, Pout=0.05, predict=false, verbose=false)
sdis(Matrix(data), group, name, stepwise, Pin, Pout, predict, verbose)
end
using RDatasets
iris = dataset("datasets", "iris")
name = names(iris)[1:4]
data = Matrix(iris[51:150, 1:4])
group = vec(iris[51:150, 5])
sdis(data, group, name=name, verbose=false)
***** 分類関数 *****
versicolor virginica 偏F値 P値
PetalWidth 1.172895 -23.599188 37.091608 < 0.000001
SepalWidth -31.942816 -20.785575 10.587491 0.001578
PetalLength 5.163281 -8.776974 24.156597 0.000004
SepalLength -30.802050 -23.689444 7.367913 0.007886
定数項 123.885865 157.212036
ウィルクスのΛ: 0.21611029704367302
等価なF値: 86.14758620895533
自由度: (4, 95.0)
P値: 9.53987626477818e-31
判別関数
5×1 Named Matrix{Float64}
A ╲ B │ Func.1
────────────┼─────────
PetalWidth │ -12.386
SepalWidth │ 5.57862
PetalLength │ -6.97013
SepalLength │ 3.5563
Constant │ 16.6631
標準化判別関数
4×1 Named Matrix{Float64}
A ╲ B │ Func.1
────────────┼─────────
PetalWidth │ -5.23483
SepalWidth │ 1.84699
PetalLength │ -5.72554
SepalLength │ 2.34542
2×1 Named Matrix{Float64}
A ╲ B │ Func.1
───────────┼──────────
マハラノビスの汎距離 │ 3.77079
理論的誤判別率 │ 0.0296881
Func.1: versicolor:virginica
判別結果
2×2 Named Matrix{Int64}
A ╲ B │ versicolor virginica
───────────┼───────────────────────
versicolor │ 48 2
virginica │ 1 49
正判別率 = 97.0
#=====
Dict{Symbol, Any} with 30 entries:
:dis => [5.29832 23.9158; 2.105 16.7657; … ; 21.4011 4.50692; 10.154…
:factor1 => ["versicolor", "virginica"]
:means => [5.936 6.588; 2.77 2.974; 4.26 5.552; 1.326 2.026]
:ncase => 100
:wlf => 86.1476
:prediction => ["versicolor", "versicolor", "versicolor", "versicolor", "ve…
:header => ["versicolor:virginica"]
:gmeans => [6.262, 2.872, 4.906, 1.676]
:wlp => 9.53988e-31
:correct => Bool[1, 1, 1, 1, 1, 1, 1, 1, 1, 1 … 1, 1, 1, 1, 1, 1, 1, 1…
:correctrate => 97.0
:dfv => [9.30873, 7.33037, 5.76261, 5.07121, 4.75754, 5.08672, 4.899…
:wldf2 => 95.0
:rownames => ["PetalWidth", "SepalWidth", "PetalLength", "SepalLength"]
:gname => ["versicolor", "virginica"]
:factor2 => ["versicolor", "virginica"]
:errorp => [0.0296881]
:dist => [3.77079]
:idf1 => 1
:correcttable => [48 2; 1 49]
:wldf1 => 4
:r => [1.0 0.48557 0.818971 0.378356; 0.48557 1.0 0.472261 0.58332…
:idf2 => 95
:wl => 0.21611
:a => [1.17289 -23.5992; -31.9428 -20.7856; 5.16328 -8.77697; -30.…
⋮ => ⋮
====#
name = names(iris)[1:4]
data = Matrix(iris[:, 1:4])
group = vec(iris[:, 5])
sdis(data, group, name=name, verbose=false)
#=====
***** 分類関数 *****
setosa versicolor virginica 偏F値 P値
PetalLength 32.861278 -10.422902 -25.533090 35.590175 < 0.000001
SepalWidth -47.175741 -14.145020 -7.370559 21.935928 < 0.000001
PetalWidth 34.796822 -12.868458 -42.158226 24.904333 < 0.000001
SepalLength -47.088333 -31.396418 -24.891698 4.721152 0.010329
定数項 170.419715 143.507990 206.539415
ウィルクスのΛ: 0.023438630650878322
等価なF値: 199.14534354008427
自由度: (8, 288.0)
P値: 1.3650058325897325e-112
判別関数
5×3 Named Matrix{Float64}
A ╲ B │ Func.1 Func.2 Func.3
────────────┼─────────────────────────────
PetalLength │ -21.6421 -29.1972 -7.55509
SepalWidth │ 16.5154 19.9026 3.38723
PetalWidth │ -23.8326 -38.4775 -14.6449
SepalLength │ 7.84596 11.0983 3.25236
Constant │ -13.4559 18.0599 31.5157
標準化判別関数
4×3 Named Matrix{Float64}
A ╲ B │ Func.1 Func.2 Func.3
────────────┼─────────────────────────────
PetalLength │ -38.0772 -51.3696 -13.2925
SepalWidth │ 7.17445 8.6459 1.47145
PetalWidth │ -18.1055 -29.2311 -11.1256
SepalLength │ 6.47528 9.15946 2.68418
2×3 Named Matrix{Float64}
A ╲ B │ Func.1 Func.2 Func.3
───────────┼──────────────────────────────────────
マハラノビスの汎距離 │ 9.47967 13.3935 4.14742
理論的誤判別率 │ 1.06946e-6 1.06568e-11 0.0190532
Func.1: setosa:versicolor
Func.2: setosa:virginica
Func.3: versicolor:virginica
判別結果
3×3 Named Matrix{Int64}
A ╲ B │ setosa versicolor virginica
───────────┼───────────────────────────────────
setosa │ 50 0 0
versicolor │ 0 48 2
virginica │ 0 1 49
正判別率 = 98.0
Dict{Symbol, Any} with 30 entries:
:dis => [0.29109 98.8847 191.789; 2.03135 80.9713 169.187; … ; 188.8…
:factor1 => ["setosa", "versicolor", "virginica"]
:means => [5.006 5.936 6.588; 3.428 2.77 2.974; 1.462 4.26 5.552; 0.24…
:ncase => 150
:wlf => 199.145
:prediction => ["setosa", "setosa", "setosa", "setosa", "setosa", "setosa",…
:header => ["setosa:versicolor", "setosa:virginica", "versicolor:virgin…
:gmeans => [5.84333, 3.05733, 3.758, 1.19933]
:wlp => 1.36501e-112
:correct => Bool[1, 1, 1, 1, 1, 1, 1, 1, 1, 1 … 1, 1, 1, 1, 1, 1, 1, 1…
:correctrate => 98.0
:dfv => NaN
:wldf2 => 288.0
:rownames => ["PetalLength", "SepalWidth", "PetalWidth", "SepalLength"]
:gname => ["setosa", "versicolor", "virginica"]
:factor2 => ["setosa", "versicolor", "virginica"]
:errorp => [1.06946e-6, 1.06568e-11, 0.0190532]
:dist => [9.47967, 13.3935, 4.14742]
:idf1 => 2
:correcttable => [50 0 0; 0 48 2; 0 1 49]
:wldf1 => 8
:r => [1.0 0.530236 0.756164 0.364506; 0.530236 1.0 0.377916 0.470…
:idf2 => 144
:wl => 0.0234386
:a => [32.8613 -10.4229 -25.5331; -47.1757 -14.145 -7.37056; 34.79…
⋮ => ⋮
=====#
