多人数一括DNAプロファイリングのための 確率計算法に関する考察

Slides:



Advertisements
Similar presentations
HBSP モデル上での 行列積を求めるアルゴリ ム 情報論理工学 吉岡健太.
Advertisements

寺尾 敦 青山学院大学社会情報学部 Fisher の直接確率法 寺尾 敦 青山学院大学社会情報学部
看護学部 中澤 港 統計学第5回 看護学部 中澤 港
点対応の外れ値除去の最適化によるカメラの動的校正手法の精度向上
多人数一括DNAプロファイリングのための 確率計算法に関する考察
確率・統計Ⅰ 第12回 統計学の基礎1 ここです! 確率論とは 確率変数、確率分布 確率変数の独立性 / 確率変数の平均
多変量解析 -重回帰分析- 発表者:時田 陽一 発表日:11月20日.
プログラミング論 I 行列の演算
Bassモデルにおける 最尤法を用いたパラメータ推定
アルゴリズムイントロダクション第5章( ) 確率論的解析
法数学勉強会(京大法医学講座) 2012/02/18 京都大学 統計遺伝学 山田
Bias2 - Variance - Noise 分解
社会心理学のStudy -集団を媒介とする適応- (仮)
ベイズ的ロジスティックモデル に関する研究
IT入門B2 ー 連立一次方程式 ー.
3章 Analysing averages and frequencies (前半 p )
確率・統計輪講資料 6-5 適合度と独立性の検定 6-6 最小2乗法と相関係数の推定・検定 M1 西澤.
最尤推定によるロジスティック回帰 対数尤度関数の最大化.
精密工学科プログラミング基礎Ⅱ 第3回資料 今回の授業で習得してほしいこと: 2次元配列の使い方 (前回の1次元配列の復習もします.)
法数学勉強会 2011/11/26 京都大学大学院医学研究科 統計遺伝学分野 山田 亮
高次元データの解析 -平均ベクトルに関する検定統計量の 漸近分布に対する共分散構造の影響-
k 個のミスマッチを許した点集合マッチング・アルゴリズム
7. 音声の認識:高度な音響モデル 7.1 実際の音響モデル 7.2 識別的学習 7.3 深層学習.
遺伝的アルゴリズムへの 統計力学的アプローチ 大阪大学 大学院理学研究科 鈴木譲 CISJ2005 於早稲田大学理工学部
決定木とランダムフォレスト 和田 俊和.
奈良女子大集中講義 バイオインフォマティクス (9) 相互作用推定
確率伝搬法と量子系の平均場理論 田中和之 東北大学大学院情報科学研究科
細胞の形と変形のための データ駆動型解析手法
独立成分分析 5 アルゴリズムの安定性と効率 2007/10/24   名雪 勲.
訓練データとテストデータが 異なる分布に従う場合の学習
大規模なこと Large scale.
ゲノム科学概論 ~ゲノム科学における統計学の役割~ (遺伝統計学)
法数学勉強会 2016年4月会 京都大学(医)統計遺伝学分野 山田 亮
法数学における ベイジアンネットワーク(2) ~成書で学ぶ~
藤田保健衛生大学医学部 公衆衛生学 柿崎 真沙子
法数学勉強会 2016/06/15 京都大学(医)統計遺伝学分野 山田 亮
法数学勉強会(京大法医学講座) 2012/02/18 京都大学 統計遺伝学 山田
東日本大震災におけるご遺体身元確認と行方不明家族捜索のためのDNA鑑定
2011/05/28 京都大学大学院 附属ゲノム医学センター統計遺伝学分野 山田 亮
5 Recursions 朴大地.
第24回日本法科学技術学会 2018/11/08 京都大学 医学研究科 統計遺伝学分野 山田 亮
確率と統計2009 第12日目(A).
尤度の比較と仮説検定とを比較する ~P値のことなど~
第3章 線形回帰モデル 修士1年 山田 孝太郎.
法数学のための 機械学習の基礎 京大(医) 統計遺伝学分野 山田 亮 2017/04/15.
「アルゴリズムとプログラム」 結果を統計的に正しく判断 三学期 第7回 袖高の生徒ってどうよ調査(3)
決断のための分布合算 京大(医)統計遺伝学分野 山田 亮.
法医学会 2013年6月26日 京都大学(医)統計遺伝学 山田 亮
最尤推定・最尤法 明治大学 理工学部 応用化学科 データ化学工学研究室 金子 弘昌.
親子鑑定に見る尤度比を 角度を変えて眺めてみる
法数学勉強会 2015/09/26 京都大学統計遺伝学分野 山田 亮
法数学勉強会 2015/09/26 京都大学統計遺伝学分野 山田 亮
東北大 情報科学 田中和之,吉池紀子 山口大 工 庄野逸 理化学研究所 岡田真人
HMM音声合成における 変分ベイズ法に基づく線形回帰
DNA鑑定を理解するために必要な数学の学び方
藤田保健衛生大学医学部 公衆衛生学 柿崎 真沙子
音響伝達特性モデルを用いた シングルチャネル音源位置推定の検討 2-P-34 高島遼一,住田雄司,滝口哲也,有木康雄 (神戸大) 研究の背景
精密工学科プログラミング基礎 第7回資料 (11/27実施)
行列 一次変換,とくに直交変換.
Chapter5 Systems of Distinct Representatives
森 裕一(岡山理科大学) 山本義郎(岡山大学自然科学研究科) 渡谷真吾,尾高好政(倉敷芸術科学大学) 垂水共之,田中 豊(岡山大学)
法数学における ベイジアンネットワーク(2) ~成書で学ぶ~
割り当て問題(assignment problem)
精密工学科プログラミング基礎Ⅱ 第2回資料 今回の授業で習得してほしいこと: 配列の使い方 (今回は1次元,次回は2次元をやります.)
Q状態イジング模型を用いた多値画像修復における 周辺尤度最大化によるハイパパラメータ推定
ガウシアングラフィカルモデルにおける一般化された確率伝搬法
混合試料の構成人数 Nuisance パラメタ
計算技術研究会 第5回 C言語勉強会 関数(function)を使う
混合ガウスモデル Gaussian Mixture Model GMM
プログラミング論 バイナリーサーチ 1.
Presentation transcript:

