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のプロットが与えられる。
おしまい。すべておしまい。