2005-10-19

MCMCによるDirichlet random generation再考

作業は統計環境『Ｒ』を用いるものとする。ガンマ分布や、Dirichlet分布からの乱数発生関数 rdirichlet (MCMCpack)を用いる(その使い方についてはこちら)

サイコロ
- 今、サイコロの「できのよさ」を調べながら、そのサイコロで博打をうつことを想定する
- サイコロであるからには、６つの目がほぼ同確率で出るようになっている。その確率はR

で言えば、次のように６つの目にそれぞれ同確率をもつベクターとして作ることができる

dice<-c(1/6,1/6,1/6,1/6,1/6,1/6)

理想的なサイコロ
- 今、サイコロを振って、何の目が出るかを予想している。サイコロが理想的であれば、このdiceベクターを信じて、それぞれの目が出る確率は $￥Large ￥frac{1}{6}$ だと信じて予想すればよい
- この「理想的なサイコロ」と「信じる」とはどういうことだろうか
  - いつサイコロを振っても、必ず、出る目は１から６まで $￥Large ￥frac{1}{6}$ である
  - 無限回、サイコロを振ると１から６の目の出る回数は $￥Large ￥frac{￥infty}{6}$ である
  - では、このように「理想的なサイコロ」を振って次に１から６のそれぞれの目が出る確率を予想しよう。もちろん、「理想的」と「知っている」ので、それぞれである。
    - これをDirichlet分布からの乱数発生としてやらせるとなるとどうなるだろうか
      - rdirichlet(1,c( $￥Large ￥frac{￥infty}{6}$ , $￥Large ￥frac{￥infty}{6}$ , $￥Large ￥frac{￥infty}{6}$ , $￥Large ￥frac{￥infty}{6}$ , $￥Large ￥frac{￥infty}{6}$ , $￥Large ￥frac{￥infty}{6}$ ))
      - この式のは、 $￥Large ￥frac{1}{6}$ に正確に一致する数値を６つ返してくる(はずである。無限大を受け付けないので、実際には、非常に大きい数値で代用して実行する)
      - この式の意味は次の通りである
        
        今、６個の事象(サイコロの目の１が出たというのが、１個の事象)が $￥Large ￥frac{￥infty}{6}$ 回ずつ観測されたとしたときに、その観測結果に基づいて、 $￥infty +1$ 回目に１、２、３、４、５、６のそれぞれの目が出る確率を順番に返せ、６つの事象のそれぞれが観測されるだろう確率だから、６つの数値を足し合わせると１になる
      - 実際にＲで乱数発生を行わせるには、k 回の試行をして、その分布をとるとしたら
        
        hist(rdirichlet(k,dice*10000)[,1],xlim(range(0,1))
      - などとしてやることで、１の目が出る確率がどのくらい $￥frac{1}{6}$ に集中するかが図示できる。dice*N のN の値を無限大らしく大きくしてやるとその集中度が増していく様子も確認できる
        
        hist(rdirichlet(k,dice*10000)[,t],xlim=range(0,1)) としてやれば、tの目が出る確率の分布が得られる
現実のサイコロ
- では、現実のサイコロで、１から６の目の出る確率が均一でないものとする。どのくらい均一でないかを知りたいとする
- よくやるとおり、何度もサイコロを振って出た目の記録を残して行く
  - その結果、18,20,16,12,18,16 (足して100回)となったとする。
- この結果を見て、このサイコロは、２の目が一番でやすい傾向があり、４の目が一番でにくい傾向があるように思われる
- ただし、それは、１の目と５の目の出る確率が、きっかり一緒だと予想することが正しくないだろうということや、それぞれの目の出る確率がそれぞれ正確に0.18,0.20,0.16,0.12,0.18,0.16であると予想するのが危険だろうことからもわかる通り、「おおまかな傾向」を示しているだけで、正確なことはわからない
- しかしながら、この100回の実験をしてみる前よりもなにがしかの情報が得られていることは確かである
- このような事前情報が得られているときに、次(101回目)に１から６の目のどれがどのくらいの確率で出ると
- 上述の「理想のサイコロ」の伝で行くと次のコマンドで、１の目が出る確率をｋ回出してみた場合の分布が得られる。
  - hist(rdirichlet(k,c(18,20,16,12,18,16))[,1],xlim=range(0,1))
- もちろん、１の目が出る確率の分布の平均値は0.18に近くなる。そのことは次のコマンドで確認できる
  - mean(rdirichlet(100000,c(18,20,16,12,18,16))[,1])
理想的なサイコロ、再び
- 今、実は、「理想的」であるようなサイコロがあるとする。残念ながら、それが理想的なサイコロかどうかは人智では知る由もない。
- 実際にこのサイコロのよさを調べるとなると、上述の「現実のサイコロ」と同様な実験をするしかない。
- 理想的なサイコロでは、目の出方はばらばらしているが、「無限回やったら $￥Large ￥frac{1}{6}$ に収束していくことは確かだが、その無限回の途中には、目の出方にはばらつきがある
- したがって、サイコロを100回投げたところ、上述の「現実のサイコロ」のような目の結果が得られる可能性もある
- 仮に、集計のたびに、きっちり、ずつの集計結果が出たとする。その集計の度に、そのとき得られている情報でそのサイコロの目の出る確率に見当をつけ、次に出る目の確率を予想するとするとその作業は次のようにrdirichlet関数を使って示すことができる
  - １の目が出る確率を、全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))
  - となる。いずれの場合も、ｋ回の試行の平均は、 $￥Large ￥frac{1}{6}$ 付近になるが、その分布の山のまとまりのよさは、60回→600回→600000回と回を増やすほど、『理想＝いつでも $￥Large ￥frac{1}{6}$ 』に近づいていくことがわかる
より現実的には
- かなり精度のよいサイコロであると信じるに足るサイコロを持っているとする。みるからに立方体だし、いかさま用に作られたとは思えない、小学2年生1月号の付録のすごろくについて来るようなサイコロのような場合である
- 「かなり精度のよいサイコロであると信じるに足る」ということはどういうことかというと
  - 『実験してみるまでもなく』dice<-c(100,100,100,100,100,100)と言っていいんじゃないかと思わせられる、ということ
  - 意訳すると、ためしに100回振ってみたら、１が出る確率が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))はエラーメッセージが出る。１の目が出た回数が０(１の目の出る事前情報が確率０)だからである。Ｒの関数的には、rdirichlet関数が呼び出しているrgamma分布の関数にこの０が渡されるのだが、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本の親ハプロタイプとその組換え体としてのハプロタイプとの関係についての尤度を計算するときに、親ハプロタイプの頻度分布を直接利用せず、頻度に単純に１を加えた値を事前dirichlet分布としている
    - これらSNPHAP,TreeLDでは１を加えているが、上述の『1月号の付録のサイコロ』に(0,4,2,1,2,1)が出たときには、それぞれ１を加えて(1,5,3,2,3,2)をrdirichlet関数に渡すのが適当か、それぞれ100を加えて(100,104,102,101,102,101)からdirichlet乱数を発生させるのがよいのかは、考えどころのように思う