裏 RjpWiki

Julia ときどき R, Python によるコンピュータプログラム,コンピュータ・サイエンス,統計学

Julia に翻訳--225 主成分回帰

2021年05月12日 | ブログラミング
#==========
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
=====#
コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

PVアクセスランキング にほんブログ村

PVアクセスランキング にほんブログ村