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

Presentation on theme: "多人数一括ＤＮＡプロファイリングのための 確率計算法に関する考察"— Presentation transcript:

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

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

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

4 個人の鑑定 と 集団の鑑定 と その違い

5 行方不明者が 身元不明者のＤＮＡ型を持つ確率
この確率が計算できる？

6 行方不明者が 身元不明者のＤＮＡ型を持つ確率
はい、できます！ 行方不明者 身元不明者 この確率が計算できる？

7 個人の鑑定 比較１ 尤度比 身元不明者　　　　　が 家系情報のない誰かである 身元不明者　　　　　が 行方不明者である

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

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

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

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

12 集団の鑑定 N人をN体に割り付ける .. .. .. .. .. .. 行方不明者 Missing 身元不明者 found Body m1
mN bN

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

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

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

16 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)

17 N!通りの確率 P(S1),P(S2),...,...,...,...,P(SN!)

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

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

20 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

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

22 N!通りを計算できない N!通りを計算しないで、最尤割り付けがわかる？ はい、わかります！

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

24 最尤推定割り付けがわかればそれが「答え」なのか？
P(Si) と P(Sj)とが第一位、第二位だとする P(Si)とP(Sj)とが等しかったら… P(Si)とP(Sj)とがほぼ等しかったら… P(Si)とP(Sj)とが数倍の違いしかなかったら… 尤度が最高の割り付けパターンを探すだけでは、解決しないかも

25 最尤推定割り付けがわかればそれが「答え」なのか？
やはりN!通りの計算が必要なのか？ 最尤割り付け、第２位割り付けだけがわかればよいのか 第１位でなくて、第２位もわかる？

26 最尤推定割り付けがわかればそれが「答え」なのか？
やはりN!通りの計算が必要なのか？ 最尤割り付け、第２位割り付けだけがわかればよいのか 第１位でなくて、第２位もわかる？ 少し大変だけれど、第 i 位割り付けは計算できる はい、わかります！

27 第１,２,…n位割り付けがわかればそれが「答え」なのか？
「僅差」の割り付けがあったら、結局、どうしたらよいのかわからない

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

29 では、どうするか…

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

31 ある家族 「我が家の行方不明者mは、遺体b1,...,bNのうちのどれか１体だと言えますか？それとも、言いかねますか？」
「mがbiである尤度」が「mがbi以外のどのbjである尤度」よりも、十分に大きければ、「m=bi」だと言える そのようなbiがなければ、「決めかねる」

32 ある家族 「mがbiである尤度」を計算しよう
「mがbiである尤度」が「mがbi以外のどのbjである尤度」よりも、十分に大きければ、「m=bi」だと言える 「mがbiである尤度」を計算しよう

33 N=3で考える ３ 体の遺体と３人の不明者を割り付ける場合の数 6通りの割り付け N=3 の場合 M1 M2 M3 仮説１ 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） 各仮説の尤度： ３つのペア全てでジェノタイプ が一致する確率

34 仮説１ 6通りの割り付け N=3 の場合 M1 M2 M3 仮説１ 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=B2）xP（M3=B3）

35 仮説2 6通りの割り付け N=3 の場合 M1 M2 M3 仮説１ 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）

36 仮説3 6通りの割り付け N=3 の場合 M1 M2 M3 仮説１ 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=B2）xP（M2=B1）xP（M3=B3）

37 仮説4 6通りの割り付け N=3 の場合 M1 M2 M3 仮説１ 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=B2）xP（M2=B3）xP（M3=B1）

38 仮説5 6通りの割り付け N=3 の場合 M1 M2 M3 仮説１ 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=B3）xP（M2=B1）xP（M3=B2）

39 仮説6 6通りの割り付け N=3 の場合 M1 M2 M3 仮説１ 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=B3）xP（M2=B2）xP（M3=B1）

40 M1=B1,それ以外の割り付けはなんでもあり 仮説1+2
6通りの割り付け N=3　の場合 M1 M2 M3 仮説１ 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）

41 M1=B1,それ以外の割り付けはなんでもあり 仮説1+2
6通りの割り付け N=3　の場合 M1 M2 M3 仮説１ 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）

42 N=3のとき 全割り付けの場合の数 3!=6 通り 「m1=b2、あとは何でもあり」の仮説 ２通り = (3-1)!

43 (N-1)!のこと M1=B1なら B1 B2 B3 M1 M2 M3 ここの割り付けは自由自在 P11 P12 P13 P21 P22

44 一般に N人 – N体 の割り付け そのうちの、１人 – １体の対応を固定する それ以外の割り付けはなんでもあり

45 (N-1)! 通りの「仮説」～「割り付け」の確率をすべて足し合わせる計算はできる？
はい、できます！ 近似法ですけど・・・

46 N=3で考える ３ 体の遺体と３人の不明者を割り付ける場合の数 6通りの割り付け N=3 の場合 M1 M2 M3 仮説１ 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） 各仮説の尤度： ３つのペア全てでジェノタイプ が一致する確率

