2005-05-28

Dirichlet分布、ぼーず確率

ぼーず確率については、ryamadaの遺伝学・遺伝統計学メモ - Coalescence発生待ち時間の算出で書いたが・・・。その発展編。

Dirichlet分布一般については、関連記事(こちらとこちら)も参照。

ながらくイメージのつかめなかった、Dirichlet分布がようやくわかったような気がする。Dirichlet分布は、シミュレーションで期待値に基づいてサンプリングをするときに避けて通れない分布でしたが、ながらく、なんとなく使ったり、回避したりしてきました。

ここから本論＞

相互に独立して発生する事象がある。それぞれある期待値をもって発生する。それらを全部併せて1000回起こすことにする。そのようなときに用いるのがDirichlet分布である(ディリクレ分布)。

今、起こりやすさが２：３：５の３つの事象が起きて、併せて1000回起きるとする。無限回数繰り返すとすると、その生起回数の比は２：３：５となるが、有限回数のときには、それがぶれる。ちょうど、魚釣り競争で、平均漁獲量が２：３：５の３人が漁獲量を競っているときに、どの人も優勝する確率があれば、どの人もびりになる確率があるのと同じことである。したがって、このように有限時間内での勝負のときには、ぼーず確率で見たように、平均釣り上げ量が組み入れられた指数関数が登場する。

ガンマ関数は、添付の図の右肩のような式をしているが、Rのガンマ関数を用いる事とする。

Rで

>rgamma(N,x)

と計算させると、漁獲量が平均ｘとなるような時間、魚釣りをしたときに、何匹釣れるかをN回分、試行してくれる。答えは有理数で返ってくる。

rgammaはこんな風に書かれている。



function (n, shape, rate = 1, scale = 1/rate) 

{

    if (any(shape <= 0)) 

        stop("shape must be strictly positive")

    .Internal(rgamma(n, shape, scale))

}

すると、今、3人が、ある制限時間内に2匹、3匹、5匹つることが期待されているときに、それぞれが、釣る魚の数はrgamma(1,2),rgamma(1,3),rgamma(1,5)でシミュレーションされる。一発勝負ならば、この３つの数値を較べれて終わり。

では、繰り返しこの勝負をするとしたときに、それぞれ勝負で、何匹対何匹対何匹になるかをN回分の勝負でデータを見たいとすると、

rgamma(N,2),rgamma(N,3),rgamma(N,5)

で得られる3つの数値のNセットがほしいデータである。

もし、釣った魚の引数ではなくて、3人の釣果の比率で比較したいとしたときに、3人の総釣果への寄与程度をN回分のデータとして得る必要があるときに、Dirichlet分布を用いる。ちなみにrgamma(1000,2)では、1000個の乱数が作成される。その平均は２に近く、rgamma(10000,2)とすれば、ますます、この10000個の乱数の平均は２に近づく

Dirichlet分布も、式を見た方が早い。1000回勝負での釣果比率を求めるとすると・・・

Rの"MCMCpack" Packageにある、Dirichletを使うとして、

まず、"MCMCpack"パッケージをダウンロードしてインストール。

Rを起動して

>library(MCMCpack)

ついて

>rdirichlet(1000,c(2,3,5))

で得られる。

この結果(1000回の試行)のそれぞれにおいて、２の人、３の人、５の人がそれぞれ、何割ずつ漁獲量を占めたかをヒストグラムにしたのが、図、である。２の人は2割前後、３の人は3割前後、５の人が5割前後になっているのがわかる。

もし、今、個々の試行の時間を100倍に伸ばし、２の人は平均200匹、３の人は300匹、５の人は500匹、釣り上げるような時間で分布をとりなおすためには、

>rdirichlet(1000,c(200,300,500))

とする。このとき、長時間かけているので、５の人が優勝する確率はずっと高くなり、２の人がびりになる確率もずっと高くなるのは、当然といえる。それが、ヒストグラムだとどのように現れるか、というと、それぞれの人のヒストグラムが期待値付近に高いピークを持ち、裾野が狭くなる。２の人は毎回毎回、ほぼ２割を釣るだけで、今の図にあるように、４割を超えることなどほとんどなくなる。

ついでに、備忘のために上記Rの処理結果をファイルに書き出してエクセルでヒストグラムまでにする過程も書いておくと

>diriout<-rdirichlet(1000,c(2,3,5))

>write(diriout,file="dirioutfile",ncolumns=1000)

こうすると、スペース区切りで1000個の有理数を１行とする3行のファイルができる。

行列変換プログラム(末尾にソース)で変換すると、３列1000行のファイルになるので、それをエクセルで開いて、列ごとにFREQUENCY関数を使って、0.05刻みのヒストグラム表示をした。

ちなみに、R MCMCpack の　rdirichlet関数の中身は



function (n, alpha) 

{

    l <- length(alpha)

    x <- matrix(rgamma(l * n, alpha), ncol = l, byrow = TRUE)

    sm <- x %*% rep(1, l)

    return(x/as.vector(sm))

}

行列置換プログラム

#! /usr/bin/perl -w
=head1 header
 ###########################################################

 #  
 #
=cut
use strict;


&main(@ARGV);

exit;
=head1 main
 ###########################################################

 #

 #      #main関数

 #      #    \arg1        $input_file : format_fixed input file

 #      #    \arg2        $output_file : format_fixed output file

 #      tab区切りファイルの行列置換

 #

 

=cut


sub main(@){

  my ($input_file,$output_file)=@_;

 open (INPUT, "$input_file");

 open (OUT, ">>$output_file");

 my @twodim;

 my $linecount=0;

 my $columncount=0;
 while(<INPUT>){

  chomp;

  my @inline = split(/\s/);

  if($columncount<@inline){

   $columncount=@inline;

  }

  for(my $i=0;$i<@inline;$i++){

   $twodim[$linecount][$i]=$inline[$i];

  }

  $linecount++;

 }
 for(my $i=0;$i<$columncount;$i++){

  for(my $j=0;$j<$linecount;$j++){

   print OUT "$twodim[$j][$i]\t";

  }

  print OUT "\n";

 }

 close OUT;

}