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

カテゴリカルな決断

決断 Decision theory カテゴリ ベータ分布 ディリクレ分布 単体
  • 昨日の記事で２つの選択肢の間の決断のことを書いた
  • そこでは、選択肢それぞれの帰結には「成功と失敗」の２つがあり、過去の記録に基づいて、その選択肢の成功の多寡に関する推測をして、より成功しやすそうな選択肢を確率的に選ぶことを書いた
  • 過去の記録に基づいて、成功・失敗の比率の推定にあたっては、ベータ分布を用いる方法を書いた
  • 一般化してみよう
  • 選択肢の数:2 →k
  • 個々の選択肢を選んだ帰結のパターン数:すべて2 → \mathbf{s} = (s_1,s_2,...,s_k)
    • ベータ分布からディリクレ分布へ
  • 選択肢の評価基準:「より成功しやすそうな」→ 「帰結パターンの生起確率ベクトルはs_i-1次元の正単体\Delta_i空間上の点であるが、それらに比較ルール\mathbf{\succeq}=\{\succeq_{\mathbf{v}_i,\mathbf{v}_j}\};i,j =1,2,...,kを入れる(推移性とか完全性は満足していなくても何とかすることにしよう。ここではk選択肢間の比較をk(k-1)/2ペアの比較にすることとして書いているが、それも、さらに一般化することは可能…ただし、生物学的にそのようにしているのかどうかは、今のところ全く検討していない)
    • ただし、帰結パターンの生起確率ベクトルをこう書くことにする
      • \mathbf{V} = (\mathbf{v}_1,\mathbf{v}_2,...,\mathbf{v}_k)
      • \mathbf{v}_i = (v_{i,1},v_{i,2},...,v_{i,s_i});i=1,2,...,k
  • 話が少しややこしくなってきているので次の二つに分割して進めよう
    • ベータ分布と一般化したディリクレ分布のこと
    • 複数の正単体空間の比較のこと

多項、ディリクレ分布：カテゴリカルな決断

決断 Decision theory カテゴリ ベータ分布 ディリクレ分布 単体
  • ベータ分布
    • \frac{1}{B(a_1,a_2)}x_1^{a_1-1}x_2^{a_2-1}
      • B(a_1,a_2) = \frac{\Gamma(a_1)\Gamma(a_2)}{\Gamma(a_1+a_2)}
      • 特にa_1,a_2が整数であれば、
      • B(a_1,a_2) = \frac{(a_1-1)!(a_2-1)!}{(a-1+a_2-1)!}
      • = \frac{(a_1-1)!(a_2-1)!}{(a_1-1+a_2-1)!(a_1+a_2-1)}
      • =\frac{1}{_{(a_1-1)+(a_2-1)}C_{a_1-1}} \times \frac{1}{a_1+a_2-1}
      • B(a_1+1,a_2+1) = \frac{1}{_(a_1+a_2)C_{a_1}}\times \frac{1}{a_1+a_2+1}
> x_1 <- 4
> x_2 <- 3
> beta(x_1,x_2)
[1] 0.01666667
> 1/choose(x_1-1+x_2-1,x_1-1) * 1/(x_1+x_2-1)
[1] 0.01666667
> beta(x_1+1,x_2+1)
[1] 0.003571429
> 1/choose(x_1+x_2,x_1) * 1/(x_1+x_2+1)
[1] 0.003571429
  • ディリクレ分布
    • \frac{1}{B(\mathbf{a})}\prod_{i=1}^K x_i^{a_i-1}=\frac{1}{B(\mathbf{a})} \mathbf{x}^{\mathbf{a}-1}
      • \mathbf{a} = (a_1,...,a_K),\mathbf{x}=(x_1,...,x_K),\sum_{i=1}^K x_i = 1
      • B(\mathbf{a}) = \frac{\prod_{i=1}^K\Gamma(a_i)}{\Gamma(\sum_{i=1}^K a_i)}
