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

分布の分布を推定したい

R 二項分布 ベータ分布 多項分布 ディリクレ分布 kNN
  • 複数の因子が多項(二項を含む)の生起確率に影響を与えるらしい
  • サンプルについて因子を観察して生起を観察し、それに基づいて、新規サンプルについて因子を測定して、生起を予測したいというような場合
  • パラメトリックにやるのもありだけれど、高次元化したりするとたちまち厄介になるのでノンパラでやってみたい
  • 高次元空間のpoint cloudの密度推定にkNNを用いたことがある
  • それを使う
  • 近隣について多項の観測度数を数え上げていく
  • 近傍に関して生起確率ベクトルが一定であるとみなすなら、k近傍までの度数から推定されるディリクレ分布の信頼区間はkを増やすことで狭くなるが、信頼区間がkによって外れることは少ないだろうことを基にする
  • もしk近傍のkが大きくなり、もはや生起確率ベクトルが一定でないとすると、そのときには、k近傍までの集計の度数からの推定がそれまでのkの推定の信頼区間から外れそうだから、それをもって、生起確率ベクトルが一定であるだろうとするkの上限とすることにする
  • このようにkの値を大きくとることで推定するディリクレ分布のピークを大きくする
  • やってみる
  • 空間に偏りをもって分布する母集団があり、その空間には多項生起確率ベクトルの場があるものとする
  • 左上は観測標本における、「真の生起確率」
  • 左中は、観測標本の「生起」
  • 左下は、知りたい点(これも母集団分布に従っている)における「真の生起確率」
  • 中上は、知りたい点の「推定生起確率」
  • 中中は、知りたい点の「真と推定の相関」
  • 中下は、知りたい点の「真と推定の確率の差」をドットの大小で表したもの
  • 右上と右中は、知りたい点の真の値と、推定値(期待値とそのベータ―分布の2.5,97.5パーセンタイル)。右上は真の値でソートしたもの、右中は推定期待値でソートしたもの

  • もろもろ関数
