主成分回帰 Principal Component Regression (PCR)

• 昨日の続き
• Rのplsパッケージのpcr()関数(PCR処理)の出力が今一つわからなかったので特異値分解との関係を確認する
• 行列式関係
  • 
  • 特異値分解
  • ,(
  • 
  • 
  • 
  • 
• 主成分回帰では、「主な」成分のみを使うので、それを取り出すことを考える
• 使用する主成分の数をとすると、上記の式のの1:a列部分(、の1:a列部分(),の1:a行ｘ1:a列部分()を使い、次のように表せて
  • 
  • 特異値分解
  • 
  • のように「部分」を取り出す
  • 
  • 
  • 
  • 
• 主成分a個のときの
  • が回帰係数行列
  • がスコア行列
  • がloading行列
• なお、pcr()関数では、の代わりにの列ベクトルの平均によって、標準化したを使っているようだが、そうすると、の特異値分解では、1個の特異値が0となる。それは、特異値を対角成分とした行列の逆行列を作れなくするので、の最大数はの行数(列数も同じ)より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<-t(t(X)-apply(X,2,mean))
Y.st<-t(t(Y)-apply(Y,2,mean))

• • を特異値分解する。標準化してあるために最後の特異値は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個からだんだんに増やしたいろいろなとそれらを用いた被説明成分()を計算する

# 閾値で成分数を決めるとすれば
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
}

• • は等しい

for(i in 1:maxa)print(predYsXB[[i]]-predYsTC[[i]])

• • 成分数を増やすに従い、はに近づき、最後は一致する

for(i in 1:maxa)print(predYsXB[[i]]-Y.st)

• • を確認する

for(i in 1:maxa)print(Bs[[i]]-Vs[[i]]%*%Cs[[i]])