主成分回帰 Principal Component Regression (PCR)

  • 昨日の続き
  • Rのplsパッケージのpcr()関数(PCR処理)の出力が今一つわからなかったので特異値分解との関係を確認する
  • 行列式関係
  • 主成分回帰では、「主な」成分のみを使うので、それを取り出すことを考える
  • 使用する主成分の数をaとすると、上記の式のUの1:a列部分(U_aVの1:a列部分(V_a),Dの1:a行x1:a列部分(D_{a,a})を使い、次のように表せて
  • 主成分a個のときの
    • B_aが回帰係数行列
    • T_aがスコア行列
    • V_aがloading行列
  • なお、pcr()関数では、Xの代わりにXの列ベクトルの平均によって、標準化したX_{st}を使っているようだが、そうすると、X_{st}特異値分解では、1個の特異値が0となる。それは、特異値を対角成分とした行列D逆行列を作れなくするので、aの最大数はDの行数(列数も同じ)より1小さくなることに注意
  • まずplsパッケージを使わずにPCRしてみる
    • YとXをランダムに作る
library(pls)
n<-6
m<-8
# n,mの大小によらず動く
X<-matrix(runif(n*m),n,m)
m2<-5
Y<-matrix(runif(n*m2),n,m2)
    • Xの特異値分解はこのままでもできるが、(pcr()の出力と合わせるために)、列ベクトルの平均で標準化する(X_{st},Y_{st})
X.st<-t(t(X)-apply(X,2,mean))
Y.st<-t(t(Y)-apply(Y,2,mean))
    • X_{st}特異値分解する。標準化してあるために最後の特異値は0になっている
svd.out<-svd(X.st)
# 最後の特異値が0
print(svd.out[[1]])
U<-svd.out[[2]]
D<-diag(svd.out[[1]])
V<-svd.out[[3]]
T<-U%*%D
# X.st = T V^Tになっている
X.st-T%*%t(V)
    • PCRでは成分の数aを指定する。それはTの1:a列のみを使うことを意味する。適当に特異値の大きさで区切って採用することもできる。今回は、0でない特異値について、1個からだんだんに増やしたいろいろなB_a,C_a,T_a,V_aとそれらを用いた被説明Y成分(predYXB,predYTC)を計算する
# 閾値で成分数を決めるとすれば
th<-0.00001 # 閾値
maxa<-length(which(svd.out[[1]]>th)) # 成分数
Bs<-Ts<-Vs<-Cs<-predYsXB<-predYsTC<-list()
for(i in 1:maxa){
	tmpa<-i
	T.2<-T[,1:tmpa]
	C<-solve(t(T.2)%*%T.2)%*%t(T.2)%*%Y.st
	if(i==1){
		tmpD<-matrix(1/svd.out[[1]][1],1,1)
		B<-as.matrix(V[,1:(tmpa)])%*%tmpD%*%t(as.matrix(U[,1:(tmpa)]))%*%Y
	}else{
		B<-V[,1:(tmpa)]%*%diag(1/svd.out[[1]][1:(tmpa)])%*%t(U[,1:(tmpa)])%*%Y
	}
	Bs[[i]]<-B
	Ts[[i]]<-T.2
	Cs[[i]]<-C
	Vs[[i]]<-V[,1:tmpa]
	predYsXB[[i]]<-X.st%*%B
	predYsTC[[i]]<-T.2%*%C
}
    • predYXB,predYTCは等しい
for(i in 1:maxa)print(predYsXB[[i]]-predYsTC[[i]])
    • 成分数を増やすに従い、predYXB,predYTCY_{st}に近づき、最後は一致する
for(i in 1:maxa)print(predYsXB[[i]]-Y.st)
    • B_a=V_aC_aを確認する
for(i in 1:maxa)print(Bs[[i]]-Vs[[i]]%*%Cs[[i]])