まとめ － 主成分分析、PLS関連
本業でPLSに関してちょっとまとめる必要があるので、計算を主成分分析から始めてたどってみた。
以下、過去記事のまとめ。
主成分分析を計算してみる
http://highschoolstudent.hatenablog.com/entry/2013/04/20/101836
主成分分析を、固有値問題、スペクトル分解で解いた。
Rのprcomp()関数を使ってみた。
主成分分析と特異値分解
http://highschoolstudent.hatenablog.com/entry/2013/04/21/125222
特異値分解とそのプロパティ。
主成分分析を特異値分解で解いた。
主成分分析をNIPALSアルゴリズムで解いてみる
http://highschoolstudent.hatenablog.com/entry/2013/04/23/192356
その他主成分分析に関する関連事項。
主成分分析を１主成分ずつ取り出す繰り返し計算で解いた。
主成分分析とPLS
http://highschoolstudent.hatenablog.com/entry/2013/04/29/144826
主成分分析とPLSモデルの比較。
PLSモデルの定義。
PLSの計算
http://highschoolstudent.hatenablog.com/entry/2013/04/30/191405
PLSをNIPALSアルゴリズムで解いた。
Rのchemometricsパッケージのpls2_nipals()関数を使ってみた。
また、pls2_nipals()での予測値がよくわからなかったので、plsパッケージのmvr()関数で予測値の確認を行った。
PLSは登場する成分が非常に多く、複雑である。誤り、追加情報などあれば、都度修正追加を行う。
PLSであるが、エンジニアリングの分野での需要は多いと思う。
PLSは線形モデルである。線形モデルは、基本的な回帰モデル、次に主成分回帰モデル、そしてPLSモデルへとだんだん複雑になって行く。
２次項など含めれば別であるが、これらはすべて変数間の直線関係をとらえる兄弟モデルである。
あてはまりの良さだけを考えるのであれば、非線形モデルであるニューラルとかの方があてはまりは良くなる。
ただし、ニューラルは変数間の関係性が不鮮明ということで、線形モデルのPLSがエンジニアリングの分野では好まれる場合があると思う。
PLSに関する文献であるが、日本語のものはあきらめた方が良いかもしれない。
英語の文献で読みやすいお勧めは以下。
Introduction to Multivariate Statistical Analysis in Chemometrics
- 作者: Kurt Varmuza,Peter Filzmoser
- 出版社/メーカー: CRC Press
- 発売日: 2009/02/20
- メディア: ハードカバー
- クリック: 1回
- この商品を含むブログを見る
タイトルにあるように、PLSはChemometrics（計量化学、ケモメトリックス）の分野で発展した手法である。横長で変数間の相関が強いデータに適用される。
他にも同様なデータを扱う官能評価系の研究者にも利用される。Sensometrics（センソメトリックス）と呼ばれる分野。
例えば、plsライブラリのoliveoilデータなどよくある例ではないだろうか。
head(oliveoil)
chemical.Acidity chemical.Peroxide chemical.K232 chemical.K270 chemical.DK G1 0.73 12.7 1.9 0.139 0.003 G2 0.19 12.3 1.678 0.116 -0.004 G3 0.26 10.3 1.629 0.116 -0.005 G4 0.67 13.7 1.701 0.168 -0.002 G5 0.52 11.2 1.539 0.119 -0.001 I1 0.26 18.7 2.117 0.142 0.001 sensory.yellow sensory.green sensory.brown sensory.glossy sensory.transp sensory.syrup G1 21.4 73.4 10.1 79.7 75.2 50.3 G2 23.4 66.3 9.8 77.8 68.7 51.7 G3 32.7 53.5 8.7 82.3 83.2 45.4 G4 30.2 58.3 12.2 81.1 77.1 47.8 G5 51.8 32.5 8 72.4 65.3 46.5 I1 40.7 42.9 20.1 67.7 63.5 52.2
str(oliveoil)
'data.frame': 16 obs. of 2 variables: $ chemical: AsIs [1:16, 1:5] 0.73 0.19 0.26 0.67 0.52 0.26 0.24 0.3 0.35 0.19 ... ..- attr(*, "dimnames")=List of 2 .. ..$ : chr "G1" "G2" "G3" "G4" ... .. ..$ : chr "Acidity" "Peroxide" "K232" "K270" ... $ sensory : AsIs [1:16, 1:6] 21.4 23.4 32.7 30.2 51.8 40.7 53.8 26.4 65.7 45 ... ..- attr(*, "dimnames")=List of 2 .. ..$ : chr "G1" "G2" "G3" "G4" ... .. ..$ : chr "yellow" "green" "brown" "glossy" ...
PLSの計算
PLSを手計算し、その後関数での結果と比べてみる。
x1 <- c(0.966,0.207,-0.552,-0.931,0.587,1.55,-0.931,0.587,-1.69,0.207)
x2 <- c(0.619,-1.009,-0.358,-0.684,0.944,1.921,0.293,0.293,-1.335,-0.684)
x3 <- c(-0.474,-0.474,-1.684,0.735,-0.474,1.849,-0.474,0.735,-0.474,0.735)
y1 <- c(1.065,-0.821,0.061,-0.957,0.311,1.813,-0.452,0.529,-1.55,0.001)
y2 <- c(-0.767,-0.067,-0.764,0.195,-1.062,1.933,-0.055,0.443,-1.093,1.236)
標準化されているので平均0、分散1。
colMeans(cbind(x1,x2,x3,y1,y2))
x1 x2 x3 y1 y2
-1.665335e-17 0.000000e+00 1.110223e-17 -5.637851e-18 -1.000000e-04
cov(cbind(x1,x2,x3,y1,y2))
x1 x2 x3 y1 y2
x1 1.0005353 0.7521066 0.4558819 0.8955430 0.4747129
x2 0.7521066 1.0000464 0.3995542 0.8875398 0.3523491
x3 0.4558819 0.3995542 0.9998569 0.4023940 0.8534377
y1 0.8955430 0.8875398 0.4023940 0.9997969 0.4593081
y2 0.4747129 0.3523491 0.8534377 0.4593081 1.0000501
引っ付ける。
X <- cbind(x1,x2,x3)
x1 x2 x3
[1,] 0.966 0.619 -0.474
[2,] 0.207 -1.009 -0.474
[3,] -0.552 -0.358 -1.684
[4,] -0.931 -0.684 0.735
[5,] 0.587 0.944 -0.474
[6,] 1.550 1.921 1.849
[7,] -0.931 0.293 -0.474
[8,] 0.587 0.293 0.735
[9,] -1.690 -1.335 -0.474
[10,] 0.207 -0.684 0.735
Y <- cbind(y1,y2)
y1 y2
[1,] 1.065 -0.767
[2,] -0.821 -0.067
[3,] 0.061 -0.764
[4,] -0.957 0.195
[5,] 0.311 -1.062
[6,] 1.813 1.933
[7,] -0.452 -0.055
[8,] 0.529 0.443
[9,] -1.550 -1.093
[10,] 0.001 1.236
この10×3のXデータと、10×2のYデータに対し、NIPALSアルゴリズムでPLSモデルを作成してみる。
第二スコアまで取り出すとする。
まず、第一スコアを取り出す。
X'Yに特異値分解を適用。
svd_res <- svd(t(X)%*%Y)
$d
[1] 14.364201 4.781158
$u
[,1] [,2]
[1,] -0.6282229 -0.2793742
[2,] -0.5787447 -0.4559675
[3,] -0.5199909 0.8450111
$v
[,1] [,2]
[1,] -0.8054403 -0.5926769
[2,] -0.5926769 0.8054403
第一Xスコア（t1）の重み（w1）。
w1 <- svd_res$u[,1]
[1] -0.6282229 -0.5787447 -0.5199909
第一Yスコア（u1）の重み（c1）。
c1 <- svd_res$v[,1]
[1] -0.8054403 -0.5926769
第一Xスコアの計算。
t1 <- X %*% w1
[,1]
[1,] -0.7186306
[2,] 0.7003870
[3,] 1.4296343
[4,] 0.5985436
[5,] -0.6686262
[6,] -3.0469773
[7,] 0.6617790
[8,] -0.9205324
[9,] 2.0807966
[10,] -0.1163741
第一Yスコアの計算。
u1 <- Y %*% c1
[,1]
[1,] -0.4032108
[2,] 0.7009759
[3,] 0.4036733
[4,] 0.6552344
[5,] 0.3789309
[6,] -2.6059077
[7,] 0.3966563
[8,] -0.6886338
[9,] 1.8962283
[10,] -0.7333540
第一Xスコアに対するX負荷量の計算。
p1_t <- solve( t(t1)%*%t1 ) %*% t(t1) %*% X
p1 <- t(p1_t)
[,1]
x1 -0.6238024
x2 -0.6037275
x3 -0.4975260
第一Yスコアに対するY負荷量の計算。
q1_t <- solve( t(u1)%*%u1 ) %*% t(u1) %*% Y
q1 <- t(q1_t)
[,1]
y1 -0.7490939
y2 -0.6692510
第一Xスコアまでで説明されない残差部分。
t1 p1' が、第一XスコアまででのXデータの近似。
Ex1 <- X - t1 %*% t(p1)
x1 x2 x3
[1,] 0.51771653 0.18514296 -0.8315374
[2,] 0.64390304 -0.58615715 -0.1255393
[3,] 0.33980926 0.50510951 -0.9727197
[4,] -0.55762710 -0.32264280 1.0327910
[5,] 0.16990943 0.54033202 -0.8066589
[6,] -0.35071160 0.08145614 0.3330496
[7,] -0.51818070 0.69253416 -0.1447477
[8,] 0.01276975 -0.26275067 0.2770112
[9,] -0.39199419 -0.07876595 0.5612504
[10,] 0.13440558 -0.75425822 0.6771009
第一Yスコアまでで説明されない残差部分を計算する。
第一YスコアまででのYデータの近似だが、u1 q1'ではない。ここで、Xスコアを説明変数、Yスコアを目的変数とした単回帰を行う。
D行列の第一対角成分（d1）を計算。（切片なしの回帰係数）
d1 <- solve(t(t1)%*%t1) %*% t(t1) %*% u1
[,1]
[1,] 0.7653233
Ey1 <- Y - as.numeric(d1) * t1 %*% t(q1)
y1 y2
[1,] 0.653009762 -1.13507785
[2,] -0.419468828 0.29173357
[3,] 0.880607987 -0.03175074
[4,] -0.613855536 0.50157007
[5,] -0.072322735 -1.40446590
以上、第一スコア関連の計算。
第二スコアの取り出しに進む。
第一スコアで行った計算の繰り返しになる。ただし、最後に計算した残差部分に対して計算を行う。
(X残差行列)'(Y残差行列) への特異値分解。
svd_res2 <- svd(t(Ex1)%*%Ey1)
$d
[1] 4.7994939 0.2475326
$u
[,1] [,2]
[1,] -0.2760840 0.7274020
[2,] -0.4590238 -0.6740562
[3,] 0.8444376 -0.1285871
$v
[,1] [,2]
[1,] -0.5199760 0.8541809
[2,] 0.8541809 0.5199760
第二Xスコア（t2）の重み（w2）。
w2 <- svd_res2$u[,1]
[1] -0.2760840 -0.4590238 0.8444376
第二Yスコア（u2）の重み（c2）。
c2 <- svd_res2$v[,1]
[1] -0.5199760 0.8541809
第二Xスコアの計算。
t2 <- Ex1 %*% w2
[,1]
[1,] -0.93009969
[2,] -0.01472129
[3,] -1.14707426
[4,] 1.17418014
[5,] -0.97610760
[6,] 0.34067511
[7,] -0.29705871
[8,] 0.35100196
[9,] 0.61831969
[10,] 0.88088466
第二Yスコアの計算。
u2 <- Ey1 %*% c2
[,1]
[1,] -1.30911119
[2,] 0.46730695
[3,] -0.48501586
[4,] 0.74762169
[5,] -1.16206186
[6,] 0.28365415
[7,] 0.28030389
[8,] -0.02499112
[9,] 0.16241458
[10,] 1.03902459
第二Xスコアに対するX負荷量の計算。
p2_t <- solve( t(t2)%*%t2 ) %*% t(t2) %*% Ex1
p2 <- t(p2_t)
[,1]
x1 -0.2978230
x2 -0.4388790
x3 0.8482805
第二Yスコアに対するY負荷量の計算。
q2_t <- solve( t(u2)%*%u2 ) %*% t(u2) %*% Ey1
q2 <- t(q2_t)
[,1]
y1 -0.3707696
y2 0.9450091
第二Xスコアまでで説明されない残差部分。
Ex2 <- Ex1 - t2 %*% t(p2)
x1 x2 x3
[1,] 0.240711453 -0.223058287 -0.0425519806
[2,] 0.639518696 -0.592618023 -0.1130514849
[3,] -0.001815838 0.001682669 0.0003209963
[4,] -0.207929247 0.192680246 0.0367568772
[5,] -0.120797870 0.111938862 0.0213541507
[6,] -0.249250712 0.230971299 0.0440615157
[7,] -0.606651611 0.562161326 0.1072413768
[8,] 0.117306205 -0.108703266 -0.0207369084
[9,] -0.207844368 0.192601591 0.0367418726
[10,] 0.396753293 -0.367656417 -0.0701364154
D行列の第二対角成分（d2）を計算。
d2 <- solve(t(t2)%*%t2) %*% t(t2) %*% u2
[,1]
[1,] 0.8001358
第二Yスコアまでで説明されない残差部分。
Ey2 <- Ey1 - as.numeric(d2) * t2 %*% t(q2)
y1 y2
[1,] 0.3770808 -0.4317964
[2,] -0.4238361 0.3028649
[3,] 0.5403100 0.8355929
[4,] -0.2655162 -0.3862693
[5,] -0.3619007 -0.6663962
[6,] 0.1672376 0.1147618
[7,] -0.1607300 0.5085755
[8,] 0.1053900 -0.2938957
[9,] -0.1736469 -0.4947633
[10,] 0.1956114 0.5103257
以上、第二スコア関連の計算。第一、第二スコアの取出し完了。
第二スコアまでの結果をひっつけ行列にする。
Xスコア（T）。
T <- cbind(t1, t2)
[,1] [,2]
[1,] -0.7186306 -0.93009969
[2,] 0.7003870 -0.01472129
[3,] 1.4296343 -1.14707426
[4,] 0.5985436 1.17418014
[5,] -0.6686262 -0.97610760
[6,] -3.0469773 0.34067511
[7,] 0.6617790 -0.29705871
[8,] -0.9205324 0.35100196
[9,] 2.0807966 0.61831969
[10,] -0.1163741 0.88088466
Yスコア（U）。
U <- cbind(u1, u2)
[,1] [,2]
[1,] -0.4032108 -1.30911119
[2,] 0.7009759 0.46730695
[3,] 0.4036733 -0.48501586
[4,] 0.6552344 0.74762169
[5,] 0.3789309 -1.16206186
[6,] -2.6059077 0.28365415
[7,] 0.3966563 0.28030389
[8,] -0.6886338 -0.02499112
[9,] 1.8962283 0.16241458
[10,] -0.7333540 1.03902459
X負荷量（P）。
P <- cbind(p1, p2)
[,1] [,2]
x1 -0.6238024 -0.2978230
x2 -0.6037275 -0.4388790
x3 -0.4975260 0.8482805
Y負荷量（Q）。
Q <- cbind(q1, q2)
[,1] [,2]
y1 -0.7490939 -0.3707696
y2 -0.6692510 0.9450091
D行列（D）。
D <- diag(c(d1,d2))
[,1] [,2]
[1,] 0.7653233 0.0000000
[2,] 0.0000000 0.8001358
X重み（W）。
W <- cbind(w1, w2)
w1 w2
[1,] -0.6282229 -0.2760840
[2,] -0.5787447 -0.4590238
[3,] -0.5199909 0.8444376
Y重み（C）。
C <- cbind(c1, c2)
c1 c2
[1,] -0.8054403 -0.5199760
[2,] -0.5926769 0.8541809
NIPALSアルゴリズムでは、最終的な回帰係数（B）は、B=W (P'W)^-1 C' と教科書では表記されていると思う。ただし、このCはY重みのことではないようである。
Y重み（C）ではないので、C*とし、B=W (P'W)^-1 C*' と書き直す。
C*'は、Xスコア（T）でYを回帰したときの係数である。
Y=T C*'+Error
よって、C*'= (T' T)^-1 T' Y
B <- W %*% solve(t(P) %*% W) %*% solve( t(T) %*% T ) %*% t(T) %*% Y
y1 y2
[1,] 0.49447906 0.1088077
[2,] 0.54069342 -0.0396553
[3,] -0.03711641 0.8233860
よって、Yの予測値はXBで計算される。
X %*% B
y1 y2
[1,] 0.8299492 -0.3097233
[2,] -0.4256093 -0.3277496
[3,] -0.4040167 -1.4324472
[4,] -0.8574749 0.5310129
[5,] 0.8182670 -0.3638494
[6,] 1.7364864 1.6149147
[7,] -0.2843437 -0.5032039
[8,] 0.4214018 0.6574398
[9,] -1.5399022 -0.5212301
[10,] -0.2947577 0.6548361
Rの関数を使用する。
chemometricsライブラリのpls2_nipals()関数を利用。この関数はNIPALSアルゴリズム。
plsライブラリのmvr()関数を使用すると、SIMPLSやKernel法を使えるとのこと。ちなみにNIPALSとKernel法は一致するようである。
library(chemometrics)
第二スコアまで取り出す。引数aはスコア数。もともと標準化しているデータなのでscale（標準化）は行わない。
pls2_res <- pls2_nipals(X, Y, a=2, scale = FALSE)
$P
[,1] [,2]
x1 0.6238024 0.2978230
x2 0.6037275 0.4388790
x3 0.4975260 -0.8482805
$T
[,1] [,2]
[1,] 0.7186306 0.93009969
[2,] -0.7003870 0.01472129
[3,] -1.4296343 1.14707426
[4,] -0.5985436 -1.17418014
[5,] 0.6686262 0.97610760
[6,] 3.0469773 -0.34067511
[7,] -0.6617790 0.29705871
[8,] 0.9205324 -0.35100196
[9,] -2.0807966 -0.61831969
[10,] 0.1163741 -0.88088466
$Q
[,1] [,2]
y1 0.7490939 0.3665540
y2 0.6692510 -0.9475753
$U
[,1] [,2]
[1,] 0.4032701 1.25693839
[2,] -0.7009166 -0.41662732
[3,] -0.4036140 0.58855239
[4,] -0.6551752 -0.70432381
[5,] -0.3788716 1.11351346
[6,] 2.6059670 -0.50458887
[7,] -0.3965970 -0.23242261
[8,] 0.6886931 -0.04181580
[9,] -1.8961691 -0.01168087
[10,] 0.7334133 -1.04754497
$D
[1] 0.7653233 0.8001358
$W
[,1] [,2]
x1 0.6282229 0.2760840
x2 0.5787447 0.4590238
x3 0.5199909 -0.8444376
$C
[,1] [,2]
y1 0.8054403 0.5199760
y2 0.5926769 -0.8541809
$B
y1 y2
x1 0.64000890 0.15218608
x2 0.69603317 -0.03463704
x3 -0.02816547 1.04246642
X負荷量。
P_r <- pls2_res$P
[,1] [,2]
x1 0.6238024 0.2978230
x2 0.6037275 0.4388790
x3 0.4975260 -0.8482805
Xスコア。
T_r <- pls2_res$T
[,1] [,2]
[1,] 0.7186306 0.93009969
[2,] -0.7003870 0.01472129
[3,] -1.4296343 1.14707426
[4,] -0.5985436 -1.17418014
[5,] 0.6686262 0.97610760
[6,] 3.0469773 -0.34067511
[7,] -0.6617790 0.29705871
[8,] 0.9205324 -0.35100196
[9,] -2.0807966 -0.61831969
[10,] 0.1163741 -0.88088466
Y負荷量。
Q_r <- pls2_res$Q
[,1] [,2]
y1 0.7490939 0.3665540
y2 0.6692510 -0.9475753
Yスコア。
U_r <- pls2_res$U
[,1] [,2]
[1,] 0.4032701 1.25693839
[2,] -0.7009166 -0.41662732
[3,] -0.4036140 0.58855239
[4,] -0.6551752 -0.70432381
[5,] -0.3788716 1.11351346
[6,] 2.6059670 -0.50458887
[7,] -0.3965970 -0.23242261
[8,] 0.6886931 -0.04181580
[9,] -1.8961691 -0.01168087
[10,] 0.7334133 -1.04754497
これの第二Yスコアが手計算からずれるが、理由はわからない。
D行列。
D_r <- diag(pls2_res$D)
[,1] [,2]
[1,] 0.7653233 0.0000000
[2,] 0.0000000 0.8001358
X重み。
W_r <- pls2_res$W
[,1] [,2]
x1 0.6282229 0.2760840
x2 0.5787447 0.4590238
x3 0.5199909 -0.8444376
Y重み。
C_r <- pls2_res$C
[,1] [,2]
y1 0.8054403 0.5199760
y2 0.5926769 -0.8541809
最終的な回帰係数。
B_r <- pls2_res$B
y1 y2
x1 0.64000890 0.15218608
x2 0.69603317 -0.03463704
x3 -0.02816547 1.04246642
この最終的な回帰係数であるが、上で手計算したものと異なる。この結果は、B=W (P'W)^-1 C' で、CをC*とせず計算したもののようである。なぜこうなっているかは不明。
当然、Yの予測値を X %*% B_r とすることはできない。
念のため、plsパッケージのmvr()関数でPLSでの予測値を出しておく。
library(pls)
NIPALSと結果が一致するとされるkernel法を使用
a1 <- mvr(Y~X, ncomp=2, method="kernelpls")
予測値。
predict(a1)
y1 y2
1 0.8299492 -0.3098233
2 -0.4256093 -0.3278496
3 -0.4040167 -1.4325472
4 -0.8574749 0.5309129
5 0.8182670 -0.3639494
6 1.7364864 1.6148147
7 -0.2843437 -0.5033039
8 0.4214018 0.6573398
9 -1.5399022 -0.5213301
10 -0.2947577 0.6547361
これは手計算に一致する。
PLSはNIPALS、SIMPLSなどのアルゴリズムによっても数値がずれるし、ソフトウェアによっても途中での基準化や制約のかけ方が異なったりで数値がずれると思う。
主成分分析とPLS
主成分分析では、以下の手順で解を求めた。
合成変数をsとする。行数nの列ベクトルである。
s = f_1 x1 + f_2 x2 + … + f_p xp = X f
f_1, f_2, …, f_p が各変数に対応する係数である。f = (f_1, f_2, …, f_p)’を行数pの列ベクトルとする。
sの分散（Var(s)）を求める。
Var(s) = s’s / n-1 = (X f)’ X f / n-1 = f’ X’ X f / n-1 = f’ (X’ X / n-1) f = f’ R f
Xは標準化されている。
R = X’X / n-1 は相関係数行列。
主成分分析では、Var(s) = f’R f が最大になるときのfを求める。
この時のfが、行列Rに対する固有値問題を解いたときの固有ベクトルを並べた行列（V）の第一列目の固有ベクトル。（第一列目の固有ベクトル（v1）が最大固有値（λ1）に対応する。）
もしくは、行列Rを特異値分解したときの特異ベクトルとなる。Rは対称行列なので、左右どちらの特異ベクトルも一致する。
また、最大化問題を解くということなので、Rに対してではなく、n-1を除いたX'Xに対して固有値分解や特異値分解を適用しても同じことである。
PLSでは、X（n×p）側の合成変数だけでなく、Y（n×q）側の合成変数も考える。
g = l_1 y1 + l_2 y2 + … + l_q yq = Y l
g = (l_1, l_2, …, l_q)’を行数qの列ベクトルとする。
sとgの共分散（Cov(s,g)）を求める。
Cov(s,g) = s' g / sqrt(n-1) = (Xf)' Yl / sqrt(n-1) = f'X'Yl / sqrt(n-1)
X'Y（p×q）に特異値分解を適用する。
X'Y = ULV'
U（p×r）、L（r×r）対角行列、V（q×r）。（r = min（rank(X), rank(Y)））
fは左特異ベクトル（U）の一行目、lは右特異ベクトル（V）の一行目となる。
（ここに関してはあいまいなので、後ほど詳しく追加修正したい）
このときのfがX重み（Wの列成分）、lがY重み（Cの列成分）、sがXスコア（Tの列成分（t1））、gがYスコア（Uの列成分（u1））となる。
主成分分析では、Xの合成変数（s）の分散（Var(s)）を最大化する。
PLSでは、Xの合成変数（s）とYの合成変数（g）の共分散（Cov(s,g)=Cor(s,g)SD(s)SD(g)）を最大化する。
（Corは相関、SDは標準偏差）
最大化の問題を解き、Xスコア（T）とYスコア（U）が求まった。
（スコアを一つずつ取り出す繰り返し計算になるがそれは後述。）
T = [t1, t2,…, tr*]、U = [u1, u2,…, ur*]。r* = min(n,p)。
T(n×r*)、U（n×r*）でサイズは同じになる。
取り出すスコアの数がおおければ、PLSモデルの複雑さは増す。
r*は取り出せる最大スコア数であるが、通常r*以下に抑えられる。
PLSモデルを記述する。スコアはkまで取り出すとする。（k < r*）
X = TP' + Ex
Y = UQ' + Ey
X、Yは上記式により、スコアと負荷量と残差に分解されており、そしてXスコアを説明変数、Yスコアを目的変数として回帰が行われる。
U = TD + H
X(n×p)、Y（n×q）がデータ。
T（n×k）がXスコア、U（n×k）がYスコア。Tは直行行列だが、Uには直行の制約がない。
P（k×p）がX負荷量、Q（k×q）がY負荷量。
Ex（n×p）がk個のXスコアでXが説明されなかった残差で、X残差と呼ばれる。Ey（n×p）がk個のYスコアでYが説明されなかった残差部分。
D（k×k）は対角行列。（呼び名は不明、D行列？）
計算では、第一Xスコアと第一Yスコアを単回帰、第二Xスコアと第二Yスコアを単回帰 … と続くので、対角行列の要素は各単回帰の回帰係数。切片は定義されない。
また、最初にXとYの線形結合式を定義したところから計算が始まるので、以下の式も計算上存在する。
T = XW
U = YC
W（p×k）がX重み、C（q×k）がY重み。
XによるYの予測といった形で書くのであれば、以下。
Y = XB + E
B（p×q）は最終的な回帰係数。E（n×q）は残差行列、Y残差と呼ばれる。
X、Y行列の予測値は以下になる。
予測X = TP'
予測Y = XB
PLSは、X側で主成分分析、Y側で主成分分析、各スコア（主成分）で単回帰、これらを同時に計算する形でパラメータを求める。
次回、PLSをNIPALSアルゴリズムで解いてみる。
主成分分析をNIPALSアルゴリズムで解いてみる
主成分分析とは、合成変数の分散を（ある制限の元）最大化する問題を解くことである。
合成変数とは元のデータの各列を説明変数とする線形結合式で、この合成変数の分散が最大になるときの線形結合式の係数が固有ベクトルであり、そのときの合成変数を主成分と呼ぶ。
データ行列をX（n×p）とする。n > pで、ランクr = pとしておく。
X = [x1 x2 … xp] と書く。x1, x2, …, xp それぞれはXを構成する変数で行数nの列ベクトルとなる。平均0、分散1に標準化されているものとする。
合成変数をsとする。行数nの列ベクトルである。
s = f_1 x1 + f_2 x2 + … + f_p xp = X f
f_1, f_2, …, f_p が各変数に対応する係数である。f = (f_1, f_2, …, f_p)’を行数pの列ベクトルとする。
sの分散（Var(s)）を求める。
Var(s) = s’s / n-1 = (X f)’ X f / n-1 = f’ X’ X f / n-1 = f’ (X’ X / n-1) f = f’ R f
x1, x2, …, xp すべてのベクトルの平均値が0なのでsの平均値も0である。
R = X’X / n-1 は相関係数行列。
Var(s) = f’R f が最大になるときのfは、行列Rに対する固有値問題を解いたときの固有ベクトルを並べた行列（V）の第一列目の固有ベクトル。（第一列目の固有ベクトル（v1）が最大固有値（λ1）に対応する。）
V（p×r）であるが、r = pと仮定しているのでV（p×p）。
V = [v1 v2 … vp]と書く。v1, v2, …, vp それぞれはVを構成する列ベクトルで行数p。
固有ベクトルv1, v2, …, vpは、固有値λ1, λ2, …,λpに対応。
sの分散（Var(s)）の最大値がλ1である。
λ1 = v1’ R v1
sの分散が最大となる線形結合式が第一主成分（t1）となる。
第p主成分までを結合した行列をT（n×p）とする。本当はT（n×r）であるが、r = pとしているのでT（n×p）。
T = [t1 t2 … tp]と書く。t1, t2, …, tp それぞれはTを構成する主成分で行数n。
t1 = X v1
前々回、スペクトル分解を用いてR = VDV’の形で分解できると書いた。（対称行列のみ可能）
これをひっくり返してD = V’RVと書くことができる。（V’V = I）
D（p×p）は、固有値λ1, λ2, …,λpを対角要素とする対角行列。本当はD（r×r）であるが、r = pとしているのでD（p×p）。
Dからλ1のみを取り出し、Vからv1のみを取り出すと上で示したように、λ1 = v1’ R v1と書くことができる。
このように、Var(s) = f’R fを最大化するところから出発した問題が、スペクトル分解を用い、λ1 = v1’ R v1と表すことができる。（Rが対象行列なのでスペクトル分解と書いたが、特異値分解でも同じこと）
上で、t1 = X v1と書いたが、すべての主成分を作成した、と考えると、T = XVと書くことができる。
これをひっくり返してX = T V’とも書くことができる。（Vは直行行列なので、V’V = VV’ = I）
Xは、T V’という形に分解できますよ、といった見方。
この場合、X（n×p）、T（n×p）でXの列数pまで主成分を作成していることになる。主成分分析の目的は次元の縮小なので、通常pより小さい次元qまでで主成分の作成はとどめる（p > q）。
この場合、X = T* V*’ + E とデータ行列Xを書くことができる。
T*（n×q）は、主成分qまでを引っ付けた行列。V*（p×q）はq番目までの固有値（λ1, λ2, …,λq）に対する固有ベクトルの行列。
E（n×p）はXがT* V*’ で説明されない部分を示す残差行列。E = X - T* V*’。
主成分分析は、スペクトル分解や特異値分解を用いると一発で解けるのだが、主成分を一つ作成し、残差行列を求め、それに対してまた主成分を一つ作成することで計算することもできる。これがNIPALSアルゴリズム（厳密には違うかもしれないが、残差を取って繰り返し処理していくという意味でNIPALS風）による計算方法である。
Xからスタート。
X
x1 x2 x3
[1,] 0.966 0.619 -0.474
[2,] 0.207 -1.009 -0.474
[3,] -0.552 -0.358 -1.684
[4,] -0.931 -0.684 0.735
[5,] 0.587 0.944 -0.474
[6,] 1.550 1.921 1.849
[7,] -0.931 0.293 -0.474
[8,] 0.587 0.293 0.735
[9,] -1.690 -1.335 -0.474
[10,] 0.207 -0.684 0.735
Xに対し特異値分解を行い、右特異ベクトル（svd(X)$v）の第一列目を第一主成分の係数（固有ベクトル）として取り出す。
v1 <- svd(X)$v[,1]
[1] -0.6254950 -0.6104604 -0.4858951
第一主成分を作成。
t1 <- X %*% v1
[,1]
[1,] -0.75178892
[2,] 0.71679136
[3,] 1.38206543
[4,] 0.64275790
[5,] -0.71312594
[6,] -3.04063179
[7,] 0.63378525
[8,] -0.90316338
[9,] 2.10236553
[10,] -0.06905545
t1 v1’ が第一主成分まででXを近似した行列。
残差（X - t1 v1’）をE1とする。
E1 <- X - t1 %*% t(v1)
x1 x2 x3
[1,] 0.49575977 0.16006263 -0.8392906
[2,] 0.65534944 -0.57142725 -0.1257146
[3,] 0.31247506 0.48569624 -1.0124612
[4,] -0.52895813 -0.29162175 1.0473129
[5,] 0.14094327 0.50866485 -0.8205044
[6,] -0.35190008 0.06481466 0.3715719
[7,] -0.53457047 0.67990081 -0.1660469
[8,] 0.02207579 -0.25834549 0.2961573
[9,] -0.37498080 -0.05158907 0.5475291
[10,] 0.16380616 -0.72615562 0.7014463
残差E1に特異値分解を行い、右特異ベクトル（svd(E1)$v）の第一列目を第二主成分の係数（固有ベクトル）として取り出す。
v2 <- svd(E1)$v[,1]
[1] -0.2845735 -0.4013624 0.8705895
第二主成分を作成。
t2 <- E1 %*% v2
[,1]
[1,] -0.93600078
[2,] -0.06659147
[3,] -1.16530044
[4,] 1.17935314
[5,] -0.95859021
[6,] 0.39761389
[7,] -0.26532069
[8,] 0.35523946
[9,] 0.60408863
[10,] 0.85550848
[t1 t2] [v1’ v2’]’ が第二主成分まででXを近似した行列。
残差（X - [t1 t2] [v1’ v2’]’）をE2とする。計算ではE1 – t2 v2’ とする。
E2 <- E1 - t2 %*% t(v2)
x1 x2 x3
[1,] 0.22939874 -0.21561292 -0.024418072
[2,] 0.63639927 -0.59815456 -0.067740753
[3,] -0.01913857 0.01798843 0.002037182
[4,] -0.19334547 0.18172628 0.020580425
[5,] -0.13184611 0.12392276 0.014034200
[6,] -0.23874970 0.22440193 0.025413424
[7,] -0.61007371 0.57341105 0.064938561
[8,] 0.12316753 -0.11576572 -0.013110419
[9,] -0.20307318 0.19086941 0.021615880
[10,] 0.40726121 -0.38278666 -0.043350429
残差E2に特異値分解を行い、右特異ベクトル（svd(E2)$v）の第一列目を第三主成分の係数（固有ベクトル）として取り出す。
v3 <- svd(E2)$v[,1]
[1] -0.72648048 0.68282230 0.07732934
第三主成分を作成。
t3 <- E2 %*% v3
[,1]
[1,] -0.31576725
[2,] -0.87600326
[3,] 0.02634424
[4,] 0.26613994
[5,] 0.18148610
[6,] 0.32863884
[7,] 0.83976615
[8,] -0.16954004
[9,] 0.27953013
[10,] -0.56059485
[t1 t2 t3] [v1’ v2’ v3’]’ が第二主成分まででXを近似した行列。（近似と言うか、Xに一致）
残差（X - [t1 t2 t3] [v1’ v2’ v3’]’）をE3とする。計算ではE2 – t3 v3’ とする。
E3 <- E2 - t3 %*% t(v3)
x1 x2 x3
[1,] 0.000000e+00 2.775558e-17 8.326673e-17
[2,] 0.000000e+00 0.000000e+00 0.000000e+00
[3,] 4.510281e-17 1.734723e-17 2.506675e-16
[4,] 0.000000e+00 -2.775558e-17 -2.393918e-16
[5,] -5.551115e-17 -6.938894e-17 5.724587e-17
[6,] 8.326673e-17 5.551115e-17 1.040834e-17
[7,] -1.110223e-16 -1.110223e-16 1.387779e-17
[8,] 1.387779e-17 2.775558e-17 -4.683753e-17
[9,] 0.000000e+00 -2.775558e-17 -2.428613e-17
[10,] -5.551115e-17 1.110223e-16 -9.714451e-17
第三主成分まで取り出すとXをすべて再現できるので、E3はもちろん0。
私は、あいまいになるほどと思っている程度で、詳しくは理解していない。
が、PLSまでたどり着きたいので記しておく。PLSでもこれと同じ方法でスコアを作成することができる。
誤りがあれば、都度修正したい。
特異値分解は、どんな行列も非常に性質の良い3つの行列（ULV’）に分解する手法である。U, L, Vの列を元の行列のランク以下までしか使用しないと、ULV’ は元の行列の最小２乗近似となる。
もやっとしか理解できていないが、多変量解析ではとにかく使いやすいと思う。
主成分分析と特異値分解
―― 特異値分解 ――
A （n×p）をランクrの行列とする。
A = ULV’ と分解することができる。
ここで、U（n×r）とV（p×r）はそれぞれ列ベクトルが直行する。U’U = V’V = I、I（r×r）は単位行列。
L（r×r）は、正の値を対角要素に取る対角行列となる。
L = diag(ψ1, ψ2, …, ψr)、ψ1 > ψ2 > … > ψr。
Uを左特異ベクトルと呼ぶ。
Vを右特異ベクトルと呼ぶ。
Lの対角要素（ψ1, ψ2, …, ψr）を特異値と呼ぶ。
――――――――――――
Uの各列ベクトルは、Aの列ベクトルが張る空間の正規直交基底である。
Vの各列ベクトルは、Aの行ベクトルが張る空間の正規直交基底である。
A’A を考える。
A’A（p×p）は対称行列となる。
A’A = ( ULV’ )’ ULV’ = VL’U’ ULV’ = VL’LV’ = V(L^2)V’
これは、行列A’Aのスペクトル分解と見ることができる。
よって、V、L^2 は、行列A’Aの固有値問題を解いた時の固有ベクトル、対角要素に固有値を持つ対角行列となる。
AA’を考える。
AA’（n×n）は対称行列となる。
AA’ = ULV’(ULV’)’ = ULV’VL’U’ = ULL’U’ = U(L^2)U’
これは、行列AA’のスペクトル分解と見ることができる。
よって、U、L^2は、行列AA’の固有値問題を解いた時の固有ベクトル、対角要素に固有値を持つ対角行列となる。
A’A、AA’それぞれにおいて、固有値（L^2の対角要素）は共通である。
L^2（r×r）の対角要素はψ1, ψ2, …, ψrの2乗。
では、主成分分析で考える。
前回、標準化された行列（X）から相関係数行列（R）を計算した。標準化されているので相関係数行列と共分散行列は一致する。
R = X’X / n-1
nはXの行数。
ここで、A = X / sqrt(n-1) と置く。
そうすると、A’A = X’X / n-1 = R となる。
特異値分解を用いると、A = ULV’ なので、V、L^2は、行列A’A = Rの固有値問題を解いた時の固有ベクトル、対角要素に固有値を持つ対角行列となる。
A <- X/sqrt(n-1)
Rではsvd()が特異値分解を行える関数である。
Svd <- svd(A)
$d
[1] 1.4452338 0.8164783 0.4950769
$u
[,1] [,2] [,3]
[1,] -0.17339499 0.38212929 0.21260483
[2,] 0.16532305 0.02718646 0.58980952
[3,] 0.31876398 0.47574259 -0.01773747
[4,] 0.14824773 -0.48147970 -0.17919096
[5,] -0.16447764 0.39135159 -0.12219387
[6,] -0.70130101 -0.16232883 -0.22127123
[7,] 0.14617825 0.10831915 -0.56541121
[8,] -0.20830848 -0.14502916 0.11415064
[9,] 0.48489629 -0.24662368 -0.18820653
[10,] -0.01592717 -0.34926771 0.37744628
$v
[,1] [,2] [,3]
[1,] -0.6254950 0.2845735 0.72648048
[2,] -0.6104604 0.4013624 -0.68282230
[3,] -0.4858951 -0.8705895 -0.07732934
V <- Svd$v
[,1] [,2] [,3]
[1,] -0.6254950 0.2845735 0.72648048
[2,] -0.6104604 0.4013624 -0.68282230
[3,] -0.4858951 -0.8705895 -0.07732934
固有ベクトルである。
L <- diag( Svd$d )
L^2
[,1] [,2] [,3]
[1,] 2.088701 0.0000000 0.0000000
[2,] 0.000000 0.6666368 0.0000000
[3,] 0.000000 0.0000000 0.2451012
固有値を対角要素とする対角行列である。
このように、主成分分析を行列Rに対する固有値問題を解き求めるように、特異値分解でも求めることが可能である。
日本語書籍で特異値分解まで取り上げているものは少ないが、以下などに記載あり。
- 作者: 岩崎学,吉田清隆
- 出版社/メーカー: 東京図書
- 発売日: 2006/05
- メディア: 単行本
- この商品を含むブログを見る
- 作者: 足立浩平,村上 隆
- 出版社/メーカー: 朝倉書店
- 発売日: 2011/08/25
- メディア: 単行本（ソフトカバー）
- この商品を含むブログを見る
主成分分析を計算してみる
主成分分析を計算してみる。
x1 <- c(0.966,0.207,-0.552,-0.931,0.587,1.55,-0.931,0.587,-1.69,0.207)
x2 <- c(0.619,-1.009,-0.358,-0.684,0.944,1.921,0.293,0.293,-1.335,-0.684)
x3 <- c(-0.474,-0.474,-1.684,0.735,-0.474,1.849,-0.474,0.735,-0.474,0.735)
ひっつけて行列にする。
X <- cbind(x1, x2, x3)
x1 x2 x3
[1,] 0.966 0.619 -0.474
[2,] 0.207 -1.009 -0.474
[3,] -0.552 -0.358 -1.684
[4,] -0.931 -0.684 0.735
[5,] 0.587 0.944 -0.474
[6,] 1.550 1.921 1.849
[7,] -0.931 0.293 -0.474
[8,] 0.587 0.293 0.735
[9,] -1.690 -1.335 -0.474
[10,] 0.207 -0.684 0.735
3ベクトルとも平均0、分散1に標準化されている。
mean(x1)
mean(x2)
mean(x3)
[1] -1.664657e-17 [1] 0 [1] 1.110223e-17
var(x1)
var(x2)
var(x3)
[1] 1.000535 [1] 1.000046 [1] 0.9998569
この10行3列の行列Xに主成分分析を適用し、3つの主成分を作成する。
最初から標準化している行列 を用いるので、相関係数行列から実施しようと、共分散行列から実施しようが結果は同じとなる。
行数（n=10）と列数（p=3）を定義しておく。
n <- nrow(X)
p <- ncol(X)
相関係数行列を求める。
cor(X)
x1 x2 x3
x1 1.0000000 0.7518879 0.4557925
x2 0.7518879 1.0000000 0.3995735
x3 0.4557925 0.3995735 1.0000000
X'X / (n-1) でも計算できる。
これは共分散行列の計算だが、標準化データなので相関係数行列と一致。
相関係数行列をRと定義
R <- t(X) %*% X / (n-1)
x1 x2 x3
x1 1.0005353 0.7521066 0.4558819
x2 0.7521066 1.0000464 0.3995542
x3 0.4558819 0.3995542 0.9998569
主成分分析は、合成変数（求める主成分）の分散を最大化する問題を解くに等しい。
これは、行列計算の性質を活かし、行列Rに対する固有値問題を解くことに等しくなる。
Eig <- eigen(R)
固有値が合成変数の分散となる。
固有値を対角要素とした対角行列を作成。
D <- diag(Eig$values)
[,1] [,2] [,3]
[1,] 2.088701 0.0000000 0.0000000
[2,] 0.000000 0.6666368 0.0000000
[3,] 0.000000 0.0000000 0.2451012
順にλ1、λ2、λ3とする。（λ1 > λ2 > λ3）主成分1、2、3の分散に対応する。
各固有値に対応する 固有ベクトル。
これは、合成変数（主成分）の係数となる。Vと定義する。
V <- Eig$vectors
[,1] [,2] [,3]
[1,] 0.6254950 -0.2845735 0.72648048
[2,] 0.6104604 -0.4013624 -0.68282230
[3,] 0.4858951 0.8705895 -0.07732934
各固有ベクトルは長さ1に基準化されている。
sum(V[,1]^2)
sum(V[,2]^2)
sum(V[,3]^2)
[1] 1 [1] 1 [1] 1
主成分スコアを計算する。Tと定義する。T = XV。
T <- X %*% V
[,1] [,2] [,3]
[1,] 0.75178892 -0.93600078 0.31576725
[2,] -0.71679136 -0.06659147 0.87600326
[3,] -1.38206543 -1.16530044 -0.02634424
[4,] -0.64275790 1.17935314 -0.26613994
[5,] 0.71312594 -0.95859021 -0.18148610
[6,] 3.04063179 0.39761389 -0.32863884
[7,] -0.63378525 -0.26532069 -0.83976615
[8,] 0.90316338 0.35523946 0.16954004
[9,] -2.10236553 0.60408863 -0.27953013
[10,] 0.06905545 0.85550848 0.56059485
各主成分の分散が固有値
var(T[,1])
var(T[,2])
var(T[,3])
[1] 2.088701 [1] 0.6666368 [1] 0.2451012
固有値問題を解いたので、RV = VD となる。
R %*% V
V %*% D
[,1] [,2] [,3]
x1 1.306472 -0.1897072 0.17806122
x2 1.275069 -0.2675629 -0.16736054
x3 1.014889 0.5803670 -0.01895351
一致。
固有値問題の性質により、対称行列の固有値（Dの対角要素）はすべて実数となり、各固有ベクトル（V）は直行となる。
直行なので V' = V^-1。
t(V)
solve(V)
[,1] [,2] [,3]
[1,] 0.6254950 0.6104604 0.48589510
[2,] -0.2845735 -0.4013624 0.87058953
[3,] 0.7264805 -0.6828223 -0.07732934
一致。
V'V = VV' = I とも書ける。
t(V) %*% V
[,1] [,2] [,3]
[1,] 1.000000e+00 -1.447410e-17 -2.558988e-16
[2,] -1.447410e-17 1.000000e+00 -1.906841e-17
[3,] -2.558988e-16 -1.906841e-17 1.000000e+00
この性質があるので、各主成分間（Tの各列）の相関が0になる（T' T の非対角要素が0になる）。
Tの相関行列は I。Iは単位行列。
cor(T)
[,1] [,2] [,3]
[1,] 1.000000e+00 1.957105e-17 5.024524e-16
[2,] 1.957105e-17 1.000000e+00 -2.224671e-16
[3,] 5.024524e-16 -2.224671e-16 1.000000e+00
T' T の非対角要素が0になる。
t(T) %*% T
[,1] [,2] [,3]
[1,] 1.879831e+01 2.218583e-16 3.128228e-15
[2,] 2.218583e-16 5.999731e+00 -7.938528e-16
[3,] 3.128228e-15 -7.938528e-16 2.205911e+00
分散とn-1で割るとIとなる。
t(T) %*% T %*% solve(D) / (n-1)
[,1] [,2] [,3]
[1,] 1.000000e+00 3.697804e-17 1.418112e-15
[2,] 1.180203e-17 1.000000e+00 -3.598754e-16
[3,] 1.664101e-16 -1.323147e-16 1.000000e+00
スペクトル分解で考えてみる。
R V = V D から変形（VV'=I）して、 R = VDV' と書くことができる。
V %*% D %*% t(V)
[,1] [,2] [,3]
[1,] 1.0005353 0.7521066 0.4558819
[2,] 0.7521066 1.0000464 0.3995542
[3,] 0.4558819 0.3995542 0.9998569
相関係数行列（R）に一致。
スペクトル分解は行列の近似（元の行列と近似行列の各要素の差の2乗和を最小にする）と考えられる。
主成分分析で言えば、相関係数行列（R）を近似する行列を与えることができる。
V = [v1 v2 v3]とする。 v1、v2、v3、は主成分1、2、3のスコアベクトル。
R = V D V' = λ1 v1' v1 + λ2 v2' v2 + λ3 v3' v3と書けるので、幾つの固有値まで採用するか（幾つの主成分まで取り出すか）で近似の精度が決まる。
すべて取り出すと、先程示したように元の行列と完全に一致する。
第二主成分までだと。R = V D V' = λ1 v1' v1 + λ2 v2' v2。
V[,1:2] %*% D[1:2,1:2] %*% t(V[,1:2])
[,1] [,2] [,3]
[1,] 0.8711773 0.8736907 0.4696512
[2,] 0.8736907 0.8857689 0.3866123
[3,] 0.4696512 0.3866123 0.9983912
Rから遠からず、といった感じである。
Rの関数を使って主成分分析を実施してみる。prcomp()関数を使う。
princomp()という関数もあるようで、これは分散の計算の際nで割るそう。対してprcomp()はn-1で割るらしい。prcomp()の方が一般的ではなかろうか、、。
上の計算でも行列Rを定義した再にn-1で割ったので、prcomp()で実施。
pca_result <- prcomp(X)
summary(pca_result)
Importance of components:
PC1 PC2 PC3
Standard deviation 1.4452 0.8165 0.49508
Proportion of Variance 0.6961 0.2222 0.08169
Cumulative Proportion 0.6961 0.9183 1.00000
固有値。
pca_result$sdev^2
[1] 2.0887007 0.6666368 0.2451012
固有ベクトル。
pca_result$rotation
PC1 PC2 PC3
x1 -0.6254950 0.2845735 0.72648048
x2 -0.6104604 0.4013624 -0.68282230
x3 -0.4858951 -0.8705895 -0.07732934
ベクトルの向きが上の計算と反対になるのは仕方がない。
以下は、JMPによる出力。素晴らしい。
