#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。
名前 pls::svdpc.fit()
ファイル名: svdpcfit.jl 関数名: svdpcfit
翻訳するときに書いたメモ
特異値分解の関数で La.svd() が使われているが，Julia では svd() で問題ないようだ。
crossprod() も単純なことしかやっていないので，適切に行列演算をするだけだ。
変数名は変更した。
==========#
using Statistics, LinearAlgebra, NamedArrays
function svdpcfit(X, Y, ncomp; dnX = [], dnY = [], center = true, stripped = false)
nobj, npred = size(X)
nresp = size(Y, 2)
length(dnX) == 0 && (dnX = ["X$i" for i = 1:npred])
length(dnY) == 0 && (dnY = ["Y$i" for i = 1:nresp])
(nresp == 1 && dnY[1] == "Y1") && (dnY[1] = "Y")
coefficients = zeros(npred, nresp, ncomp);
if !stripped
fittedvalues = zeros(nobj, nresp, ncomp);
end
if center
Xmeans = vec(mean(X, dims=1));
X = (X' .- Xmeans)';
Ymeans = vec(mean(Y, dims=1));
Y = (Y' .- Ymeans)';
else
Xmeans = zeros(npred)
Ymeans = zeros(nresp)
end
u, d, vt = svd(X);
D = d[1:ncomp];
scores = u[:, 1:ncomp] * diagm(D);
loadings = vt[:, 1:ncomp];
Xvar = D .^ 2;
Xtotvar = sum(X .^ 2)
tQ = (scores' * Y) ./ Xvar;
Yloadings = transpose(tQ)
for a = 1:ncomp
coefficients[:, :, a] = loadings[:, 1:a] * tQ[1:a, :]
if !stripped
fittedvalues[:, :, a] = scores[:, 1:a] * tQ[1:a, :]
end
end
if stripped
Dict(:coefficients => coefficients, :Xmeans => Xmeans, :Ymeans => Ymeans)
else
residuals = Y .- fittedvalues;
fittedvalues = Ymeans' .+ fittedvalues
compnames = ["Comp $i" for i = 1:ncomp]
nCompnames = ["$i comps" for i = 1:ncomp]
Dict(:nobj => nobj, :npred => npred, :nresp => nresp,
:dnX => dnX, :dnY => dnY, :ncomp => ncomp,
:compnames => compnames, :nCompnames => nCompnames,
:coefficients => coefficients, :scores => scores, :loadings => loadings,
:Yloadings => Yloadings, :projection => loadings,
:Xmeans => Xmeans, :Ymeans => Ymeans,
:fittedvalues => fittedvalues, :residuals => residuals,
:Xvar => Xvar, :Xtotvar => Xtotvar)
end
end
function svdpcprint(obj)
ncomp = obj[:ncomp]
nobj = obj[:nobj]
dnX = obj[:dnX]
dnY = obj[:dnY]
compnames = obj[:compnames]
println("coefficients\n")
for i = 1:ncomp
println(", , $(obj[:nCompnames][i])\n",
NamedArray(obj[:coefficients][:, :, i], (dnX, dnY)))
end
println("\nscores\n", NamedArray(obj[:scores], (1:nobj, compnames)))
println("\nloadings\n", NamedArray(obj[:loadings], (dnX, compnames)))
SSloadings = vec(sum(obj[:loadings] .^ 2, dims=1))
ProportionVar = SSloadings ./ size(obj[:loadings], 1)
CumulativeVar = cumsum(ProportionVar)
println(NamedArray(hcat(SSloadings, ProportionVar, CumulativeVar)',
(["SS loadings", "Proportion Var", "Cumulative Var"], compnames)))
println("\nYloadings\n", NamedArray(obj[:Yloadings], (dnY, compnames)))
SSloadings = vec(sum(obj[:Yloadings] .^ 2, dims=1))
ProportionVar = SSloadings ./ size(obj[:Yloadings], 1)
CumulativeVar = cumsum(ProportionVar)
println(NamedArray(hcat(SSloadings, ProportionVar, CumulativeVar)',
(["SS loadings", "Proportion Var", "Cumulative Var"], compnames)))
println("\nprojection\n", NamedArray(obj[:projection], (dnX, compnames)))
println("\nXmeans\n", NamedArray(reshape(obj[:Xmeans], 1, :), (["Xmeans"], dnX)))
if obj[:nresp] == 1
println("\nYmean\n", obj[:Ymeans][1])
else
println("\nYmeans\n", NamedArray(reshape(obj[:Ymeans], 1, :), (["Ymeans"], dnY)))
end
println("\nfitted values\n")
for i = 1:ncomp
println(", , $(obj[:nCompnames][i])\n",
NamedArray(obj[:fittedvalues][:, :, i], (1:nobj, dnY)))
end
println("\nresiduals\n")
for i = 1:ncomp
println(", , $(obj[:nCompnames][i])\n",
NamedArray(obj[:residuals][:, :, i], (1:nobj, dnY)))
end
println("\nXvar\n", NamedArray(reshape(obj[:Xvar], 1, :), (["Xvar"], compnames)))
println("\nXtotvar = $(obj[:Xtotvar])")
end
X = [61.9 53.7 41.3 47.7 56.7
70.5 47.2 53.5 48.1 40.8
61.4 42.7 48.4 65.3 48.5
54.6 42.8 73.3 18.6 48.8
46.0 42.2 40.0 35.6 58.4
43.4 40.2 50.7 52.2 45.0
59.8 52.6 32.6 47.9 60.3
49.2 45.3 41.2 62.9 39.4
52.7 41.2 47.5 61.7 51.8
48.1 52.8 58.1 50.3 77.7
46.6 65.2 38.8 38.6 39.7
41.6 32.2 46.6 54.3 47.0
41.1 20.7 47.6 46.9 58.5
53.3 57.8 65.3 44.5 20.8
55.9 73.1 54.1 43.6 49.8];
Y = [43.3, 76.1, 64.2, 53.2, 62.8, 36.2, 66.9, 54.1,
40.6, 62.8, 57.3, 46.0, 60.6, 42.7, 54.6];
obj = svdpcfit(X, Y, 3);
svdpcprint(obj)
#=====
coefficients
, , 1 comps
5×1 Named Matrix{Float64}
A ╲ B │ Y
──────┼───────────
X1 │ -0.0120428
X2 │ -0.0346016
X3 │ -0.019896
X4 │ 0.0222783
X5 │ 0.0279437
, , 2 comps
5×1 Named Matrix{Float64}
A ╲ B │ Y
──────┼───────────
X1 │ -0.0241079
X2 │ -0.0530765
X3 │ 0.00703593
X4 │ -0.0196489
X5 │ 0.0524697
, , 3 comps
5×1 Named Matrix{Float64}
A ╲ B │ Y
──────┼───────────
X1 │ 0.0381823
X2 │ 0.145092
X3 │ -0.0978055
X4 │ -0.0585286
X5 │ 0.281049
scores
15×3 Named Matrix{Float64}
A ╲ B │ Comp 1 Comp 2 Comp 3
──────┼────────────────────────────────
1 │ -0.347277 4.41909 -13.1949
2 │ 9.79227 5.45379 4.10098
3 │ -7.78349 13.4306 3.5913
4 │ 18.6041 -32.0181 7.03149
5 │ -7.50786 -10.9465 -6.27208
6 │ -5.3772 0.233441 10.1319
7 │ -6.56632 6.23862 -17.3912
8 │ -5.82498 17.0669 8.08649
9 │ -11.1815 7.73628 3.13631
10 │ -9.59517 -13.01 -18.9783
11 │ 14.9809 6.6044 -7.27666
12 │ -14.1683 -0.102886 12.8477
13 │ -24.012 -14.1205 11.3212
14 │ 28.6131 5.6228 18.2251
15 │ 20.3738 3.39205 -15.3592
loadings
5×3 Named Matrix{Float64}
A ╲ B │ Comp 1 Comp 2 Comp 3
──────┼────────────────────────────────
X1 │ 0.219312 0.20188 -0.189628
X2 │ 0.630129 0.309135 -0.60328
X3 │ 0.362326 -0.450645 0.319166
X4 │ -0.405709 0.701556 0.11836
X5 │ -0.508881 -0.410387 -0.695858
3×3 Named Adjoint{Float64, Matrix{Float64}}
A ╲ B │ Comp 1 Comp 2 Comp 3
───────────────┼───────────────────────
SS loadings │ 1.0 1.0 1.0
Proportion Var │ 0.2 0.2 0.2
Cumulative Var │ 0.2 0.4 0.6
Yloadings
1×3 Named Transpose{Float64, Matrix{Float64}}
A ╲ B │ Comp 1 Comp 2 Comp 3
──────┼───────────────────────────────────
Y │ -0.054912 -0.0597631 -0.328486
3×3 Named Adjoint{Float64, Matrix{Float64}}
A ╲ B │ Comp 1 Comp 2 Comp 3
───────────────┼───────────────────────────────────
SS loadings │ 0.00301533 0.00357163 0.107903
Proportion Var │ 0.00301533 0.00357163 0.107903
Cumulative Var │ 0.00301533 0.00658696 0.11449
projection
5×3 Named Matrix{Float64}
A ╲ B │ Comp 1 Comp 2 Comp 3
──────┼────────────────────────────────
X1 │ 0.219312 0.20188 -0.189628
X2 │ 0.630129 0.309135 -0.60328
X3 │ 0.362326 -0.450645 0.319166
X4 │ -0.405709 0.701556 0.11836
X5 │ -0.508881 -0.410387 -0.695858
Xmeans
1×5 Named Matrix{Float64}
A ╲ B │ X1 X2 X3 X4 X5
───────┼────────────────────────────────────────────
Xmeans │ 52.4067 47.3133 49.2667 47.88 49.5467
Ymean
54.760000000000005
fitted values
, , 1 comps
15×1 Named Matrix{Float64}
A ╲ B │ Y
──────┼────────
1 │ 54.7791
2 │ 54.2223
3 │ 55.1874
4 │ 53.7384
5 │ 55.1723
6 │ 55.0553
7 │ 55.1206
8 │ 55.0799
9 │ 55.374
10 │ 55.2869
11 │ 53.9374
12 │ 55.538
13 │ 56.0785
14 │ 53.1888
15 │ 53.6412
, , 2 comps
15×1 Named Matrix{Float64}
A ╲ B │ Y
──────┼────────
1 │ 54.515
2 │ 53.8964
3 │ 54.3847
4 │ 55.6519
5 │ 55.8265
6 │ 55.0413
7 │ 54.7477
8 │ 54.0599
9 │ 54.9117
10 │ 56.0644
11 │ 53.5427
12 │ 55.5442
13 │ 56.9224
14 │ 52.8528
15 │ 53.4385
, , 3 comps
15×1 Named Matrix{Float64}
A ╲ B │ Y
──────┼────────
1 │ 58.8493
2 │ 52.5492
3 │ 53.2051
4 │ 53.3422
5 │ 57.8868
6 │ 51.7131
7 │ 60.4605
8 │ 51.4036
9 │ 53.8814
10 │ 62.2985
11 │ 55.9329
12 │ 51.3239
13 │ 53.2036
14 │ 46.8661
15 │ 58.4838
residuals
, , 1 comps
15×1 Named Matrix{Float64}
A ╲ B │ Y
──────┼──────────
1 │ -11.4791
2 │ 21.8777
3 │ 9.01259
4 │ -0.538414
5 │ 7.62773
6 │ -18.8553
7 │ 11.7794
8 │ -0.979861
9 │ -14.774
10 │ 7.51311
11 │ 3.36263
12 │ -9.53801
13 │ 4.52145
14 │ -10.4888
15 │ 0.958767
, , 2 comps
15×1 Named Matrix{Float64}
A ╲ B │ Y
──────┼─────────
1 │ -11.215
2 │ 22.2036
3 │ 9.81525
4 │ -2.45191
5 │ 6.97353
6 │ -18.8413
7 │ 12.1523
8 │ 0.040109
9 │ -14.3117
10 │ 6.73559
11 │ 3.75733
12 │ -9.54416
13 │ 3.67757
14 │ -10.1528
15 │ 1.16149
, , 3 comps
15×1 Named Matrix{Float64}
A ╲ B │ Y
──────┼──────────
1 │ -15.5493
2 │ 23.5508
3 │ 10.9949
4 │ -0.142171
5 │ 4.91324
6 │ -15.5131
7 │ 6.43952
8 │ 2.6964
9 │ -13.2814
10 │ 0.501486
11 │ 1.36705
12 │ -5.32387
13 │ 7.39641
14 │ -4.1661
15 │ -3.8838
Xvar
1×3 Named Matrix{Float64}
A ╲ B │ Comp 1 Comp 2 Comp 3
──────┼──────────────────────────
Xvar │ 3117.67 2220.14 2047.35
Xtotvar = 9027.521333333332
=====#
Y2 = [ 46.0 58.3
53.0 47.7
58.2 66.4
44.8 72.4
45.2 41.8
32.2 48.0
50.8 53.9
65.9 56.0
61.6 49.7
41.8 29.5
55.2 48.1
51.4 60.7
60.5 45.1
60.3 50.3
47.5 43.0 ];
obj2 = svdpcfit(X, Y2, 3);
svdpcprint(obj2)
#=====
coefficients
, , 1 comps
5×2 Named Matrix{Float64}
A ╲ B │ Y1 Y2
──────┼─────────────────────────
X1 │ -0.00683921 0.01963
X2 │ -0.0196505 0.0564012
X3 │ -0.0112991 0.0324309
X4 │ 0.012652 -0.036314
X5 │ 0.0158694 -0.0455487
, , 2 comps
5×2 Named Matrix{Float64}
A ╲ B │ Y1 Y2
──────┼───────────────────────
X1 │ 0.0597084 0.0233661
X2 │ 0.0822523 0.0621222
X3 │ -0.159849 0.0240911
X4 │ 0.243912 -0.0233308
X5 │ -0.11941 -0.0531434
, , 3 comps
5×2 Named Matrix{Float64}
A ╲ B │ Y1 Y2
──────┼───────────────────────
X1 │ 0.00973666 -0.042046
X2 │ -0.0767267 -0.145979
X3 │ -0.075741 0.134187
X4 │ 0.275103 0.0174976
X5 │ -0.302786 -0.293179
scores
15×3 Named Matrix{Float64}
A ╲ B │ Comp 1 Comp 2 Comp 3
──────┼────────────────────────────────
1 │ -0.347277 4.41909 -13.1949
2 │ 9.79227 5.45379 4.10098
3 │ -7.78349 13.4306 3.5913
4 │ 18.6041 -32.0181 7.03149
5 │ -7.50786 -10.9465 -6.27208
6 │ -5.3772 0.233441 10.1319
7 │ -6.56632 6.23862 -17.3912
8 │ -5.82498 17.0669 8.08649
9 │ -11.1815 7.73628 3.13631
10 │ -9.59517 -13.01 -18.9783
11 │ 14.9809 6.6044 -7.27666
12 │ -14.1683 -0.102886 12.8477
13 │ -24.012 -14.1205 11.3212
14 │ 28.6131 5.6228 18.2251
15 │ 20.3738 3.39205 -15.3592
loadings
5×3 Named Matrix{Float64}
A ╲ B │ Comp 1 Comp 2 Comp 3
──────┼────────────────────────────────
X1 │ 0.219312 0.20188 -0.189628
X2 │ 0.630129 0.309135 -0.60328
X3 │ 0.362326 -0.450645 0.319166
X4 │ -0.405709 0.701556 0.11836
X5 │ -0.508881 -0.410387 -0.695858
3×3 Named Adjoint{Float64, Matrix{Float64}}
A ╲ B │ Comp 1 Comp 2 Comp 3
───────────────┼───────────────────────
SS loadings │ 1.0 1.0 1.0
Proportion Var │ 0.2 0.2 0.2
Cumulative Var │ 0.2 0.4 0.6
Yloadings
2×3 Named Transpose{Float64, Matrix{Float64}}
A ╲ B │ Comp 1 Comp 2 Comp 3
──────┼───────────────────────────────────
Y1 │ -0.0311849 0.329639 0.263525
Y2 │ 0.0895075 0.0185063 0.344949
3×3 Named Adjoint{Float64, Matrix{Float64}}
A ╲ B │ Comp 1 Comp 2 Comp 3
───────────────┼───────────────────────────────────
SS loadings │ 0.00898409 0.109004 0.188435
Proportion Var │ 0.00449204 0.054502 0.0942176
Cumulative Var │ 0.00449204 0.0589941 0.153212
projection
5×3 Named Matrix{Float64}
A ╲ B │ Comp 1 Comp 2 Comp 3
──────┼────────────────────────────────
X1 │ 0.219312 0.20188 -0.189628
X2 │ 0.630129 0.309135 -0.60328
X3 │ 0.362326 -0.450645 0.319166
X4 │ -0.405709 0.701556 0.11836
X5 │ -0.508881 -0.410387 -0.695858
Xmeans
1×5 Named Matrix{Float64}
A ╲ B │ X1 X2 X3 X4 X5
───────┼────────────────────────────────────────────
Xmeans │ 52.4067 47.3133 49.2667 47.88 49.5467
Ymeans
1×2 Named Matrix{Float64}
A ╲ B │ Y1 Y2
───────┼─────────────────
Ymeans │ 51.6267 51.3933
fitted values
, , 1 comps
15×2 Named Matrix{Float64}
A ╲ B │ Y1 Y2
──────┼─────────────────
1 │ 51.6375 51.3622
2 │ 51.3213 52.2698
3 │ 51.8694 50.6967
4 │ 51.0465 53.0585
5 │ 51.8608 50.7213
6 │ 51.7944 50.912
7 │ 51.8314 50.8056
8 │ 51.8083 50.872
9 │ 51.9754 50.3925
10 │ 51.9259 50.5345
11 │ 51.1595 52.7342
12 │ 52.0685 50.1252
13 │ 52.3755 49.2441
14 │ 50.7344 53.9544
15 │ 50.9913 53.2169
, , 2 comps
15×2 Named Matrix{Float64}
A ╲ B │ Y1 Y2
──────┼─────────────────
1 │ 53.0942 51.444
2 │ 53.1191 52.3707
3 │ 56.2967 50.9452
4 │ 40.4921 52.466
5 │ 48.2524 50.5187
6 │ 51.8713 50.9164
7 │ 53.8879 50.9211
8 │ 57.4342 51.1878
9 │ 54.5255 50.5357
10 │ 47.6373 50.2937
11 │ 53.3366 52.8565
12 │ 52.0346 50.1233
13 │ 47.7208 48.9828
14 │ 52.5879 54.0585
15 │ 52.1095 53.2797
, , 3 comps
15×2 Named Matrix{Float64}
A ╲ B │ Y1 Y2
──────┼─────────────────
1 │ 49.617 46.8925
2 │ 54.1998 53.7854
3 │ 57.243 52.184
4 │ 42.3451 54.8915
5 │ 46.5996 48.3552
6 │ 54.5413 54.4113
7 │ 49.3049 44.922
8 │ 59.5652 53.9772
9 │ 55.352 51.6175
10 │ 42.636 43.7472
11 │ 51.419 50.3464
12 │ 55.4203 54.5551
13 │ 50.7042 52.888
14 │ 57.3906 60.3452
15 │ 48.0619 47.9816
residuals
, , 1 comps
15×2 Named Matrix{Float64}
A ╲ B │ Y1 Y2
──────┼─────────────────────
1 │ -5.6375 6.93775
2 │ 1.6787 -4.56981
3 │ 6.33061 15.7033
4 │ -6.2465 19.3415
5 │ -6.6608 -8.92132
6 │ -19.5944 -2.91203
7 │ -1.03144 3.0944
8 │ 14.0917 5.12805
9 │ 9.62464 -0.692502
10 │ -10.1259 -21.0345
11 │ 4.04051 -4.63423
12 │ -0.668503 10.5748
13 │ 8.12452 -4.14408
14 │ 9.56563 -3.65442
15 │ -3.49131 -10.2169
, , 2 comps
15×2 Named Matrix{Float64}
A ╲ B │ Y1 Y2
──────┼─────────────────────
1 │ -7.0942 6.85597
2 │ -0.119074 -4.67074
3 │ 1.90335 15.4548
4 │ 4.30789 19.934
5 │ -3.0524 -8.71874
6 │ -19.6713 -2.91635
7 │ -3.08793 2.97895
8 │ 8.46578 4.8122
9 │ 7.07446 -0.835671
10 │ -5.83729 -20.7937
11 │ 1.86344 -4.75645
12 │ -0.634588 10.5767
13 │ 12.7792 -3.88276
14 │ 7.71214 -3.75848
15 │ -4.60946 -10.2797
, , 3 comps
15×2 Named Matrix{Float64}
A ╲ B │ Y1 Y2
──────┼─────────────────────
1 │ -3.61703 11.4075
2 │ -1.19978 -6.08537
3 │ 0.956954 14.216
4 │ 2.45492 17.5085
5 │ -1.39956 -6.55519
6 │ -22.3413 -6.41134
7 │ 1.49508 8.97803
8 │ 6.33479 2.02277
9 │ 6.24797 -1.91754
10 │ -0.836035 -14.2472
11 │ 3.78102 -2.24638
12 │ -4.02027 6.14493
13 │ 9.79578 -7.78799
14 │ 2.90939 -10.0452
15 │ -0.561928 -4.98156
Xvar
1×3 Named Matrix{Float64}
A ╲ B │ Comp 1 Comp 2 Comp 3
──────┼──────────────────────────
Xvar │ 3117.67 2220.14 2047.35
Xtotvar = 9027.521333333332
=====#