多人数一括DNAプロファイリングのための 確率計算法に関する考察 法数学勉強会 2011/09/10 京都大学大学院医学研究科 統計遺伝学分野 奈良原舞子 山田 亮

状況 大災害が発生 多数の行方不明者 多数の身元不明遺体 外見や所持品などから身元が特定された遺体はすでに遺族に返還されている。 残っている遺体の手がかりは、主に遺伝情報

使えるデータ 遺体 ジェノタイプ 発見時の状況など 行方不明者 家系情報 家族のジェノタイプ

個人の鑑定 簡単に描くと 行方不明者 Missing 身元不明者 found Body m1 b1 集団の誰か

集団の鑑定 .. .. .. .. .. .. 行方不明者 Missing 身元不明者 found Body m1 b1 m2 b2 m3 mN bN

集団の鑑定 N人をN体に割り付ける N! = N×(N-1) ×(N-2) ×... ×2×1 通り (m1,m2,...,mN)=(b(s1),b(s2),...,b(sN)) 割り付け方:Si=(s1,s2,...,sN) がN!通り

Si=(i1,i2,...,iN)という割り付け m1 = b(s1) m2 = b(s2) m1 = b(s2) ... ... mN = b(sN)

Si=(i1,i2,...,iN)という割り付けを 観察する確率は? P(m1=b(s1)) P(m2=b(s2)) P(m3=b(s3)) × P(m2=b(s2)) × P(m3=b(s3)) m1 = b(s1) m2 = b(s2) m1 = b(s2) × × × P(mN=b(sN)) ... ... mN = b(sN)

N!通りの確率 P(S1),P(S2),...,...,...,...,P(SN!) 最も大きなP(Si)となるSiは最尤推定割り付けがある

N!通りの確率 P(S1),P(S2),...,...,...,...,P(SN!) 計算する...? 1!=1 2!=2 3!=6 4!=24 5!=120 6!=720 7!=5040 8!=40320 9!=362880 10!=3,628,800 3百万 11!=39,916,800 12!=479,001,600 4.8億 15!=1.3 x 1012 20!=2.4 x 1018

多すぎてN!通りを計算できない N!通りを計算しないで、最尤割り付けがわかる? 重みづけ最適化・重みづけマッチング ハンガリアン・アルゴリズムなど