> library(MCMCpack)
> K <- 2
> a <- c(5,10)
> x <- runif(K)
> x <- x/sum(x)
> dbeta(x[1],a[1],a[2])
[1] 0.7927549
> ddirichlet(x,a)
[1] 0.7927549
> K <- 4
> a <- sample(10:20,K)
> x <- runif(K)
> x <- x/sum(x)
> ddirichlet(x,a)
[1] 8.641452e-09

正単体空間の比較・対応付け、カテゴリカルな決断

決断 Decision theory カテゴリ ベータ分布 ディリクレ分布 単体
  • 簡単のためにk=2と限定して考えてみる
  • s_1 = s_2 = 2の場合は、どんなことをしているのかがわかっているつもりに慣れるので、それに限局して考えてみる
    • それぞれの正単体は\Delta_1,\Delta_2はいずれも１次元空間の長さ１の線分である
    • この２つの正単体の直積空間は単位正方形である
    • ２つの選択肢に、同じ帰結「成功・失敗」があり、それを用いて１次元空間にスカラーが与えられるので
      • ここで２つの選択肢から比較が可能なスカラーが得られているが、ここを一般化するときには、そのような好条件に限定しないような枠組みが必要であることを注記する
    • 第１の選択肢の正単体をx軸方向の長さ１の線分、第２の選択肢の正単体をy軸方向の長さ１の線分、２つの選択肢のスカラーをz軸に取れば
    • 第１の選択肢に関しては、(x,y,z)3次元単位立方体において、(0,0,0),(1,0,1),(0,1,0),(1,1,1)を通る平面が
    • 第２の選択肢に関しては、同様に(0,0,0),(0,1,1),(1,0,0),(1,1,1)を通る平面が得られ、
    • この２平面の交線(0,0,0)と(1,1,1)とを結ぶ直線が得られる
    • この交線をxy平面に射影すると、この交線のx軸側とy軸側とは、それぞれ、第１選択肢優勢、第２選択肢優勢の領域を表している
    • xy平面には、第１選択肢、第２選択肢それぞれの生起確率の積(第１，２選択同士は独立なので)である生起確率が対応しているから、第１、２選択肢優性の確率が求められる
  • 以上を、kの値を任意としs_iの値も任意とできるように一般化したい
    • ３つに問題が分けられそうだ
      • １つ目。ある一つのオプションに関してs_i-1次元正単体空間に推移性や完全性のある順序がある場合とない場合があるだろうから、それらを統一的に取り扱う方法(ルール？）
      • ２つ目。ある二つのオプションに関してs_i-1,s_j-1次元正単体空間のそれぞれに、何かしらの値？値のセット？比較ルール？があったとして、２つの異なるルールの比較の方法
      • ３つ目。ある任意の数のオプションに関して、ペアワイズな比較ではない、比較を取り入れる方法

帰結の正単体、カテゴリカルな決断

決断 Decision theory カテゴリ ベータ分布 ディリクレ分布 単体
  • もう少し考えよう
  • ２つの選択肢があって、それぞれに３つの帰結がある。その３つの帰結のうち２つは共通であるとすれば、帰結は全部で４つある
  • この帰結４つを正四面体の４点に置く
  • こうすれば、ある選択肢の結果はこの正四面体のある面上の分布になる
  • ２つの選択肢の良しあしは、正四面体の隣接する２面に関する評価関数と、２面の面上にある確率密度分布(これはディリクレ分布)となる

奇病の治療法の確立、正単体空間の比較・対応付け、カテゴリカルな決断

決断 Decision theory カテゴリ ベータ分布 ディリクレ分布 単体
  • 奇病が発生した。この奇病にかかると、音楽と美術とスポーツに関する能力が破壊されるという。音楽は演奏できなくなるし鑑賞もできなくなる。美術は創作活動ができなくなるし鑑賞もできなくなる。スポーツは自らプレイすることができなくなるし、鑑賞もできなくなる。
  • 今、３つの治療法、M,A,Sが開発されたと言うが、効果は不明だという。この３つの治療法は一つしか選べない(併用すると効果は無になることが分かっている)とする。
  • また、治療は即刻開始しなければ、無効であり、副作用はないことにする。
  • M,A,Sの効果はそれぞれ、音楽系、美術系・スポーツ系にのみ効果があり、他の系には全く無効である
  • また、効果の出方は４通りある。
    • 「著効」
    • 「半効」
    • 「無効」
  • それぞれの効果をM1,M2,M3,A1,A2,A3,S1,S2,S3とする
  • 効果は立場によって優劣が異なる。例を挙げる。
    • ホタル界ではA1>A2=S1>S2>M1=M2=M3=A3=S3と評価された
    • ウグイス界ではM1=S1>S2>M2=A1>A2>M3=A3=S3と評価された
    • ちなみに、選択肢２つで、その帰結が「成功と失敗」の２つの場合には、X1=Y1>X2=Y2という優劣情報になっている。

# オプション数
k <- 3
# 各オプションの帰結の種類数
s <- sample(2:5,k,replace=TRUE)
# 各オプション３個にしてみる
s <- rep(3,k)
s.list <- list()
for(i in 1:k){
	s.list[[i]] <- 1:s[i]
}
# 各オプションごとに帰結の生起確率
true.prob <- list()
for(i in 1:k){
	true.prob[[i]] <- rdirichlet(1,rep(10,s[i]))
}
# ちょっと固定してみる
# オプション２とオプション３とが競り合うように
#true.prob[[1]] <- c(0.5,0.5)
#true.prob[[2]] <- c(0.45,0.55)
true.prob[[3]] <- c(0.5,0.1,0.4)
true.prob[[2]] <- c(0.0,0.5,0.5)
# 起こりうる帰結の組合せの網羅
cat.combin <- expand.grid(s.list)

# どの帰結を好むかの行列を作る
#preference <- sample(1:k,prod(s),replace=TRUE)
preference <- matrix(1,prod(s),k)
for(i in 1:prod(s)){
	tmp.s <- sample(0:(k-1),1)
	preference[i,sample(1:k,tmp.s)] <- 0
	preference[i,] <- preference[i,]/sum(preference[i,])
}
#preference <- matrix(c(1,1,0,1,1,0,1,1),ncol=2)
# ホタルのモデル
scores <- matrix(c(1,4,3,1,3,2,1,1,1),byrow=TRUE,3,3)
for(i in 1:length(preference[,1])){
	tmp <- c(scores[cat.combin[i,1],1],scores[cat.combin[i,2],2],scores[cat.combin[i,3],3])
	preference[i,] <- as.numeric(tmp==max(tmp))
}
# 和が１となるように標準化
for(i in 1:length(preference[,1])){
	preference[i,] <- preference[i,]/sum(preference[i,])
}
# 生起確率は全く不明なところからスタート
pre <- list()
for(i in 1:k){
	pre[[i]] <- sample(0:0,s[i],replace=TRUE)
}
# そのつど、確率的によさげなオプションを選んでは、その結果をhistoryに記録し
# 記録に応じて、次の選択を判断する
n.history <- 100
history <- matrix(0,n.history,sum(s))
# 結果からディリクレ乱数で「生起確率」を推定する
# 作成するディリクレ乱数の個数
n.iter <- 1000
#N.iter <- 1

for(h in 1:n.history){
	tmp <- c()
	for(j in 1:k){
			tmp <- c(tmp,pre[[j]])
	}
	history[h,] <- tmp
	chosen <- rep(0,k)
	#for(ii in 1:N.iter){
		this.time.chosen <- rep(0,k)
		rs <- list()
		for(j in 1:k){
			rs[[j]] <- rdirichlet(n.iter,pre[[j]]+1)
		}

		# 推定した生起確率とディリクレ乱数から
		# もっとも好みに沿った結果が得られるオプションを
		# 確率的に選ぶ
		for(i in 1:n.iter){
			tmp.rs <- list()
			for(j in 1:k){
				tmp.rs[[j]] <- rs[[j]][i,]
			}
			probs <- expand.grid(tmp.rs)
			tmp.prob <- apply(probs,1,prod)
			tmp2 <- preference * tmp.prob
			tmp3 <- apply(tmp2,2,sum)
			tmp.max <- which(tmp3 == max(tmp3))
			if(length(tmp.max)==1){
				choice <- tmp.max
			}else{
				choice <- sample(tmp.max,1)
			}
			#choice <- tmp.max[sample(seq(tmp.max),1)]
			this.time.chosen[choice] <- this.time.chosen[choice]+1
			#print(choice)
			#tmp <- matrix(NA,N.iter,k)
			#for(j in 1:k){
			#	tmp[,j] <- sample(1:s[j],replace=TRUE,prob=rs[[j]])
			#}
		}
		#tmp.max <- which(this.time.chosen==max(this.time.chosen))
		#if(length(tmp.max)==1){
		#	selected.chosen <- tmp.max
		#}else{
		#	selected.chosen <- sample(tmp.max,1)
		#}
		#selected.chosen <- tmp.max[sample(seq(tmp.max),1)]
		chosen[selected.chosen] <- chosen[selected.chosen]+1
	#}
	#tmp.max <- which(chosen==max(chosen))
	#if(length(tmp.max)==1){
	#	really.chosen <- tmp.max
	#}else{
	#	really.chosen <- sample(tmp.max,1)
	#}
	# 選んだオプションの結果は、隠された真の生起確率によって観察される
	really.chosen <- sample(1:k,1,prob=this.time.chosen/sum(this.time.chosen))
	#really.chosen <- tmp.max[sample(seq(tmp.max),1)]
	really.happen <- sample(1:s[really.chosen],1,prob=true.prob[[really.chosen]])
	pre[[really.chosen]][really.happen] <- pre[[really.chosen]][really.happen]+1

}

matplot(history,type="l")

決断関数の構造、カテゴリカルな決断

決断 Decision theory カテゴリ ベータ分布 ディリクレ分布 単体
  • 上の例では、複数の選択肢から発生するすべての帰結に順序を入れた。さらに、その順序に応じて、なるべく、「順序が前」の帰結が得られるように「単純な」処理を繰り返した
  • 「順序〜大小関係」のみを用いるというのは、結局のところ、不等号には重み１の差をつけ、等号には重みの差をつけずに、「線形な面」をつくっているのと(多分)同じ
  • そうすると、「治療選択における、治癒と死亡」ただし、「治癒にも副作用ありとなしがあって、死亡にも副作用ありなしがある」のような場合には、「副作用の有無によらず治癒」は「副作用の有無によらず死亡」の上位におきたい
  • すべての帰結を同一尺度に並べ、ただし、不等号に持たせる重みに軽重をつける、というのも一つのやり方だろう
  • また別のやり方としては、「順序〜大小関係」にヒエラルキーを入れることになるだろう
  • いずれにしろ、その「関数・重みづけ係数」がわかっていれば、簡単なわけで、それがわからないながら、前進したい、と、それをする方便を入れないといけない
  • そんな方向としては、モンテカルロ比較をしながら、「選択肢Aだと、こうなって、選択肢Bだとこうなる、と仮定したら、AとBとどっちをとりますか？」だけでなく(これだと順序を入れただけ)、どのくらい強く(主観的で酔い)Aをとりたいですか？」とインターラクティブにモンテカルロを回す、というのもありかもしれない。
  • ひとまず、２治療、副作用持ち込み編(順序だけで評価)の書き散らしソースを
# オプション数
k <- 2
# 各オプションの帰結の種類数
s <- sample(2:5,k,replace=TRUE)
# 各オプション３個にしてみる
s <- rep(4,k)
s.list <- list()
for(i in 1:k){
	s.list[[i]] <- 1:s[i]
}
# 各オプションごとに帰結の生起確率
true.prob <- list()
for(i in 1:k){
	true.prob[[i]] <- rdirichlet(1,rep(10,s[i]))
}
# ちょっと固定してみる
# オプション２とオプション３とが競り合うように
#true.prob[[1]] <- c(0.5,0.5)
#true.prob[[2]] <- c(0.45,0.55)
#true.prob[[3]] <- c(0.5,0.1,0.4)
#true.prob[[2]] <- c(0.0,0.5,0.5)
true.prob[[1]] <- c(0.0,0.6,0.1,0.3)
true.prob[[2]] <- c(0.4,0.15,0.15,0.3)
# 起こりうる帰結の組合せの網羅
cat.combin <- expand.grid(s.list)

# どの帰結を好むかの行列を作る
#preference <- sample(1:k,prod(s),replace=TRUE)
preference <- matrix(1,prod(s),k)
for(i in 1:prod(s)){
	tmp.s <- sample(0:(k-1),1)
	preference[i,sample(1:k,tmp.s)] <- 0
	preference[i,] <- preference[i,]/sum(preference[i,])
}
#preference <- matrix(c(1,1,0,1,1,0,1,1),ncol=2)
# ホタルのモデル
scores <- matrix(c(7,4,3,1,6,5,3,2),ncol=2)
for(i in 1:length(preference[,1])){
	tmp <- c(scores[cat.combin[i,1],1],scores[cat.combin[i,2],2])
	preference[i,] <- as.numeric(tmp==max(tmp))
}
# 和が１となるように標準化
for(i in 1:length(preference[,1])){
	preference[i,] <- preference[i,]/sum(preference[i,])
}
# 生起確率は全く不明なところからスタート
pre <- list()
for(i in 1:k){
	pre[[i]] <- sample(0:0,s[i],replace=TRUE)
}
# そのつど、確率的によさげなオプションを選んでは、その結果をhistoryに記録し
# 記録に応じて、次の選択を判断する
n.history <- 1000
history <- matrix(0,n.history,sum(s))
# 結果からディリクレ乱数で「生起確率」を推定する
# 作成するディリクレ乱数の個数
n.iter <- 100
#N.iter <- 1

for(h in 1:n.history){
	tmp <- c()
	for(j in 1:k){
			tmp <- c(tmp,pre[[j]])
	}
	history[h,] <- tmp
	chosen <- rep(0,k)
	#for(ii in 1:N.iter){
		this.time.chosen <- rep(0,k)
		rs <- list()
		for(j in 1:k){
			rs[[j]] <- rdirichlet(n.iter,pre[[j]]+1)
		}

		# 推定した生起確率とディリクレ乱数から
		# もっとも好みに沿った結果が得られるオプションを
		# 確率的に選ぶ
		for(i in 1:n.iter){
			tmp.rs <- list()
			for(j in 1:k){
				tmp.rs[[j]] <- rs[[j]][i,]
			}
			probs <- expand.grid(tmp.rs)
			tmp.prob <- apply(probs,1,prod)
			tmp2 <- preference * tmp.prob
			tmp3 <- apply(tmp2,2,sum)
			tmp.max <- which(tmp3 == max(tmp3))
			if(length(tmp.max)==1){
				choice <- tmp.max
			}else{
				choice <- sample(tmp.max,1)
			}
			#choice <- tmp.max[sample(seq(tmp.max),1)]
			this.time.chosen[choice] <- this.time.chosen[choice]+1
			#print(choice)
			#tmp <- matrix(NA,N.iter,k)
			#for(j in 1:k){
			#	tmp[,j] <- sample(1:s[j],replace=TRUE,prob=rs[[j]])
			#}
		}
		tmp.max <- which(this.time.chosen==max(this.time.chosen))
		if(length(tmp.max)==1){
			selected.chosen <- tmp.max
		}else{
			selected.chosen <- sample(tmp.max,1)
		}
		selected.chosen <- tmp.max[sample(seq(tmp.max),1)]
		chosen[selected.chosen] <- chosen[selected.chosen]+1
	#}
	#tmp.max <- which(chosen==max(chosen))
	#if(length(tmp.max)==1){
	#	really.chosen <- tmp.max
	#}else{
	#	really.chosen <- sample(tmp.max,1)
	#}
	# 選んだオプションの結果は、隠された真の生起確率によって観察される
	really.chosen <- sample(1:k,1,prob=this.time.chosen/sum(this.time.chosen))
	#really.chosen <- tmp.max[sample(seq(tmp.max),1)]
	really.happen <- sample(1:s[really.chosen],1,prob=true.prob[[really.chosen]])
	pre[[really.chosen]][really.happen] <- pre[[really.chosen]][really.happen]+1

}

matplot(history,type="l")

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


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