47 仮説１ 6通りの割り付け N=3 の場合 M1 M2 M3 仮説１ 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=B2）xP（M3=B3）

48 仮説2 6通りの割り付け N=3 の場合 M1 M2 M3 仮説１ 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）

49 仮説3 6通りの割り付け N=3 の場合 M1 M2 M3 仮説１ 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=B2）xP（M2=B1）xP（M3=B3）

50 仮説4 6通りの割り付け N=3 の場合 M1 M2 M3 仮説１ 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=B2）xP（M2=B3）xP（M3=B1）

51 仮説5 6通りの割り付け N=3 の場合 M1 M2 M3 仮説１ 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=B3）xP（M2=B1）xP（M3=B2）

52 仮説6 6通りの割り付け N=3 の場合 M1 M2 M3 仮説１ 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=B3）xP（M2=B2）xP（M3=B1）

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

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

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

56 パーマネントの計算方法 今日は、割愛 気になる方は
Wikipedia から情報の入手は可能です 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= ){ 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)

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

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

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

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

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

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

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

64 NxN の確率行列 割り付けの計算のために正方行列がほしい。 全ての遺体と行方不明者が１対１対応すると仮定
仮定できなければ、足りない分を一般集団で補う 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）

65 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

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

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

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

69 シミュレーションとその結果 奈良原さんからゲット

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

71 家系の形 siblings の数 status の割り付け の割合は、1:10:1
genotype data 提供者とmissing が全員 siblings の関係（1世代家系） genotype data 提供者とmissing の少なくとも一組が親子の関係にある（2世代家系） genotype data 提供者とmissing の少なくとも一組が祖父母孫の関係にある（3世代家系） の割合は、1:10:1 (実際のデータの割合をもとに決めた) siblings の数 siblings の数が多くなると、計算量が多くなるので、1世代家系と2世代家系では3人まで、3世代家系では2人までと制限している。 status の割り付け 全てのstatus pattern を同確率で割り付け 結果、DNA提供者数が実際のデータより多めになってしまった感がある。（特に3世代家系で）

72 確率計算後 結局 全てのセルが０の家系を除外 全てのセルが０でない家系は除外しなかった
 残り３３家系 全てのセルが０でない家系は除外しなかった  ５４遺体は全て残っている 複数の行方不明者がいる家系のために2つのダミー遺体を追加 NxN にするために23のダミー家系を追加 ダミーの確率は正確法で計算したもの 結局 (33+23) 家系　vs (54+2) 遺体 でpermanent 計算へ

73 Simulation　結果

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

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

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

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

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

79 「事前確率」を取り込める 性別・所持品等の情報を容易に取り込める
Missings 血縁関係情報 血縁者のＤＮＡ型 付加情報 Bodies 遺体のＤＮＡ型 付加情報 確率行列 パーマネント計算 尤度割合行列 引き取り・引き渡し判定

80 申請のない行方不明者と発見されていない遺体のとりあつかい
?

81 正方行列を使う計算 N人のmissings, N人のbodies
NxN正方行列であるということは… すべての行方不明者が、家族に探されていて すべての行方不明者の遺体が発見されている ことが前提 行列を使った計算をしているということは… すべての行と列とが「対等」「独立」なことが前提 確率行列の計算でＤＮＡ多型のみを使っているということは… ＤＮＡ多型以外の情報は利用しないことが前提

82 NxN正方行列であるということは… 対応 すべての行方不明者が、家族に探されていて すべての行方不明者の遺体が発見されている ことが前提
捜索している家族のいない行方不明者も取り扱う 一般集団のジェノタイプ確率を使う 発見されていない遺体も、あるものとして取り扱う 行列は正方になる 行列が巨大化する 新たな工夫を検討中

83 確率行列の計算でＤＮＡ多型のみを使っているということは…
ＤＮＡ多型以外の情報は利用しないことが前提 対応 ありえない場合分けができれば行列を分割できる 処理は圧倒的に速くなる 例えば、性別を利用すると行列は半分になる 「ありえない」は事前確率が０である場合だが、「事前確率の重み」を取り込みたいときには、「確率行列」のセルにその重みを反映させることができる 発見場所や年齢情報、衣服情報などを使うことが対応する

84 「一般化」 個人鑑定を同じ枠組みで考えることが(おそらく)可能
B ? M

85 行列を使った計算をしているということは…
１家系に複数の行方不明者がいる場 非独立な確率・・・行列を使った計算が苦手 行列を使った計算をしているということは… すべての行と列とが「対等」「独立」なことが前提 １家族が複数のmissingsを探しているときには、P(mi=bj),P(mi’=b=j’)が独立ではないので、 確率行列がNxNにならない 確率行列から尤度割合行列を計算するときに、積による計算(パーマネントの前提)が使えない 対応 適当な対応策はない 独立を仮定して計算した上で、同一家系内の複数のmissingsに複数のbodiesを割り付けることが有望となったときには、その場合を個別に計算する →この計算はできる

86 ご協力 法医学教室 玉木先生 真鍋さん 問題の設定・シミュレーションデータの作成・計算結果の検算等、全面的にご教授・ご支援いただいております