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

Slides:



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

5 章 標本と統計量の分布 湯浅 直弘. 5-1 母集団と標本 ■ 母集合 今までは確率的なこと これからは,確率や割合がわかっていないとき に, 推定することが目標. 個体:実験や観測を行う 1 つの対象 母集団:個体全部の集合  ・有限な場合:有限母集合 → 1つの箱に入っているねじ.  ・無限な場合:無限母集合.
到着時刻と燃料消費量を同時に最適化する船速・航路計画
点対応の外れ値除去の最適化によるカメラの動的校正手法の精度向上
確率・統計Ⅰ 第12回 統計学の基礎1 ここです! 確率論とは 確率変数、確率分布 確率変数の独立性 / 確率変数の平均
多変量解析 -重回帰分析- 発表者:時田 陽一 発表日:11月20日.
近似アルゴリズム 第10章 終了時刻最小化スケジューリング
Pattern Recognition and Machine Learning 1.5 決定理論
Bassモデルにおける 最尤法を用いたパラメータ推定
アルゴリズムイントロダクション第5章( ) 確率論的解析
分子生物情報学 動的計画法に基づく配列比較法 (ペアワイズアライメント法)
Extremal Combinatrics Chapter 4
「データ学習アルゴリズム」 第3章 複雑な学習モデル 3.1 関数近似モデル ….. … 3層パーセプトロン
法数学勉強会(京大法医学講座) 2012/02/18 京都大学 統計遺伝学 山田
上坂吉則 尾関和彦 文一総合出版 宮崎大輔2003年6月28日(土)
Bias2 - Variance - Noise 分解
社会心理学のStudy -集団を媒介とする適応- (仮)
ベイズ的ロジスティックモデル に関する研究
IT入門B2 ー 連立一次方程式 ー.
確率・統計輪講資料 6-5 適合度と独立性の検定 6-6 最小2乗法と相関係数の推定・検定 M1 西澤.
最尤推定によるロジスティック回帰 対数尤度関数の最大化.
第3章 補足:パラメータが極小値に収束する例
多人数一括DNAプロファイリングのための 確率計算法に関する考察
法数学勉強会 2011/11/26 京都大学大学院医学研究科 統計遺伝学分野 山田 亮
高次元データの解析 -平均ベクトルに関する検定統計量の 漸近分布に対する共分散構造の影響-
遺伝的アルゴリズムへの 統計力学的アプローチ 大阪大学 大学院理学研究科 鈴木譲 CISJ2005 於早稲田大学理工学部
正規分布における ベーテ近似の解析解と数値解 東京工業大学総合理工学研究科 知能システム科学専攻 渡辺研究室    西山 悠, 渡辺澄夫.
決定木とランダムフォレスト 和田 俊和.
対立仮説下でのみ存在する 遺伝形式という母数を持つ 2x3分割表検定に関する考察 ~SNPによるケース・コントロール関連検定~
モデルの逆解析 明治大学 理工学部 応用化学科 データ化学工学研究室 金子 弘昌.
Probabilistic Method 2007/11/02 岩間研究室M1 市場孝之.
確率伝搬法と量子系の平均場理論 田中和之 東北大学大学院情報科学研究科
細胞の形と変形のための データ駆動型解析手法
独立成分分析 5 アルゴリズムの安定性と効率 2007/10/24   名雪 勲.
訓練データとテストデータが 異なる分布に従う場合の学習
ゲノム科学概論 ~ゲノム科学における統計学の役割~ (遺伝統計学)
法数学勉強会 2016年4月会 京都大学(医)統計遺伝学分野 山田 亮
法数学における ベイジアンネットワーク(2) ~成書で学ぶ~
25. Randomized Algorithms
法数学勉強会 2016/06/15 京都大学(医)統計遺伝学分野 山田 亮
法数学勉強会(京大法医学講座) 2012/02/18 京都大学 統計遺伝学 山田
東日本大震災におけるご遺体身元確認と行方不明家族捜索のためのDNA鑑定
2011/05/28 京都大学大学院 附属ゲノム医学センター統計遺伝学分野 山田 亮
第24回日本法科学技術学会 2018/11/08 京都大学 医学研究科 統計遺伝学分野 山田 亮
確率と統計2009 第12日目(A).
尤度の比較と仮説検定とを比較する ~P値のことなど~
第3章 線形回帰モデル 修士1年 山田 孝太郎.
法数学のための 機械学習の基礎 京大(医) 統計遺伝学分野 山田 亮 2017/04/15.
「アルゴリズムとプログラム」 結果を統計的に正しく判断 三学期 第7回 袖高の生徒ってどうよ調査(3)
情報経済システム論:第13回 担当教員 黒田敏史 2019/5/7 情報経済システム論.
決断のための分布合算 京大(医)統計遺伝学分野 山田 亮.
法医学会 2013年6月26日 京都大学(医)統計遺伝学 山田 亮
最尤推定・最尤法 明治大学 理工学部 応用化学科 データ化学工学研究室 金子 弘昌.
親子鑑定に見る尤度比を 角度を変えて眺めてみる
第16章 動的計画法 アルゴリズムイントロダクション.
法数学勉強会 2015/09/26 京都大学統計遺伝学分野 山田 亮
法数学勉強会 2015/09/26 京都大学統計遺伝学分野 山田 亮
パターン認識 ークラスタリングとEMアルゴリズムー 担当:和田 俊和 部屋 A513
藤田保健衛生大学医学部 公衆衛生学 柿崎 真沙子
音響伝達特性モデルを用いた シングルチャネル音源位置推定の検討 2-P-34 高島遼一,住田雄司,滝口哲也,有木康雄 (神戸大) 研究の背景
精密工学科プログラミング基礎 第7回資料 (11/27実施)
行列 一次変換,とくに直交変換.
Chapter5 Systems of Distinct Representatives
法数学における ベイジアンネットワーク(2) ~成書で学ぶ~
割り当て問題(assignment problem)
参考:大きい要素の処理.
精密工学科プログラミング基礎Ⅱ 第2回資料 今回の授業で習得してほしいこと: 配列の使い方 (今回は1次元,次回は2次元をやります.)
ガウシアングラフィカルモデルにおける一般化された確率伝搬法
混合試料の構成人数 Nuisance パラメタ
混合ガウスモデル Gaussian Mixture Model GMM
Presentation transcript:

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

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

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

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