library(MCMCpack)
# ベータ分布のzの周りの分散
my.beta.2moment <- function(a,b,z){
	a*(a+1)/(a+b)/(a+b+1)-2*z*a/(a+b)+z^2
}
my.dirichlet.mean <- function(as){
	as/sum(as)
}
my.dirichlet.var <- function(as){
	sum.as <- sum(as)
	as * (sum.as-as)/(sum.as^2 * (sum.as+1))
}
my.dirichlet.2moment <- function(as,t=0){
	m <- my.dirichlet.mean(as)
	my.dirichlet.var(as)+m^2-2*m*t+t^2
}
my.Dirichlet.knn.conf <- function(D,X,Y,Y.level,mm="est",r.n=100,q=0.95,n.mc=100,est=FALSE,self=TRUE){
	# ベータ分布クオンタイル
	qs <- sort(c((1-q)/2,1-(1-q)/2))
	# Xの説明変数の数
	n.x <- length(X[1,])
	# Yのレベル数
	n.y <- length(Y.level)
	# サンプル数
	n.sample <- length(Y)
	# 予想対象数
	n.cnt <- length(D[,1])
	# Xの距離行列
	D.mat <- t(apply(D,1,my.dist.mat,X))
	rs <- seq(from=0,to=max(D.mat),length=r.n)
	if(!self){
		D.mat[which(D.mat==0)] <- max(D.mat)+1
	}
	D.ord <- t(apply(D.mat,1,order))
	D.level.ord <- array(0,c(n.cnt,n.sample,n.y+1))
	for(i in 1:n.cnt){
		for(j in 1:n.y){
			D.level.ord[i,which(Y[D.ord[i,]]==Y.level[j]),j] <- 1
		}
		D.level.ord[i,,n.y+1] <- 1
		D.level.ord[i,,] <- apply(D.level.ord[i,,],2,cumsum)
	}
	Est <- Mode <- array(0,c(n.cnt,n.sample,n.y))
	for(i in 1:n.cnt){
		Est[i,,] <- (D.level.ord[i,,1:n.y]+1)/(D.level.ord[i,,n.y+1]+n.y)
		Mode[i,,] <- (D.level.ord[i,,1:n.y])/D.level.ord[i,,n.y+1]
	}
	K <- rep(1,n.cnt)
	if(mm=="mode"){
		MM <- Mode
	}else{
		MM <- Est
	}
	for(i in 1:n.cnt){
		loop <- TRUE
		for(j in 2:n.sample){
			for(j2 in 1:j){
				if(n.y==2){
					tmp.LU <- qbeta(qs,D.level.ord[i,j2,1]+1,D.level.ord[i,j2,2]+1)
					#print(D.level.ord[i,j2,1:n.y])
					#print(tmp.LU)
					#print(Est[i,j,1])
					if(prod(tmp.LU-MM[i,j,1])>0){
						#print("break")
						loop <- FALSE
						break
						
					}else{
						K[i] <- j2
					}
				}else{
					tmp.q <- my.qdirichlet(MM[i,j,],D.level.ord[i,j2,1:n.y]+1,n.mc)
					if(tmp.q>q){
						break
					}else{
						K[i] <- j2
					}
				}
				# EstでなくModeか？
				#K[i] <- j2
			}
			if(!loop){
				break
			}
		}
	}
	return(list(X=X,D=D,Y=Y,n.sample=n.sample,n.cnt=n.cnt,D.level.ord=D.level.ord,Est=Est,Mode=Mode,K=K))
}
my.qdirichlet <- function(v,a,n=10000){
	r <- rdirichlet(n,a)
	d <- ddirichlet(r,a)
	d0 <- ddirichlet(v,a)
	length(which(d>=d0))/n
}
my.Dirichlet.knn <- function(X,Y,Y.level,r.n=100,est=FALSE,self=TRUE){
	# Xの説明変数の数
	n.x <- length(X[1,])
	# Yのレベル数
	n.y <- length(Y.level)
	# サンプル数
	n.sample <- length(Y)
	# Xの距離行列
	D.mat <- as.matrix(dist(X))
	rs <- seq(from=0,to=max(D.mat),length=r.n)
	if(!self){
		D.mat[which(D.mat==0)] <- max(D.mat)+1
	}
	D.ord <- t(apply(D.mat,1,order))
	D.level.ord <- array(0,c(n.sample,n.sample,n.y+1))
	for(i in 1:n.sample){
		for(j in 1:n.y){
			D.level.ord[i,which(Y[D.ord[i,]]==Y.level[j]),j] <- 1
		}
		D.level.ord[i,,n.y+1] <- 1
		D.level.ord[i,,] <- apply(D.level.ord[i,,],2,cumsum)
	}
	
	
	# 距離が定める仮説の対数尤度
	LL.r <- rep(0,length(rs))
	D.within <- matrix(0,length(rs),n.sample)
	for(i in 1:length(rs)){
		tmp <- apply((D.mat <= rs[i]),2,sum)
		D.within[i,] <- tmp
	}
	# rを基準にカウント数
	cnt.r <- array(0,c(n.sample,length(rs),n.y+1))
	for(i in 1:n.sample){
		for(j in 1:length(rs)){
			tmp <- D.within[j,i]
			if(tmp>0){
				cnt.r[i,j,] <- D.level.ord[i,D.within[j,i],]
			}
			
		}
	}
	for(i in 1:length(rs)){
		for(j in 1:n.sample){
			tmp.level <- which(Y.level==Y[j])
			#LL.r[i] <- LL.r[i] + log((D.level.ord[j,D.within[i,j],tmp.level]+1)/(D.level.ord[j,D.within[i,j],n.y+1]+n.y))
			#LL.r[i] <- LL.r[i] + log((cnt.r[j,i,tmp.level]+1)/(cnt.r[j,i,n.y+1]+n.y))
			#LL.r[i] <- LL.r[i] + sum(cnt.r[j,i,1:n.y]*log((cnt.r[j,i,1:n.y]+1)/(cnt.r[j,i,n.y+1]+n.y)))+lgamma(cnt.r[j,i,n.y+1]+1)-sum(lgamma(cnt.r[j,i,1:n.y]+1))
			v <- cnt.r[j,i,1:n.y]
			LL.r[i] <- LL.r[i] + log(dmultinom(v/sum(v),prob=(v+1)/sum(v+length(v))))
		}
	}
	#LL.r.max <- which(LL.r == max(LL.r))
	#if(length(LL.r.max)==1 & LL.r.max[1] == 1){
	#	L.r <- c(rep(0,length(rs)-1),1)
	#}else{
		# 事前確率を考慮？
		#LL.r <- LL.r + c(rep(1/(length(rs)-1),length(rs)-1),1)/2
		# 仮説の相対確率
		#LL.r[1] <- min(LL.r[1])
		L.r <- exp(LL.r-max(LL.r))
		#L.r[1] <- 0
		L.r <- L.r/sum(L.r)
	#}
	# 期待値

	Est <- Var <- c()
	if(est){
		Est <- matrix(0,n.sample,n.y)

		for(i in 1:n.sample){
			for(j in 1:length(rs)){
				Est[i,] <- Est[i,] + L.r[j] * my.dirichlet.mean(cnt.r[i,j,1:n.y]+1)
			}
		}
		# 期待値の周囲の分散
		Var <- matrix(0,n.sample,n.y)
		for(i in 1:n.sample){
			for(j in 1:length(rs)){
				Var[i,] <- Var[i,] + L.r[j] * my.dirichlet.2moment(cnt.r[i,j,1:n.y]+1,Est[i,])
			}
		}

	}
	return(list(X=X,Y=Y,L.r=L.r,rs=rs,Est=Est,Var=Var,cnt.r=cnt.r,n.x=n.x,n.y=n.y,n.sample=n.sample))
}
my.dist.mat <- function(x,X){
	d <- length(x)
	tmp <- t(t(X)-x)
	sqrt(apply(tmp^2,1,sum))
}

