確率質量分布で表すRandom Exchangeable Partitions - ryamadaの遺伝学・遺伝統計学メモ

Kingmanのpaintbox〜単位線分タイリングがrandom exchangeable partitionsの表現であることがわかったが、実際、どんなrandom exchangeale partitionsが現れるのかは確率事象なので
確率事象として生成されるときにどんな確率事象なのか、そのときにどのような分布になるのか、ということが気にかかる
前の記事では有限個線分の個数を指定し、一様ディリクレ乱数によって線分の長さを定めて生成した。この方法だと、有限個線分の指定個数がかなり大きな制約条件とはなっているが、これも、一つの生成確率事象ルールではある
以下では、いくつかの方法を記述する
どんな条件があるかと言えば、 $\sum_{i=0}^\infty p_i=1;p_i \ge 0$ 。確率が0より大とすれば確実に無限個に分割。確率が0のものが無限個あるとすれば、有限個への分割もこの定義でカバーされる
自然数・非負整数を大とする確率質量分布は使える
- 有限個の項に対応する多項分布の確率質量に比例した単位線分分割は、Kingmanのタイリングの無限小分割を持たないバージョンとなる

k <- 10
sort(rdirichlet(1,rep(1,k)))

- ポアソン分布の確率質量は0,1,...,無限大のそれぞれに正の確率質量が与えられその和は1

plot(dpois(0:100,3),type="h")
plot(sort(dpois(0:100,3)),type="h")

- ポアソン分布よりも分散が大きい確率質量分布に負の二項分布がある

plot(dnbinom(0:100,3,0.3),type="h")
plot(sort(dnbinom(0:100,3,0.3)),type="h")

- ディリクレ分布、ポアソン分布、負の二項分布のいずれも、ガンマ分布を使って次のように表現できる
  - k項のディリクレ乱数は、尺度パラメターが1のガンマ分布に従うk個の乱数 $\Gamma_1,...,\Gamma_k$ を発生し $\frac{\Gamma_i}{\sum_{i=1}^k \Gamma_i}$ とすることで生成できる。このとき $\sum_{i=1}^k \Gamma_i$ は有限の値を取ることが知られている。このときのガンマ分布の形状パラメタベクトルがディリクレ分布の形状パラメタになっている(らしい)
  - ポアソン分布の場合は1個のガンマ分布乱数をその期待値として指定することでガンマ分布と紐づく
  - 負の二項分布の場合は、ガンマ乱数を発生させ、その乱数に対応するポアソン分布から乱数を一つ取り出すことを繰り返すとその結果が負の二項分布となると知られています。負の二項分布のsizeとprobとのパラメタとガンマ分布のshape,rateパラメタとの関係はshape = size, rate = prob/(1-prob)です)こちら)

N <- 10^4
a <- c(0.1,0.2,0.3)
Rdirichlet <- rdirichlet(N,a)
Rgamma <- cbind(rgamma(N,a[1],1),rgamma(N,a[2],1),rgamma(N,a[3],1))
Rgamma <- Rgamma/apply(Rgamma,1,sum)
plot(sort(Rdirichlet[,1]),sort(Rgamma[,1]))

N <- 10^5
sz <- 4
pr <- 0.3
sp <- sz
rt <- pr/(1-pr)
Rnb <- rnbinom(N,sz,pr)
RgammaPois <- rep(0,N)
for(i in 1:N){
	tmp <- rgamma(1,sp,rt)
	RgammaPois[i] <- rpois(1,tmp)
}

plot(sort(Rnb),sort(RgammaPois))