読者です 読者をやめる 読者になる 読者になる

INFINITY -数学とかプログラミングとか-

統計とプログラムを使って役に立たせたい

TeX用コマンド入力を支援するための辞書をご利用ください。
sanctuary's blogは,適当なことが書いてあります。

PCA

#Copyright(c) Hiroshi Tashiro
data<-source("なんとか.なんとか")$value

data


data=t(data)
data

#
#相関行列
data.cor=cor(t(data))
data.cor

X=data
X
Xvar=function(X){
	return((1/n)*rowSums(X))#各行の和
}

xvar=Xvar(X)
xvar#標本平均ベクトル

n=ncol(data)
n
p=nrow(data)
p

#標本分散共分散行列の初期化
S=matrix(0,nrow=p,ncol=p)
#S
#(X[,i]-xvar)%*%t(X[,i]-xvar)
for(i in 1:n){
	S=S+(X[,i]-xvar)%*%t(X[,i]-xvar)
}
S=(1/n)*S
S

#相関行列を求める
#標本分散共分散行列の初期化
R=matrix(0,nrow=p,ncol=p)

Z=matrix(0,nrow=p,ncol=n)
Z

#標本分散共分散行列の初期化
#基準化

for(i in 1:n){
	for(j in 1:p){
		Z[j,i]=(X[j,i]-xvar[j])/sqrt(S[j,j])
	}
}

Z


for(j in 1:p){
	for(k in 1:p){
		R[j,k]=S[j,k]/(sqrt(S[j,j])*sqrt(S[k,k]))
	}
}
R



#Rの固有値
Revalues=eigen(R)$value
Revalues

#Rの固有ベクトル
Revectors=eigen(R)$vectors
Revectors[,1]=-Revectors[,1]
Revectors



#寄与率
#i番目の寄与率
for(i in 1:p){
	cat(Revalues[i]/sum(Revalues),",")
}


#累積寄与率(cumulative contribution ratio)
cum=0

for(j in 1:p){
	for(i in 1:j){
		cum=cum+Revalues[i]/sum(Revalues)
	}
	cat(cum,",")
	cum=0
}



X

y1=matrix(0,nrow=n,ncol=1)
for(i in 1:n){
	y1[i,]=t(Revectors[,1])%*%Z[,i]
	y2[i,]=t(Revectors[,2])%*%Z[,i]
	cat(y1[i,],y2[i,],"\n")

}
Z
y1

Revectors[,2]
y2=matrix(0,nrow=n,ncol=1)
for(i in 1:n){
	y2[i,]=t(Revectors[,2])%*%Z[,i]
}

y2

y3=matrix(0,nrow=n,ncol=1)
for(i in 1:n){
	y3[i,]=t(Revectors[,3])%*%Z[,i]
}

y3

注意として、固有ベクトルを求めて、主成分の解釈を与えるときに、マイナスを施したほうが解釈を与えやすい時がある。


いまや当たり前だが、PCAを行なってくれる関数がもちろんある。

#主成分分析
data.pc<-princomp(data,cor=TRUE)
data.pc

一発だ!(・∀・)イイネ!!

Call:
princomp(x = data, cor = TRUE)

Standard deviations:
   Comp.1    Comp.2    Comp.3    Comp.4    Comp.5    Comp.6 
1.5819456 1.4247218 1.1509526 0.6719099 0.4384829 0.1855998 
 8  variables and  168 observations.


こんな感じで返してくれる。

また、

summary(data.pc,loadings=TRUE)

とすると、寄与率や固有ベクトルを返してくれる。



y1,y2を第1,2主成分とする。このとき、

plot(y1,y2,main="PC1-PC2 of Data",
	xlim=c(min(y1),max(y1)),
	ylim=c(min(y2),max(y2)),
	xlab="第1主成分",
	ylab="第2主成分",type="n",lwd=2)

text(y1,y2,
labels=abbreviate(row.names(t(data))),cex=0.7,lwd=2)

などとすれば、PC1とPC2のプロットが与えられる。


おしまい。すべておしまい。