2014-04-29

ディリクレ分布をいじる

kカテゴリあって、全部でN回サンプリングしたところ、それぞれが $\mathbf{a}=(a_1,...,a_k)$ 回ずつ観測されたとする。ただし $\sum_{i=1}^k a_i=N$
このとき、サンプリング元のkカテゴリの割合を $\mathbf{x}=(x_1,...,x_k);\sum_{i=1}^k x_i=1$ は確率変数であって、そのような確率は $P(\mathbf{x}|\mathbf{a})=\frac{1}{B(\mathbf{a}+1=(a_1+1,a_2+1,...,a_k+1))}\prod_{i=1}^k x_i^{a_i}$ に従うと想定するとよいことが知られていて、この分布をディリクレ分布と言う。
観測数を $(a_1-1,a_2-1,...,a_k-1)$ と置くことの方が一般的で、この場合、0回の観測を $a_i-1=0$ から $a_i=1$ に対応づけるが、こちらを採用すると $P(\mathbf{x}|\mathbf{a})=\frac{1}{B(\mathbf{a})} \prod_{i=1}^k x_i^{a_i-1}$ と表される
例を扱っておこう
- k=3タイプある場合に、タイプ1が３個、タイプ2が2個、タイプ3が１個、観察されたとすれば、3タイプの割合が0.5,0.3,0.2であるのは

library(MCMCpack)
obs <- c(3,2,1)
pr <- c(0.5,0.3,0.2)
ddirichlet(pr,obs+1)

次のようにしてもよい

my.beta <- function(a){
	exp(sum(lgamma(a))-lgamma(sum(a)))
}
1/my.beta(obs+1) * prod(pr^obs)