最尤推定割り付けがわかればそれが「答え」なのか? P(Si) と P(Sj)とが第一位、第二位だとする P(Si)とP(Sj)とが等しかったら… P(Si)とP(Sj)とがほぼ等しかったら… P(Si)とP(Sj)とが数倍の違いしかなかったら…

第1,2,…n位割り付けがわかればそれが「答え」なのか? 「僅差」の割り付けがあったら、結局、どうしたらよいのかわからない 尤度が高い割り付けパターンを探すだけでは、解決しないかも

ある家族 「我が家の行方不明者mは、遺体b1,...,bNのうちのどれか1体だと言えますか?それとも、言いかねますか?」 ある遺体を保管しているところ 「この遺体bは、探されている行方不明者m1,...,mNのだれか1人だと言えますか?それとも、言いかねますか?」

N=3で考える 3 体の遺体と3人の不明者を割り付ける場合の数 6通りの割り付け N=3 の場合 M1 M2 M3 仮説1 B1 B2 B3 仮説2 仮説3 仮説4 仮説5 仮説6 確率行列 B1 B2 B3 M1 P(M1=B1) P(M1=B2) M2 P(M2=B1) P(M2=B2) P(M2=B3) M3 P(M3=B1) P(M3=B2) P(M3=B3) 各仮説の尤度: 3つのペア全てでジェノタイプ が一致する確率

仮説2 6通りの割り付け N=3 の場合 M1 M2 M3 仮説1 B1 B2 B3 仮説2 仮説3 仮説4 仮説5 仮説6 確率行列 B1 P(M1=B1) P(M1=B2) P(M1=B3) M2 P(M2=B1) P(M2=B2) P(M2=B3) M3 P(M3=B1) P(M3=B2) P(M3=B3) この仮説の尤度=P(M1=B1)xP(M2=B3)xP(M3=B2)

M1=B1,それ以外の割り付けはなんでもあり 仮説1+2 6通りの割り付け N=3 の場合 M1 M2 M3 仮説1 B1 B2 B3 仮説2 仮説3 仮説4 仮説5 仮説6 確率行列 B1 B2 B3 M1 P(M1=B1) P(M1=B2) P(M1=B3) M2 P(M2=B1) P(M2=B2) P(M2=B3) M3 P(M3=B1) P(M3=B2) P(M3=B3) この仮説の尤度=P(M1=B1)xP(M2=B2)xP(M3=B3)

M1=B1,それ以外の割り付けはなんでもあり 仮説1+2 6通りの割り付け N=3 の場合 M1 M2 M3 仮説1 B1 B2 B3 仮説2 仮説3 仮説4 仮説5 仮説6 確率行列 B1 B2 B3 M1 P(M1=B1) P(M1=B2) P(M1=B3) M2 P(M2=B1) P(M2=B2) P(M2=B3) M3 P(M3=B1) P(M3=B2) P(M3=B3) この仮説の尤度=P(M1=B1)xP(M2=B3)xP(M3=B2)

行列式(Determinant) 割り付けの場合ごとに掛け算をする「加える」要素と「引く」要素がある Wikipedia

パーマネント 割り付けの場合ごとに掛け算をする 全部を「加える」

パーマネント 行列式の計算は簡単で正確 この6通りの確率の和が3次正方行列のパーマネント パーマネントを求めるためのいくつかの方法がある 近似的 近似法を使うことでだいぶ速く計算できる。