my.dirichlet.est <- function(X.r,d.knn.out){
	Est <- Var <- c()
	X <- X.r
	n.sample <- d.knn.out$n.sample
	n.pnt <- length(X[,1])
	n.y <- d.knn.out$n.y
	rs <- d.knn.out$rs
	L.r <- d.knn.out$L.r
	cnt.r <- d.knn.out$cnt.r
	Est <- matrix(0,n.pnt,n.y)
	# Xの距離行列
	#D.mat <- as.matrix(dist(X))
	D.mat <- t(apply(X,1,my.dist.mat,d.knn.out$X))
	
	D.ord <- t(apply(D.mat,1,order))
	D.level.ord <- array(0,c(n.pnt,n.sample,n.y+1))
	for(i in 1:n.pnt){
		for(j in 1:n.y){
			D.level.ord[i,which(Y[D.ord[i,]]==Y.level[j]),j] <- 1
		}
		D.level.ord[i,,n.y+1] <- 1
		D.level.ord[i,,] <- apply(D.level.ord[i,,],2,cumsum)
	}
	#rs <- seq(from=0,to=max(D.mat),length=r.n)
	# 距離が定める仮説の対数尤度
	LL.r <- rep(0,length(rs))
	D.within <- matrix(0,length(rs),n.pnt)
	for(i in 1:length(rs)){
		tmp <- apply((D.mat <= rs[i]),1,sum)
		D.within[i,] <- tmp
	}
	# rを基準にカウント数
	cnt.r <- array(0,c(n.pnt,length(rs),n.y+1))
	for(i in 1:n.pnt){
		for(j in 1:length(rs)){
			tmp <- D.within[i,j]
			if(tmp>0){
				cnt.r[i,j,] <- D.level.ord[i,D.within[i,j],]
			}
			
		}
	}

	for(i in 1:n.pnt){
		for(j in 1:length(rs)){
			Est[i,] <- Est[i,] + L.r[j] * my.dirichlet.mean(cnt.r[i,j,1:n.y]+1)
		}
	}
	# 期待値の周囲の分散
	Var <- matrix(0,n.pnt,n.y)
	for(i in 1:n.pnt){
		for(j in 1:length(rs)){
			Var[i,] <- Var[i,] + L.r[j] * my.dirichlet.2moment(cnt.r[i,j,1:n.y]+1,Est[i,])
		}
	}
	return(list(Est=Est,Var=Var))

}
  • お試し
