分割表の「ひらき」 - ryamadaの遺伝学・遺伝統計学メモ

昨日の続き
NxMの分割表を考える
すべてのサンプルについて、(サンプル数)x(N+M)の表にデータを格納することとする
行数がサンプル数
すべての行には２つの１が立っていて、他は０
N列のうちのいずれかに１が一つ、残りのM列のうちのいずれかに１が一つ、という制約がある
行について足し合わせると、すべてのサンプルで、その値は２
すべてが２で等しいから、N+M列の和に関して分散は０
分散は、各列の分散と、列が作る共分散との和に分解できるから、N+M個の項目に関する分散共分散の和は、『N+M列の和に関する分散』となるが、これが０
また、N列のみ、M列のみに着目すると、N列についての和は１、M列に関しても和が１であって、これらはサンプル間での分散が０だから、N列に関する分散共分散行列の和も０であり、M列に関する共分散行列の和も０
したがって、N列とM列との共分散が作るNxM「共分散行列」の和も０

N<-4
M<-3
T<-matrix(rpois(N*M,50),N,M)
D<-matrix(0,sum(T),(N+M))
counter<-1
for(i in 1:N){
	for(j in 1:M){
		#for(k in 1:T[i,j]){
		if(T[i,j]>1){
			for(k in 1:T[i,j]){
				D[counter,]<-c(rep(0,(i-1)),1,rep(0,N-i),rep(0,(j-1)),1,rep(0,M-j))
				counter<-counter+1
			}
		}

	}
}
D
cov(D)
image(cov(D))
sum(cov(D))
sum((cov(D)[1:N,(N+1):(N+M)])^2)*sum(T)/(N*M)
sum(cov(D)[1:N,1:N])
sum(cov(D)[(N+1):(N+M),(N+1):(N+M)])

これに前日と同様に固有値分解をしてやる

M<-D
wholemean<-mean(M)
Ns<-nrow(M)
Nm<-ncol(M)

M<-M-wholemean # 全平均が0になるように
mu<-apply(M,2,mean) # 列平均
M<-t(t(M)-mu) # 列平均が0になるように
# 固有値分解
svdout<-svd(M)
M2<-svdout$u%*%diag(svdout$d) # 分解後データ行列
par(mfcol=c(1,2))
# 固有値分解前後をimage()プロット
image(1:sum(Ns),1:Nm,M,xlab="サンプル(大集団→小集団)",ylab="項目")
image(1:sum(Ns),1:Nm,M2,xlab="サンプル(大集団→小集団)",
ylab="PCA後eigen項目")
df1<-as.data.frame(M);df2<-as.data.frame(M2) # データフレーム化
L<-1:5;par(mfcol=c(1,1))
plot(df1[,L]) # ５軸がなす軸ペアでサンプルをプロット。分離しない
plot(df2[,L]) # 固有値分解後に分離力のあるトップ５軸でのプロットは分離する
vM1<-apply(M,2,var) # 項目ごとの分散(分解前)
vM2<-apply(M2,2,var) # 項目ごとの分散(分解後)
ylim<-c(min(vM1,vM2),max(vM1,vM2))
# 固有値分解前の各項目の分散はどれも同じ程度だが
# 固有値分解には、分散の大きいものと小さいものとのコントラストが大きくなっている
plot(vM1,ylim=ylim,type="b")
par(new=TRUE)
plot(vM2,ylim=ylim,type="b",col="red")
library(rgl)
cD<-cov(M) # 分解前の行列の分散共分散行列
cD2<-cov(M2) # 分解後の行列の分散共分散行列
# 対角成分の和(分散の和)は等しい
sum(diag(cD)) 
sum(diag(cD2))
# 非対角成分の和(共分散の和)は違う
sum(cD)-sum(diag(cD))
sum(cD2)-sum(diag(cD2))

open3d()
plot3d(row(cD),col(cD),cD,type="l")
open3d()
plot3d(row(cD2),col(cD2),cD2,type="l")

par(mfcol=c(1,2))
persp(cov(M))
persp(cov(M2))
par(mfcol=c(1,1))