今、ある観察をすると、それに応じて、サンプリング元の割合の確率密度分布がディリクレ分布として得られるわけであるが、じゃあ、そんな確率密度分布から、実際に観察データが得られる確率の期待値を得る、となったらどうなるだろうか
$\int_{\math{x}\in \Omega_k} \frac{1}{B(\mathbf{a})} \prod_{i=1}^k x_i^{a_i-1} d \mathbf{x}=1$ は、確率密度関数の定義から言える。ただし $\Omega_k$ はkカテゴリの比率が取りうるk-1正単体空間である。
それぞれの $\mathbf{x}$ において、 $n=\sum_{i=1}^k (a_i-1)$ 回のうち、 $\mathbf{a}=(a_1-1,...,a_k-1)$ と観測する確率は $\frac{n!}{\prod_{i=1}(a_i-1)!}\prod_{i=1}^k x_i^{a_i-1}=C(\mathbf{a}) \prod_{i=1}^k x_i^{a_i-1}$ であるから
$\int_{\mathbf{x}\in \Omega_k} \frac{1}{B(\mathbf{a})} \prod_{i=1}^k x_i^{a_i-1} C(\mathbf{a}) \prod_{i=1}^k x_i^{a_i-1} d \mathbf{x}=\int_{\math{x}\in \Omega_k} \frac{C(\mathbf{a})}{B(\mathbf{a})}\prod_{i=1}^k x_i^{(2a_i -1)-1} d\mathbf{x}$
これを少しいじると
[tex:\frac{C(\mathbf{a})}{B(\mathbf{a})}B*1\int_{\math{x}\in \Omega_k} \frac{1}{B*2}\prod_{i=1}^k x_i^{(2a_i -1)-1} d\mathbf{x}]となり、この積分部分は、ディリクレ分布の $\Omega_k$ 全体の積分なので1であるから、積分は消えて
$\frac{C(\mathbf{a})}{B(\mathbf{a})}B(2\mathbf{a}-1)$ となる
今、ある観察データに関して $C(\mathbf{a})$ は共通なので、この値の比を問題にする場合には無視できる
ここで、カテゴリ数がの場合との場合とを考えることとする。今、確かにkカテゴリは観察されているが、それ以上のカテゴリが観察していないとすると、
- $\mathbf{x}=(x_1,...,x_k)$ というk次元の確率ベクトルと $\mathbf{y}=(y_1,...,y_k,y_{k+1},...,y_K)$ という確率ベクトルとの2つを考えることとなる。
ここで先ほどの $C(\mathbf{a})$ を無視した $\frac{B(2\mathbf{a}-1)}{B(\mathbf{a})}$ を、kの場合とKの場合とで比較しよう
ここで、観察サンプル総数をn、その観察内訳を $\mathbf{a}=(a_1-1,...,a_k-1);a_i -1 > 0$ とすれば、kカテゴリのすべてが少なくとも１回以上、観察されたことになっている。ここで、K>kカテゴリを想定しているから、 $\mathbf{b}=(a_1-1,...,a_k-1,b_{k+1}-1=0,...,b_{K}=0)$ とする。これに対応する、 $\frac{B(2\mathbf{b}-1)}{B(\mathbf{b})}$ は次のようになる
$\frac{\prod_{i=1}^K \Gamma(2b_i-1)}{\Gamma(\sum_{i=1}^K (2b_i-1))}\frac{\Gamma(\sum_{i=1}^K b_i)}{\prod_{i=1}^K \Gamma(b_i)}$
ここで $\Gamma(1)=1$ 、 $\Gamma(2\times 1-1)=\Gamma(1)=1$ であることに注意する。これは観測されているカテゴリはkカテゴリであるが、1回も観測されていないK-kカテゴリについての項である。このK-kカテゴリが、今考えている比に与える影響は、 $\Gamma(b_i)$ も $\Gamma(2 b_i-1)$ もその値は1であって、影響を与えない、ということを意味する。
また、 $\Gamma(a_i)=\Gamma(b_i=a_i); i \le k$ であることにも注意すれば、 $\prod_{i=1}^k \Gamma(a_i)=\prod_{i=1}^K \Gamma(b_i)$ であり、 $\prod_{i=1}^k \Gamma(2a_i-1)=\prod_{i=1}^K \Gamma(2b_i-1)$ である。
したがって
$\frac{\frac{B(2\mathbf{a}-1)}{B(\mathbf{a})}}{\frac{B(2\mathbf{b}-1)}{B(\mathbf{b})}}$ は、観測総数と考慮するカテゴリ数との和とを用いて表されるガンマ関数4つだけが残る。この意味するところは、「観測総数」と「カテゴリ数」とが決めるのであって、観測カテゴリの個々のカテゴリ別観測数は、この比に影響しないことを意味している。
$\frac{\Gamma(n+K)}{\Gamma(2n+K)}\frac{\Gamma(2n+k)}{\Gamma(n+k)}$ となる

my.beta <- function(a){
	exp(sum(lgamma(a))-lgamma(sum(a)))
}

$\frac{(\sum a_i)!}{\prod a_i!}$

my.comb <- function(a){
	exp(lgamma(sum(a)+1)-sum(lgamma(a+1)))
}

$C(\mathbf{a})\frac{B(2\mathbf{a}-1)}{B(\mathbf{a}}$

my.CBB <- function(a){
	my.like.dir(a)*exp(lgamma(sum(a)-length(a)+1)-sum(lgamma(a)))
}

組合せ項は比を取る時は不要なので取り除いたもの

my.like.dir <- function(a){
	my.beta(2*a-1)/my.beta(a)
}
# 観測度数bを使って比を上記の比を計算
my.like.dir.obs <- function(b){
	if(length(b)==1){
		return(1)
	}
	a <- b+1
	my.like.dir(a)
}
# サンプル数n、観察カテゴリ数k、想定カテゴリ数Kの場合のk, K 頂点正単体空間の発生確率総積分
my.ratio.dir.integral <- function(n,k,K){
	exp(lgamma(n+K)-lgamma(2*n+K)-lgamma(n+k)+lgamma(2*n+k))
}

*1:2a_1-1,2a_2-1,...,2a_k-1

*2:2a_1-1,2a_2-1,...,2a_k-1