多次元確率分布　ディレクレ分布

各成分が正の$${d}$$次元のベクトル$${{\bf \alpha}=(\alpha_1,\cdots, \alpha_d)^T}$$を用いたガンマ分布に従う確率変数$${{\bf y}=(y_1,\cdots,y_d)^T}$$を、$${y_{tot}=\Sigma^d_{j=1}y_j}$$で正規化した確率変数$${{\bf x}=(x_1,\cdots,x_d)^T=\displaystyle{(\frac{y_1}{y_{tot}},\cdots,\frac{y_d}{y_{tot}})^T}}$$が従う分布が、ディレクレ分布である。
定義から明らかに、$${x_i \geq 0: i=1,\cdots, d }$$であり、$${\Sigma^d_{j=1}x_j=1}$$である。
このディレクレ分布に従う確率変数は、$${d}$$面の偏ったサイコロと考えられ、その確率密度関数は、$${d}$$次元のベータ関数、
$${B_d({\bf \alpha})=\displaystyle{\int_{{\bf x}} d{\bf x}\Pi^d_{j=1}x_j^{\alpha_j-1}}}$$
を用いて、
$${f(x)=\displaystyle{\frac{\Pi^d_{j=1}x_j^{\alpha_j-1}}{B_d({\bf \alpha})}}}$$
と与えられている。
　この$${d}$$次元のベータ関数は、元の一次元のベータ関数$${B(\alpha, \beta)=\displaystyle{\int_0^1x^{\alpha-1}(1-x)^{\beta-1}dx}}$$において、積分範囲を$${0}$$から$${p (0< p< 1)}$$として、
$${\displaystyle{\int_0^px^{\alpha-1}(p-x)^{\beta-1}dx=\int_0^{\frac{\pi}{2}}p^{\alpha+\beta-2}\sin^{2(\alpha-1)}x \cos^{2(\beta-1)}x2p\sin x\cos xdx}}$$
ここで、$${x=p\sin^2x}$$と置いている。
$${=\displaystyle{2p^{\alpha+\beta-1}\int_0^{\frac{\pi}{2}}\sin^{2\alpha-1}x \cos^{2\beta-1}x dx=p^{\alpha+\beta-1}B(\alpha, \beta)}}$$
となることを用いる。
各$${x_j}$$の積分範囲を明確に表示すると、
$${B_d({\bf \alpha})=\displaystyle{\int_0^1 x_1^{\alpha_1-1}dx_1 \int_0^{1-x_1} x_2^{\alpha_2-1}dx_2 \cdots \int_0^{1-\sum^{d-2}_{j=1} x_{j}}x_{d-1}^{\alpha_{d-1}-1} (1-\sum^{d-1}_{j=1}x_j)^{\alpha_d-1} dx_{d-1}}}$$
となる。
ここで、最後の$${d-1}$$についての積分で、$${p=1-\sum^{d-2}_{j=1}x_j}$$とおけば、
$${\displaystyle{ \int_0^{1-\sum^{d-2}_{j=1} x_{j}}x_{d-1}^{\alpha_{d-1}-1} (1-\sum^{d-1}_{j=1}x_j)^{\alpha_d-1} dx_{d-1}=\int_0^p x_{d-1}^{\alpha_{d-1}-1} (p-x_{d-1})^{\alpha_d-1} dx_{d-1}}}$$
$${\displaystyle{=(1-\sum^{d-2}_{j=1}x_j)^{\alpha_{d-1}+\alpha_d-1}B(\alpha_{d-1}, \alpha_d)}}$$
これを積分式に戻して、さらに$${p=1-\sum^{d-3}_{j=1}x_j}$$とすれば、
$${\displaystyle{ \int_0^{1-\sum^{d-3}_{j=1} x_{j}}x_{d-2}^{\alpha_{d-2}-1} (1-\sum^{d-2}_{j=1}x_j)^{\alpha_d+\alpha_{d-1}-1} dx_{d-2}}}$$
$${\displaystyle{=(1-\sum^{d-3}_{j=1}x_j)^{\alpha_{d-2}+\alpha_{d-1}+\alpha_d-1}B(\alpha_{d-2}, \alpha_{d-1}+\alpha_d)}}$$
これを順次適応して行って、
$${\displaystyle{B_d({\bf \alpha})=\Pi^{d-1}_{j=1} B(\alpha_j, \sum^{d}_{j'=j+1}\alpha_{j'}) }}$$が得られる。
ここに、ガンマ関数、
$${ B(\alpha,\beta)= \displaystyle{ \frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha+\beta)} }}$$を適用すれば、
$${\displaystyle{B_d({\bf \alpha})=\Pi^{d-1}_{j=1}\frac{\Gamma(\alpha_j)\Gamma(\sum^{d}_{j'=j+1}\alpha_{j'})}{\Gamma(\sum^{d}_{j'=j}\alpha_{j'})}=\frac{\Gamma(\alpha_1)\cdots\Gamma(\alpha_d)}{\Gamma(\sum^{d}_{j'=1}\alpha_{j'})}}}$$
$${\sum^{d}_{j'=1}\alpha_{j'}=\alpha_{tot}}$$として、最終的に、
$${\displaystyle{B_d({\bf \alpha})=\frac{\Pi_{j=1}^d\Gamma(\alpha_j)}{\Gamma(\alpha_{tot})}}}$$
が得られる。
明らかに、$${d=2}$$の時は、$${B_d({\bf \alpha})=B(\alpha_1,\alpha_2 )}$$のベータ関数になり、確率密度は、
$${f(x)=\displaystyle{\frac{x_1{^\alpha_1-1}(1-x_1)^{\alpha_2-1}}{B(\alpha_1,\alpha_2)}}}$$のベータ分布に従う確率変数と同じとなることから、ディリクレ分布はベータ分布の拡張と見ることができる。