パーマネントの計算方法 今日は、割愛 気になる方は Wikipedia http://en.wikipedia.org/wiki/Permanent から情報の入手は可能です library(sets) library(gregmisc) # Ryser PermanentRyser<-function(A){ n<-length(A[,1]) p<-0 s<-as.set(1:n) for(i in 1:n){ scomb<-set_combn(s,i) #print(length(scomb)) c<-(-1)^i sumSameLen<-0 for(j in scomb){ #print(j) tmpprod<-1 for(k in 1:n){ tmpsum<-0 #print(k) for(l in j){ #print(A[k,l]) tmpsum<-tmpsum+A[k,l] } tmpprod<-tmpprod*tmpsum sumSameLen<-sumSameLen+tmpprod p<-p+c*sumSameLen p*(-1)^n # calculate based om difinition of permutation PermanentBeta<-function(A){ perms<-permutations(n,n) for(i in 1:length(perms[,1])){ t<-1 for(j in 1:n){ t<-t*A[j,perms[i,j]] p<-p+t p # Sinkhorn Balancing Scale sinkhorn<-function(A,epsilon=0.00001){ B<-A scaleX<-scaleY<-rep(1,length(A[,1])) row_sum<-apply(A,1,sum) while(max(abs(t(row_sum)-1))>epsilon){ B<-B/row_sum scaleX<-scaleX/row_sum col_sum<-apply(B,2,sum) scaleY<-scaleY/col_sum B<-t(t(B)/col_sum) row_sum <- apply(B,1,sum) list(mat=B,scaleX=scaleX,scaleY=scaleY) # Huber/Law factor HuberLaw<-function(x){ x<-x x[which(x==0)]<-x[which(x==0)]+0.5 hl<-x mores<-which(hl>1) unMores<-which(hl<=1) hl[mores]<-hl[mores]+0.5*log(hl[mores])+exp(1)-1 hl[unMores]<-hl[unMores]*(exp(1)-1)+1 hl/(exp(1)) Minc<-function(x){ toosmall<-which(x<10^(-10)) x2<-x x2[toosmall]<-10^(-10) exp(lgamma(x2+1)*(1/x2)) # Approximation permanentEstimation<-function(A,iter=10000,epsilon=10^(-4),factor="hl"){ n<-length(A[1,]) shOut<-sinkhorn(A,epsilon) #print(shOut) B<-shOut[[1]] x<-shOut[[2]] y<-shOut[[3]] row_scale<-1/apply(B,1,max) C<-B*row_scale save_C<-C number_of_successes<-0 for(i in 1:iter){ column<-1 C<-save_C row_sums<-apply(C,1,sum) while(column<=n){ #print(column) if(factor=="hl"){ hl<-prod(HuberLaw(row_sums)) }else if(factor=="minc"){ hl<-prod(Minc(row_sums)) tmphl2<-HuberLaw(row_sums-C[,column]) tmphl2<-Minc(row_sums-C[,column]) hl2<-prod(tmphl2) #print(hl) #print(hl2) #print(tmphl2) row_probs<-C[,column]/tmphl2*hl2/hl #print(row_probs) # 重みづけサンプリング cumsumprob<-cumsum(row_probs) rrand<-runif(1) row_pick<-sum(cumsumprob<rrand)+1 #print(row_pick) if(row_pick==(n+1)){ column<-n+2 }else{ row_sums<-row_sums-C[,column] C[row_pick,]<-rep(0,n) column<-column+1 row_sums[row_pick]<-0 if(column==(n+1)){ number_of_successes=number_of_successes+1 hl_C<-prod(HuberLaw(row_sums)) hl_C<-prod(Minc(row_sums)) per_estimate<-hl_C*number_of_successes/iter per_estimate<-per_estimate/prod(row_scale)/prod(x)/prod(y) per_estimate ### MAIN FUNCTION ### # Calculate permanent CalcPermanent <- function(M,approx=T,epsilon=0.0001,iter=1000,factor="hl",beta=F, maxnfor.exact=16){ dimM <- dim(M) if(length(dimM)!=2) stop("M must have two dimensions") if(dimM[1]!=dimM[2]) stop("matrix must be a square matrix") n <- nrow(M) if(!approx){ if(n > maxnfor.exact) stop("n is too large for exact method. Use approx=T") else if(beta && n < 10) p <- PermanentBeta(M) else if(!beta || n >=10) p <- PermanentRyser }else if(approx){ if(any(M< 0)) stop("all elements of M must be non-negative values") p <- permanentEstimation(M,iter,epsilon,factor)

 近似法は速い 計算にかかる時間 (sec) 近似法 正確法 (N: 行列サイズ)

Fluctuation of estimation もとがNxN行列なら2xN 個の L(All) が出る 各列和・各行和 そのL(All) は推定値なので、誤差がある 最大値と最小値の差は0.102 注: scale 調整後

パーマネント計算 近似法の精度 「真のパーマネント」を大きな行列で計算するのは非現実的なので、「真の精度」を評価するのは難しいのですが パーマネント計算 近似法の精度 「真のパーマネント」を大きな行列で計算するのは非現実的なので、「真の精度」を評価するのは難しいのですが Approximating the Permanent with Belief Propagation, by Bert Huang and Tony Jebara @ http://www.cs.columbia.edu/~bert/permanentTR.pdf これは別なパーマネント近似法ですが…