行方不明者が 身元不明者のDNA型を持つ確率 この確率が計算できる?

行方不明者が 身元不明者のDNA型を持つ確率 はい、できます! 行方不明者 身元不明者 この確率が計算できる?

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

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

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

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

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

集団の鑑定 N人をN体に割り付ける .. .. .. .. .. .. 行方不明者 Missing 身元不明者 found Body m1 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(mN=b(sN)) ... ... 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!)

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

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!通りを計算しないで、最尤割り付けがわかる? 重みづけ最適化・重みづけマッチング ハンガリアン・アルゴリズムなど

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

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

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

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

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

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

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

では、どうするか…

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

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

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

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つのペア全てでジェノタイプ が一致する確率

仮説1 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=B2)xP(M3=B3)

仮説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)

仮説3 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=B2)xP(M2=B1)xP(M3=B3)

仮説4 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=B2)xP(M2=B3)xP(M3=B1)

仮説5 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=B3)xP(M2=B1)xP(M3=B2)

仮説6 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=B3)xP(M2=B2)xP(M3=B1)

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)

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

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

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

(N-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つのペア全てでジェノタイプ が一致する確率

仮説1 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=B2)xP(M3=B3)

仮説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)

仮説3 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=B2)xP(M2=B1)xP(M3=B3)

仮説4 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=B2)xP(M2=B3)xP(M3=B1)

仮説5 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=B3)xP(M2=B1)xP(M3=B2)

仮説6 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=B3)xP(M2=B2)xP(M3=B1)

行列式(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: 行列サイズ)

パーマネント計算 近似法の精度 「真のパーマネント」を大きな行列で計算するのは非現実的なので、「真の精度」を評価するのは難しいのですが パーマネント計算 近似法の精度 「真のパーマネント」を大きな行列で計算するのは非現実的なので、「真の精度」を評価するのは難しいのですが 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人をランダムに選択  家系も遺体も足りない状況

家系の形 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世代家系で)

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

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人 正解ペア

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

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

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

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

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

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

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

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

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

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

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