ディレクレ分布に従う確率変数$${x_j}$$の期待値は、
$${E[x_j]=\displaystyle{\int_{x_j} dx_j \frac{x_j\Pi^d_{j'=1}x_{j'}^{\alpha_j'-1}}{B_d({\bf \alpha})} =\frac{B_d(\alpha_1, \cdots, \alpha_j+1, \alpha_{j+1},\cdots,\alpha_d)}{B_d(\alpha_1, \cdots, \alpha_d)} }}$$
$${=\displaystyle{\frac{\Gamma(\alpha_1)\cdots\Gamma(\alpha_j+1)\cdots\Gamma(\alpha_d)\cdot \Gamma^{-1}(\alpha_{tot}+1)}{\Gamma^{-1}(\alpha_{tot})\Pi_{j'=1}^d\Gamma(\alpha_j')}=\frac{\Gamma(\alpha_{tot})\Gamma(\alpha_j+1)}{\Gamma(\alpha_{tot}+1)\Gamma(\alpha_j)}}}$$
$${=\displaystyle{\frac{\alpha_j}{\alpha_{tot}} }}$$
最後に、$${\Gamma(\alpha+1)=\alpha\Gamma(\alpha)}$$を使っている。
分散共分散は以下の計算から得られる。
$${E[x_jx_j]=\displaystyle{\int_{x_j,x_i}dx_j dx_i \frac{x_ix_j\Pi^d_{j'=1}x_{j'}^{\alpha_j'-1}}{B_d({\bf \alpha})} =\frac{B_d(\alpha_1, \cdots,\alpha_i+(1-\delta_{ij}), \cdots, \alpha_j+1+\delta_{ij}, \cdots,\alpha_d)}{B_d(\alpha_1, \cdots, \alpha_d)} }}$$
$${=\displaystyle{\frac{\Gamma(\alpha_1)\cdots\Gamma(\alpha_i+1-\delta_{ij})\cdots \Gamma(\alpha_j+1+\delta_{ij})\cdots\Gamma(\alpha_d)\cdot \Gamma^{-1}(\alpha_{tot}+2)}{\Gamma^{-1}(\alpha_{tot})\Pi_{j'=1}^d\Gamma(\alpha_j')} }}$$
$${\displaystyle{ = \left\{\begin{array}{ll}\frac{\alpha_i\alpha_j}{\alpha_{tot}(\alpha_{tot}+1)} & (i \neq j)\\ \\ \frac{(\alpha_j+1)\alpha_j}{\alpha_{tot}(\alpha_{tot}+1)} & (i = j)\end{array}\right.}}$$
これから、
$${\displaystyle{ V[x_i, x_j] = \left\{\begin{array}{ll}\frac{\alpha_i\alpha_j}{\alpha_{tot}(\alpha_{tot}+1)}- \frac{\alpha_i\alpha_j}{\alpha_{tot}^2}& =-\frac{\alpha_i\alpha_j}{\alpha_{tot}^2(\alpha_{tot}+1)} & (i \neq j)\\ \\ \frac{\alpha_j(\alpha_j+1)}{\alpha_{tot}^2(\alpha_{tot}+1)}- \frac{\alpha_j^2}{\alpha_{tot}^2}&=\frac{\alpha_j(\alpha_{tot}-\alpha_j)}{\alpha_{tot}^2(\alpha_{tot}+1)} & (i =j)\end{array}\right.}}$$
となる。

ディレクレ分布に従う乱数のpythonでの実装は、numpy.randomのdirichletで以下のように行える。ここでは、$${{\bf \alpha}}$$を３次元のベクトル$${(10, 5, 3)}$$として、$${{\bf x}}$$20組を生成している。

import numpy as np
from numpy import random as rand
import matplotlib.pyplot as plt
rand.seed(0)

a=(10, 5, 3)
x = rand.dirichlet(alpha=a, size=20).transpose()
plt.barh(range(20), x[0])
plt.barh(range(20), x[1], left=x[0], color='g')
plt.barh(range(20), x[2], left=x[0]+x[1], color='r')
plt.title("Dirchlet sum of each x vector");