尤度のNxN行列 「mi=bj、あとは何でもあり」に対応する(N-1)!仮説の確率を合算する この行列の各行の和は、どの行も等しい 各行の和は以下の和 「mi=b1、あとは何でもあり」 「mi=b2、あとは何でもあり」 … 「mi=bN、あとは何でもあり」 これは「miも何でもあり、他も何でもあり」だから

尤度のNxN行列 「mi=bj、あとは何でもあり」に対応する(N-1)!仮説の確率を合算する この行列の各列の和もやはり等しい 各列の和は、各行の和とも等しい

行を列に入れ替えても同じこと 家族が知りたいことにも、遺体保管者が知りたいことにも、答えられる ある家族 「我が家の行方不明者mは、遺体b1,...,bNのうちのどれか1体だと言えますか?それとも、言いかねますか?」 列 ある遺体を保管しているところ 「この遺体bは、探されている行方不明者m1,...,mNのだれか1人だと言えますか?それとも、言いかねますか?」

尤度割合のNxN行列 尤度のNxN行列の各列の和、各行の和はすべて等しいので、その値で、尤度のNxN行列のすべての成分を割ってやる 各行、各列の和は、すべて1

2つのNxN行列 確率行列 尤度割合行列 P(mi = bj) の行列 「m1=b2、あとは何でもあり」に対応する(N-1)!仮説の確率を合算

NxN の確率行列 割り付けの計算のために正方行列がほしい。 全ての遺体と行方不明者が1対1対応すると仮定 仮定できなければ、足りない分を一般集団で補う B1 B2 B3 ・・・・ BN M1 P(M1=B1) P(M1=B2) M2 P(M2=B1) P(M2=B2) P(M2=B3) M3 P(M3=B1) P(M3=B2) P(M3=B3) : MN P(MN=B1) P(MN=B2) P(MN=B3)

NxN尤度割合行列 L(mi=bj) / L(ALL) の NxN 行列ができる 基準 v を満たしたペアの割り付けが決定する ・・・・ BN M1 0.9998 0.0000 M2 0.9953 M3 0.7 M4 0.0030 M5 0.0001 0.3 : MN

行を列に入れ替えても同じこと 家族が知りたいことにも、遺体保管者が知りたいことにも、答えられる ある家族 「我が家の行方不明者mは、遺体b1,...,bNのうちのどれか1体だと言えますか?それとも、言いかねますか?」 列 ある遺体を保管しているところ 「この遺体bは、探されている行方不明者m1,...,mNのだれか1人だと言えますか?それとも、言いかねますか?」

遺体引き取り と 遺体引き渡し NxN尤度割合行列のセルの値を使って、「遺体引き取り」「遺体引き渡し」の判断ができるだろう 閾値は…

処理フロー Missings Bodies 血縁関係情報 遺体のDNA型 血縁者のDNA型 確率行列 パーマネント計算 尤度割合行列 引き取り・引き渡し判定

simulation data 想定した全体 手元のデータ 全ての行方不明者:104人 100家系 + 重複4家系 100家系 + 重複4家系 重複:行方不明者が複数いる家系 手元のデータ 100家系のうち68家系(incl. 重複3家系) 本当は全部使う予定だったが時間の関係で途中まで計算したところで割り付けをしたので家系も足りない状況になった。 104人の不明者のうち、54人をランダムに選択  家系も遺体も足りない状況

Simulation 結果

正解のペア全:30ペア v > 0.999 : 27 ペア ☆ v > 0.9 : 1 ペア ! 家系として0.999 を満たす 全て正解 v > 0.9 : 1 ペア ! 正解 家系として0.999 を満たす 1家系ー2遺体 が該当 0.999で感度 27/30  ハズレなし

permanent計算後の確率(尤度比)が0でなかったペアの値をsort してplot したもの 正解ペア はずれペア 同一家系の2人 正解ペア

尤度比の自然対数 vs 1/(1-p) の自然対数 (p: 尤度割合)

可能性と課題 可能性 課題 「事前確率」を取り込める 「一般化」 申請のない行方不明者と発見されていない遺体のとりあつかい 性別・所持品等の情報を容易に取り込める 「一般化」 個人鑑定を同じ枠組みで考えることが(おそらく)可能 課題 申請のない行方不明者と発見されていない遺体のとりあつかい 1家系に複数の行方不明者がいる場合 非独立な確率・・・行列を使った計算が苦手