# 説明変数次元
n.x <- 3
# 被説明変数に影響を及ぼす変数数
n.x.a <- 2
# 被説明変数レベル数
n.y <- 2
Y.level <- 1:n.y
# 標準的な分布を定める
# 説明変数空間分布を定める
# 構成分布の数
n.d <- 3
# 分布の構成比
library(MCMCpack)
r <- rdirichlet(1,rep(1,n.d))
# 分布は正規分布
k.ctr <- 10
k.sd <- k.ctr/2
ctrs <- matrix(rnorm(n.d*n.x)*k.ctr,ncol=n.x)
sds <- matrix(runif(n.d*n.x)*k.sd,ncol=n.x)

# 確率密度分布関数を定める
# 原点からの距離Lに応じた三角関数 k1 * (sin(k2 * L)+1)/2 + k3 ; k1 <=1, k1+k3<=1
k1 <- runif(1)
K1 <- 0.8
k2 <- 1
k3 <- runif(1) * (1-k1)

my.calc.p <- function(x,k1,k2,k3){
	k1 * (sin(k2*sqrt(sum(x^2)))+1)/2+k3
}

# サンプリングする
# サンプル総数
n.s <- 500

X <- matrix(0,0,n.x)
Y <- rep(0,n.s)
r.d <- c(rmultinom(1,size=n.s,prob=r))
for(i in 1:n.d){
	tmp <- matrix(0,r.d[i],0)
	for(j in 1:n.x){
		tmp <- cbind(tmp,rnorm(r.d[i],ctrs[i,j],sds[i,j]))
	}
	X <- rbind(X,tmp)
}
P <- apply(X[,1:n.x.a],1,my.calc.p,k1,k2,k3)
for(i in 1:n.s){
	Y[i] <- sample(1:n.y,1,prob=c(P[i],1-P[i]))
}
plot(X[,1:2],col=Y)
# 推定する
est <- FALSE
#est.out <- my.Dirichlet.knn(X,Y,Y.level,r.n=100,est)
#plot(est.out$L.r,type="l")

# 座標を与える
n.r.pt <- 800

X.random <- matrix(0,0,n.x)
r.d <- c(rmultinom(1,size=n.r.pt,prob=r))
for(i in 1:n.d){
	tmp <- matrix(0,r.d[i],0)
	for(j in 1:n.x){
		tmp <- cbind(tmp,rnorm(r.d[i],ctrs[i,j],sds[i,j]))
	}
	X.random <- rbind(X.random,tmp)
}
#X.random <- X[sample(1:n.s,n.r.pt),]
#X.random <- jitter(X.random)
#X.range <- apply(X,2,range)
#n.random.pt <- 100
#X.random <- matrix(0,n.random.pt,0)
#for(i in 1:n.x){
#	X.random <- cbind(X.random,runif(n.random.pt,min=X.range[1,i],max=X.range[2,i]))
#}

# 正解生起確率を計算する
true.pr <- apply(X.random[,1:n.x.a],1,my.calc.p,k1,k2,k3)
# 推定生起確率を計算する
#est.pr.var <- my.dirichlet.est(X.random,est.out)
#est.pr <- est.pr.var$Est[,1]
#est.var <- est.pr.var$Var[,1]
# 推定値を評価する

#sd <- (true.pr-est.pr)/sqrt(est.var)

#plot(c(true.pr),c(sd))
#length(which(abs(sd)<=2))/length(sd)

#quantile(abs(sd),c(0.5,0.9,0.95,0.99))

#ord <- order(true.pr)

matplot(cbind(c(true.pr),c(est.pr),c(est.pr+2*sqrt(est.var)),c(est.pr-2*sqrt(est.var)))[ord,],type="p",cex=0.1)

# conf評価
conf.out <- my.Dirichlet.knn.conf(X.random,X,Y,Y.level,r.n=100,q=0.75,n.mc=100,est=FALSE,self=TRUE)

