MCMCによるDirichlet random generation再考
作業は統計環境『R』を用いるものとする。ガンマ分布や、Dirichlet分布からの乱数発生関数 rdirichlet (MCMCpack)を用いる(その使い方についてはこちら)
- サイコロ
- 今、サイコロの「できのよさ」を調べながら、そのサイコロで博打をうつことを想定する
- サイコロであるからには、6つの目がほぼ同確率で出るようになっている。その確率はR
で言えば、次のように6つの目にそれぞれ同確率をもつベクターとして作ることができる
- dice<-c(1/6,1/6,1/6,1/6,1/6,1/6)
- 理想的なサイコロ
- 今、サイコロを振って、何の目が出るかを予想している。サイコロが理想的であれば、このdiceベクターを信じて、それぞれの目が出る確率は だと信じて予想すればよい
- この「理想的なサイコロ」と「信じる」とはどういうことだろうか
- いつサイコロを振っても、必ず、出る目は1から6まで である
- 無限回、サイコロを振ると1から6の目の出る回数は である
- では、このように「理想的なサイコロ」を振って次に1から6のそれぞれの目が出る確率を予想しよう。もちろん、「理想的」と「知っている」ので、それぞれ である。
- これをDirichlet分布からの乱数発生としてやらせるとなるとどうなるだろうか
- rdirichlet(1,c(,,,,,))
- この式のは、に正確に一致する数値を6つ返してくる(はずである。無限大を受け付けないので、実際には、非常に大きい数値で代用して実行する)
- この式の意味は次の通りである
- 今、6個の事象(サイコロの目の1が出たというのが、1個の事象)が 回ずつ観測されたとしたときに、その観測結果に基づいて、 回目に1、2、3、4、5、6のそれぞれの目が出る確率を順番に返せ、6つの事象のそれぞれが観測されるだろう確率だから、6つの数値を足し合わせると1になる
- 実際にRで乱数発生を行わせるには、k 回の試行をして、その分布をとるとしたら
- hist(rdirichlet(k,dice*10000)[,1],xlim(range(0,1))
- などとしてやることで、1の目が出る確率がどのくらい に集中するかが図示できる。dice*N のN の値を無限大らしく大きくしてやるとその集中度が増していく様子も確認できる
- hist(rdirichlet(k,dice*10000)[,t],xlim=range(0,1)) としてやれば、tの目が出る確率の分布が得られる
- これをDirichlet分布からの乱数発生としてやらせるとなるとどうなるだろうか
- 現実のサイコロ
- では、現実のサイコロで、1から6の目の出る確率が均一でないものとする。どのくらい均一でないかを知りたいとする
- よくやるとおり、何度もサイコロを振って出た目の記録を残して行く
- その結果、18,20,16,12,18,16 (足して100回)となったとする。
- この結果を見て、このサイコロは、2の目が一番でやすい傾向があり、4の目が一番でにくい傾向があるように思われる
- ただし、それは、1の目と5の目の出る確率が、きっかり一緒だと予想することが正しくないだろうということや、それぞれの目の出る確率がそれぞれ正確に0.18,0.20,0.16,0.12,0.18,0.16であると予想するのが危険だろうことからもわかる通り、「おおまかな傾向」を示しているだけで、正確なことはわからない
- しかしながら、この100回の実験をしてみる前よりもなにがしかの情報が得られていることは確かである
- このような事前情報が得られているときに、次(101回目)に1から6の目のどれがどのくらいの確率で出ると
- 上述の「理想のサイコロ」の伝で行くと次のコマンドで、1の目が出る確率をk回出してみた場合の分布が得られる。
- hist(rdirichlet(k,c(18,20,16,12,18,16))[,1],xlim=range(0,1))
- もちろん、1の目が出る確率の分布の平均値は0.18に近くなる。そのことは次のコマンドで確認できる
- mean(rdirichlet(100000,c(18,20,16,12,18,16))[,1])
- 理想的なサイコロ、再び
- 今、実は、「理想的」であるようなサイコロがあるとする。残念ながら、それが理想的なサイコロかどうかは人智では知る由もない。
- 実際にこのサイコロのよさを調べるとなると、上述の「現実のサイコロ」と同様な実験をするしかない。
- 理想的なサイコロでは、目の出方はばらばらしているが、「無限回やったら に収束していくことは確かだが、その無限回の途中には、目の出方にはばらつきがある
- したがって、サイコロを100回投げたところ、上述の「現実のサイコロ」のような目の結果が得られる可能性もある
- 仮に、集計のたびに、きっちり、 ずつの集計結果が出たとする。その集計の度に、そのとき得られている情報でそのサイコロの目の出る確率に見当をつけ、次に出る目の確率を予想するとするとその作業は次のようにrdirichlet関数を使って示すことができる
- 1の目が出る確率を、全60回の実験後、全600回の実験後、...、全600000回の実験後に予想して分布をとるとすると
- hist(rdirichlet(k,dice*100)[,1],xlim=range(0,1))
- hist(rdirichlet(k,dice*1000)[,1],xlim=range(0,1))
- hist(rdirichlet(k,dice*1000000)[,1],xlim=range(0,1))
- となる。いずれの場合も、k回の試行の平均は、 付近になるが、その分布の山のまとまりのよさは、60回→600回→600000回と回を増やすほど、『理想=いつでも 』に近づいていくことがわかる
- 1の目が出る確率を、全60回の実験後、全600回の実験後、...、全600000回の実験後に予想して分布をとるとすると
- より現実的には
- かなり精度のよいサイコロであると信じるに足るサイコロを持っているとする。みるからに立方体だし、いかさま用に作られたとは思えない、小学2年生1月号の付録のすごろくについて来るようなサイコロのような場合である
- 「かなり精度のよいサイコロであると信じるに足る」ということはどういうことかというと
- 『実験してみるまでもなく』dice<-c(100,100,100,100,100,100)と言っていいんじゃないかと思わせられる、ということ
- 意訳すると、ためしに100回振ってみたら、1が出る確率がhist(rdirichlet(100,c(100,100,100,100,100,100)[,1]))で得られるくらいには正確なサイコロ、だろうと信じている、ということ
- このようなサイコロだけれども、やはり、本当はどれくらい「よくできたサイコロ」かを知りたいから、実験してみるとする
- 100回投げたら、(18,20,16,12,18,16)と出た。このときは、rdirichlet(k,c(18,20,16,12,18,16))で101回目からの目の出る確率を予想してよさそうに思われる
- では、10回だけ投げて(0,4,2,1,2,1)と出たとしよう。このときはどうだろうか
- そもそもrdirichlet(k,c(0,4,2,1,2,1))はエラーメッセージが出る。1の目が出た回数が0(1の目の出る事前情報が確率0)だからである。Rの関数的には、rdirichlet関数が呼び出しているrgamma分布の関数にこの0が渡されるのだが、rgamma関数の引数shapeが0で定義されていないからである。数学的にガンマ分布の第1引数は0より大であるという定義であるからとも言い換えられる
- この問題は、統計遺伝学でMCMCを用いるときに事前分布から事後分布を作るときにも発生する問題で、その例としては phase-unknownなSNPのdiploid genotypeデータからのハプロタイプ推定にもあり、それについてのコメントは、アプリケーションsnphapの解説文にもある。当該解説文はこちら。
- 該当箇所の抜粋は以下の通り
- If the prior is a Dirichlet distribution with constant df on all possible haplotypes*, the full conditional posterior distribution is also Dirichlet. To obtain the set of df parameters for this posterior Dirichlet, we simply add the constant prior df to the number of chromosomes currently assigned to each haplotype.
- 該当箇所の抜粋は以下の通り
- 本メモの記事で扱っているZollnerらのMCMCベースのARG構築アルゴリズムTreeLDでも同様にMCMCにて事後分布を発生させるにあたり、2本の親ハプロタイプとその組換え体としてのハプロタイプとの関係についての尤度を計算するときに、親ハプロタイプの頻度分布を直接利用せず、頻度に単純に1を加えた値を事前dirichlet分布としている
- これらSNPHAP,TreeLDでは1を加えているが、上述の『1月号の付録のサイコロ』に(0,4,2,1,2,1)が出たときには、それぞれ1を加えて(1,5,3,2,3,2)をrdirichlet関数に渡すのが適当か、それぞれ100を加えて(100,104,102,101,102,101)からdirichlet乱数を発生させるのがよいのかは、考えどころのように思う