ryamadaの遺伝学・遺伝統計学メモ このページをアンテナに追加 RSSフィード

数学・コンピュータ関連の姉妹ブログ『ryamadaのコンピュータ・数学メモ』
京都大学大学院医学研究科ゲノム医学センター統計遺伝学分野のWiki
講義・スライド
医学生物学と数学とプログラミングの三重学習を狙う学習ツール
駆け足で読む○○シリーズ
ぱらぱらめくるシリーズ
カシオの計算機
オンライン整数列大辞典

2012-11-19 分散と共分散

[][][]共分散

  • こちらで分散・共分散と遺伝率が話題になった
  • 共分散ってなに?と口語で説明したいときに:
  • 「平均は計算したい」
  • 「平均を計算したら、平均を中心としたばらつきは計算したい」→分散
  • 変数が2つあるとき:
  • それぞれの変数に平均は計算できる
  • 2次元空間的な平均は「2次元空間の重心」
  • 「重心は計算したい」
  • 「重心を計算したら、重心を中心としたばらつきは計算したい」→??
  • 個々の点と重心との離れ具合は
    • 距離で測るのもよい
    • 個々の点と重心とを対角線の両端とする長方形の面積とするのもよい
  • 個々の点と重心と距離は、2つの変数のそれぞれの重心の和(ピタゴラスの定理)だから、個々の変数の分散を求めておけば、わかる
  • 個々の点と重心との対角線の両端とする長方形の面積は、別途、計算しておく必要がある
  • これが共分散

2010-11-17 カテゴリカルデータ分割表の固有値分解

ryamada222010-11-17

[][][][][]分割表の「ひらき」

  • 昨日の続き
  • NxMの分割表を考える
  • すべてのサンプルについて、(サンプル数)x(N+M)の表にデータを格納することとする
  • 行数がサンプル数
  • すべての行には2つの1が立っていて、他は0
  • N列のうちのいずれかに1が一つ、残りのM列のうちのいずれかに1が一つ、という制約がある
  • 行について足し合わせると、すべてのサンプルで、その値は2
  • すべてが2で等しいから、N+M列の和に関して分散は0
  • 分散は、各列の分散と、列が作る共分散との和に分解できるから、N+M個の項目に関する分散共分散の和は、『N+M列の和に関する分散』となるが、これが0
  • また、N列のみ、M列のみに着目すると、N列についての和は1、M列に関しても和が1であって、これらはサンプル間での分散が0だから、N列に関する分散共分散行列の和も0であり、M列に関する共分散行列の和も0
  • したがって、N列とM列との共分散が作るNxM「共分散行列」の和も0
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]) # 5軸がなす軸ペアでサンプルをプロット。分離しない
plot(df2[,L]) # 固有値分解後に分離力のあるトップ5軸でのプロットは分離する
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))

2010-11-16 固有値分解による分散共分散行列の変化

ryamada222010-11-16

[][][][]分散が集約されて共分散が小さくなる

  • ryamada本のR7-5.Rでは、固有値分解を実施している
  • 固有値分解では、オリジナルの軸が説明する分散が10人並みなのに対して、固有値分解を施すことにより、説明する分散が大きい方から軸をとっている
  • その図がこれだが

http://www.genome.med.kyoto-u.ac.jp/func-gen-photo/albums/StatGenetTextbook/PartIII-014.jpeg

http://www.genome.med.kyoto-u.ac.jp/func-gen-photo/albums/StatGenetTextbook/PartIII-015.jpeg

  • そのようにして軸を取り直すことで、分散共分散行列が
    • 対角成分(分散)の和は変わらないが、大きな値を持つものと、ほとんど値を持たないものとに分けなおして、大きい順に並べていること
    • 非対角成分(共分散)が小さくなっていることが見て取れる
#偏った集団構成(100人規模の亜集団4つと10人規模の亜集団を20個)で
#100項目のデータを作成
Nm<-100 #項目数
# 亜集団別の人数発生(100人くらいの4亜集団と20人くらいの10亜集団)
Ns<-c(rpois(4,100),rpois(20,10)) 
Npop<-length(Ns) #亜集団数
M<-NULL #全ジェノタイプデータを納める行列
#亜集団別に平均を振ってシミュレーション
for(j in 1:Npop){
 tmpM<-matrix(rep(0,Nm*Ns[j]),ncol=Nm)
 for(i in 1:Nm){ # 項目ごとのループ
  af<-rnorm(1) # 項目の頻度のおよその値
  tmpM[,i]<-rnorm(Ns[j],af) # 亜集団別の頻度
 }
 #全データ行列に格納
 M<-rbind(M,tmpM)
}
# データを標準化
wholemean<-mean(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]) # 5軸がなす軸ペアでサンプルをプロット。分離しない
plot(df2[,L]) # 固有値分解後に分離力のあるトップ5軸でのプロットは分離する
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))