conf.out$K
conf.ok <- rep(0,conf.out$n.cnt)
est.pr <- est.var <- est.L <- est.U <- conf.ok 
#par(ask=TRUE)
for(i in 1:conf.out$n.cnt){
	plot(ps,dbeta(ps,conf.out$D.level.ord[i,conf.out$K[i],1]+1,conf.out$D.level.ord[i,conf.out$K[i],2]+1))
	abline(v=true.pr[i],col=2)
	conf.ok[i] <- pbeta(true.pr[i],conf.out$D.level.ord[i,conf.out$K[i],1]+1,conf.out$D.level.ord[i,conf.out$K[i],2]+1)
	est.pr[i] <- (conf.out$D.level.ord[i,conf.out$K[i],1]+1)/(conf.out$D.level.ord[i,conf.out$K[i],1]+1+conf.out$D.level.ord[i,conf.out$K[i],2]+1)
	est.var[i] <- (conf.out$D.level.ord[i,conf.out$K[i],1]+1)*(conf.out$D.level.ord[i,conf.out$K[i],2]+1)/((conf.out$D.level.ord[i,conf.out$K[i],1]+1)+(conf.out$D.level.ord[i,conf.out$K[i],2]+1))^2*((conf.out$D.level.ord[i,conf.out$K[i],1]+1)+(conf.out$D.level.ord[i,conf.out$K[i],1]+1)+1)
	est.L[i] <- qbeta(0.025,conf.out$D.level.ord[i,conf.out$K[i],1]+1,conf.out$D.level.ord[i,conf.out$K[i],2]+1)
	est.U[i] <- qbeta(0.925,conf.out$D.level.ord[i,conf.out$K[i],1]+1,conf.out$D.level.ord[i,conf.out$K[i],2]+1)

}
ord <- order(true.pr)
matplot(cbind(true.pr[ord],est.pr[ord],est.L[ord],est.U[ord]),type="l")

par(ask=FALSE)

plot(true.pr,est.pr)
plot(true.pr,est.pr,xlim=c(0,1),ylim=c(0,1))
hist(conf.ok)
plot(sort(conf.ok))
length(which(conf.ok<0.025))
length(which(conf.ok>0.975))
length(conf.ok)

par(mfcol=c(3,3))
plot(X,col=rgb(P,0,1-P),pch=20,main="データ点の真の確率")
plot(X,col=Y,pch=20,main="データ点の生起0/1")
plot(X.random,col=rgb(true.pr,0,1-true.pr),pch=20,main="知りたい点の真の生起確率")
plot(X.random,col=rgb(est.pr,0,1-est.pr),pch=20,main="知りたい点の推定生起確率")
plot(true.pr,est.pr,main="真の生起確率と推定確率",xlim=c(0,1),ylim=c(0,1))
tmp <- abs((true.pr-est.pr))
plot(X.random,cex=tmp*5,pch=20,main="真の確率と推定確率の差")

#plot(X.random,col=rgb(est.pr,(est.U-est.L),1-est.pr),pch=20)
ord <- order(true.pr)
matplot(cbind(true.pr[ord],est.pr[ord],est.L[ord],est.U[ord]),type="l",main="真の確率と推定確率とその推定信頼区間１")
ord <- order(est.pr)
matplot(cbind(true.pr[ord],est.pr[ord],est.L[ord],est.U[ord]),type="l",main="真の確率と推定確率とその推定信頼区間１")

par(mfcol=c(1,1))

logo

図書室

日記全体の目次
(記事一覧)

疾患遺伝子解析を学ぶためのカタログ(H18 3月改訂)の目次

 教科書

京都大学大学院医学研究科集中講義 全８限Part I(2007)
京都大学大学院医学研究科集中講義 全８限Part II(2008)
京都大学大学院医学研究科統計遺伝学講義 全８限(2009)

以下の関連サイトへのリンク道案内ページ
遺伝子多型解析Wiki
東大医科研HGCゲノム機能解析分野(旧)
理研リウマチチーム解析資料ページ


疾患遺伝子解析を学ぶためのカタログ(H18 3月改訂)
遺伝統計学のための統計学基礎 全９限
Haploviewを用いた連鎖不平衡マッピング解析実習 全９限
Arlequinを用いた集団遺伝学実習 全７限
一から始める遺伝子多型・形質関連解析
多数のGWAS/GWASメタアナリシスとその周辺解析用のアプリケーションへのリンク

複合遺伝性疾患とは

----------News----------
20万アクセス
2008/07/20
40万アクセス
2010/09/12


-------------------------