디리클레분포(Dirichlet distribution)

디리클레분포의 정의

$\mathbf{X} = (X_1,\dots,X_K) \sim Dirichlet(\alpha_1, \dots, \alpha_K)$일 때 $\mathbf{X}$는 $K-1$ simplex 위의 support(정의역) 

$$S_K = \{ \mathbf{x} : 0 \le x_i \le 1, \sum_{i=1}^{K}{x_i} = 1 \}$$

에서 정의되며, 확률밀도함수(pdf)는 다음과 같이 주어진다.

 

$$f(x_1,\dots,x_K) = \frac{\Gamma(\alpha_0)}{\Gamma(\alpha_1) \cdots \Gamma(\alpha_K)} \prod_{i=1}^{K}x_i^{\alpha_i-1}\mathrm{I}_{\mathbf{x} \in S_k}(\mathbf{x})$$

이다. (단, $\displaystyle \alpha_0 = \sum_{i=1}^{K}{\alpha_i}$이다.)

 

특히, $K=2$인 경우, $\mathbf{X}$는 베타분포 $Beta(\alpha_1, \alpha_2)$를 따르는 일변수 확률변수가 된다. (2변수 확률변수여야 할 것 같지만, $x_1+x_2=1$의 제약조건이 있으므로 1변수로 보는 것이 맞다.)

 

한편, $\displaystyle \int_{\mathbf{x} \in S_k} \prod_{i=1}^{K}x_i^{\alpha_i-1}{\rm d}\mathbf{x}   $의 적분값을 1로 만들어주는 normalizing constant $\displaystyle \frac{\Gamma(\alpha_0)}{\Gamma(\alpha_1) \cdots \Gamma(\alpha_K)} $은 다변수베타함수(multivariate beta function)을 이용하여 표현할 수 있다. 다변수베타함수는 다음과 같이 감마함수로 표현할 수 있다.

 

$${\rm B}(\alpha_1,\dots,\alpha_K) = \frac{\Gamma(\alpha_1)\cdots \Gamma(\alpha_K)}{\Gamma(\alpha_0)}, \quad \alpha_0 = \sum_{i=1}^{K}{\alpha_i} $$

 

따라서 디리클레분포의 pdf는

 

$$f(x_1,\dots,x_K) = \frac{1}{{\rm B}(\alpha_1,\dots,\alpha_K)  } \prod_{i=1}^{K}x_i^{\alpha_i-1}\mathrm{I}_{\mathbf{x} \in S_k}(\mathbf{x})$$

 

로 다시 쓸 수 있다.

 

즉, 이와 같은 디리클레분포의 정의로 미루어볼 때, 디리클레함수는 베타분포의 다변수 version이라고 할 수 있다.

 

 

 

디리클레분포의 성질

 

확률변수 $(X_1, \dots, X_{K})$가 각각 $Gamma(\alpha_1), \dots, Gamma(\alpha_{K})$을 따를 때, 1-1변환

$$ u(X_1, \cdots, X_K) = \begin{cases}  Y_1 &= \displaystyle \frac{X_1}{X_1+ \cdots + X_K} \\ &\; \vdots \\  Y_{K-1} &= \displaystyle \frac{X_{K-1}}{X_1+\cdots+X_K} \\ Y_K &= X_1+ \cdots + X_K \end{cases}$$

을 생각하자.

 

그 역변환 $u^{-1}$은

$$u^{-1}(y_1,\cdots,y_K) = \begin{cases} x_1 = y_1 y_K \\ \vdots \\ x_{K-1} = y_{K-1} y_K \\ x_K = y_K(1-y_1 - \cdots - y_{K-1})   \end{cases} $$

이고, 각 $X_i, i=1,\dots,K$는 positive random variable이므로, $Y_1,\dots,Y_K$의 support는 

$$S_Y = \{(y_1,\dots,y_K) | y_1>0, \dots, y_K>0, 0<y_1 + \cdots + y_{K-1}<1  \} $$

가 된다.

$u^{-1}$의 Jacobian 행렬식은

$$\begin{align*} J_{u^{-1}} &= \begin{vmatrix} y_K & 0    & \cdots & 0   & -y_K \\ 0    & y_K &           & 0   & y_K \\ \vdots &  & \ddots &     &   \vdots  \\ 0 & 0 & & y_K & -y_K  \\ y_1 & y_2 & \cdots & y_{K-1} & 1-\sum_{i=1}^{K-1}{y_i}\end{vmatrix} \\ \\  &= \begin{vmatrix} y_K & 0    & \cdots & 0   & 0 \\ 0    & y_K &           & 0   & 0 \\ \vdots &  & \ddots &     &   \vdots  \\ 0 & 0 & & y_K & 0  \\ y_1 & y_2 & \cdots & y_{K-1} & 1\end{vmatrix} = {y_K}^{K-1} \end{align*}$$

 

$(X_1, \dots, X_K)$의 joint pdf는

$$f(x_1,\dots,x_k) = \frac{1}{\Gamma(\alpha_1)\cdots\Gamma(\alpha_K)\beta^{\alpha_0}} \prod_{i=1}^{K}{x_i^{\alpha_i-1} } \exp{\left(-\frac{x_1+\dots+x_K}{\beta}\right)} {\rm I}_{x_1,\dots,x_K>0}(x_1,\dots,x_K)$$

따라서 $(Y_1,\dots,Y_K)$의 joint pdf는

$$\begin{align*} f(y_1,\dots,y_K) &= \frac{1}{\Gamma(\alpha_1)\cdots\Gamma(\alpha_K)\beta^{\alpha_0} }(y_1 y_K)^{\alpha_1-1}\cdots (y_{K-1} y_K)^{\alpha_{K-1}-1}\{y_K(1-\textstyle \sum_{i=1}^{K-1}{y_i}) \}^{\alpha_K-1} \displaystyle \exp{\left(-\frac{y_K}{\beta}\right)} \cdot J_{u^{-1}} \cdot {\rm I}_{S_Y}(y_1,\dots,y_K) \\ &= \frac{1}{\Gamma(\alpha_1)\cdots\Gamma(\alpha_K)\beta^{\alpha_0}} \prod_{i=1}^{K-1}{({y_i}^{\alpha_i-1})} {y_K}^{\alpha_0 - K} (1-\textstyle \sum_{i=1}^{K-1}{y_i})^{\alpha_K-1} \displaystyle \exp{\left(-\frac{y_K}{\beta}\right)} {y_K}^{K-1} {\rm I}_{S_Y}(y_1,\dots,y_K) \\ &= \frac{1}{\Gamma(\alpha_1)\cdots\Gamma(\alpha_K)\beta^{\alpha_0}} \prod_{i=1}^{K-1}{({y_i}^{\alpha_i-1})} {y_K}^{\alpha_0 - 1} (1-\textstyle \sum_{i=1}^{K-1}{y_i})^{\alpha_K-1} \displaystyle \exp{\left(-\frac{y_K}{\beta}\right)}{\rm I}_{S_Y}(y_1,\dots,y_K) \end{align*}$$

 

$(Y_1,\dots,Y_{K-1})$의 joint pdf를 계산하기 위해 $y_K$에 대해 적분하면

 

$$\begin{align*} f(y_1,\dots,y_{K-1}) &= \int_0^{\infty}{ f(y_1,\dots,y_K)  {\rm d}y_K  } \\ &= \int_0^{\infty}{ \frac{1}{\Gamma(\alpha_1)\cdots\Gamma(\alpha_K)\beta^{\alpha_0}} \prod_{i=1}^{K-1}{({y_i}^{\alpha_i-1})} {y_K}^{\alpha_0 - 1} (1-\textstyle \sum_{i=1}^{K-1}{y_i})^{\alpha_K-1} \displaystyle \exp{\left(-\frac{y_K}{\beta}\right)} {y_K}^{K-1} {\rm d}{y_K} } \\ &= \frac{\Gamma(\alpha_0)}{\Gamma(\alpha_1)\cdots\Gamma(\alpha_K)} \prod_{i=1}^{K-1}{ {y_i}^{\alpha_i-1} } (1-\textstyle \sum_{i=1}^{K-1}{y_i})^{\alpha_K-1} \displaystyle \int_0^{\infty}{ \frac{1}{\Gamma{(\alpha_0)}\beta^{\alpha_0}} {y_K}^{\alpha_0 - 1} \displaystyle \exp{\left(-\frac{y_K}{\beta}\right)} {\rm d}{y_K}} \\ &=  \frac{\Gamma(\alpha_0)}{\Gamma(\alpha_1)\cdots\Gamma(\alpha_K)} \prod_{i=1}^{K-1}{ {y_i}^{\alpha_i-1} } (1-\textstyle \sum_{i=1}^{K-1}{y_i})^{\alpha_K-1}\\ &= \frac{\Gamma(\alpha_0)}{\Gamma(\alpha_1)\cdots\Gamma(\alpha_K)} \prod_{i=1}^{K}{ {y_i}^{\alpha_i-1}}, \quad \text{where} \; y_K = 1-y_1-\cdots-y_{K-1} \end{align*}$$

 

이는 $Dirichlet(\alpha_1,\dots,\alpha_K)$를 따르는 pdf이다.

 

평균, 최빈값(mode), 분산, 공분산

 

디리클레분포의 각 $X_j, 1 \le j \le K$에 대해 평균, 최빈값, 분산, 공분산을 계산하면 다음과 같다.

 

평균

$$\begin{align*}  E(X_j) &= \int_{\mathbf{x} \in S_K}{ \frac{x_j}{{\rm B}(\alpha_1,\dots,\alpha_K)} \prod_{i=1}^{K}{ x_i^{\alpha_i -1} } {\rm d}\mathbf{x}} \\ &= \int_{\mathbf{x} \in S_K}{ \frac{1}{{\rm B}(\alpha_1,\dots,\alpha_K)} x_j^{\alpha_j} \prod_{1\le i \ne j \le K}{x_i^{\alpha_i -1}  } {\rm d}\mathbf{x}}  \\ &= \frac{{\rm B}(\alpha_1,\dots,\alpha_j +1 ,\dots, \alpha_K)}{{\rm B}(\alpha_1, \dots, \alpha_K)} \boxed{\int_{\mathbf{x} \in S_K}{ \frac{1}{{\rm B}(\alpha_1,\dots,\alpha_j +1 ,\dots, \alpha_K)} x_j^{\alpha_j} \prod_{1\le i \ne j \le K}{x_i^{\alpha_i -1}  }  {\rm d}\mathbf{x}} } \\ &= \frac{{\rm B}(\alpha_1,\dots,\alpha_j +1 ,\dots, \alpha_K)}{{\rm B}(\alpha_1, \dots, \alpha_K)} \\ &= \frac{\Gamma(\alpha_0)}{\Gamma(\alpha_1)\cdots\Gamma(\alpha_K)} \frac{\Gamma(\alpha_1)\cdots\Gamma(\alpha_j +1)\cdots\Gamma(\alpha_K)}{\Gamma(\alpha_0 + 1)} = \frac{\alpha_j}{\alpha_0} \end{align*}$$

 

여기서 박스 친 부분은 $Dirichlet(\alpha_1, \cdots, \alpha_j +1 ,\cdots, \alpha_K)$의 pdf의 적분임을 이용하여 계산하였다.

 

최빈값(mode)

최빈값을 계산하기 위해서는 pdf의 값이 극대인 $\mathbf{x} = (x_1, \cdots, x_K)$의 값을 계산하면 된다.

 

pdf가 곱함수의 형태이므로 pdf에 자연로그를 취해주는 것이 편하다. (MLE를 계산할 때의 논리와 비슷하다) 

 

즉, 최빈값 mode를 찾는 것은 다음 극대화 문제를 푸는 것과 같다.

 

$$\begin{align*} \underset{\mathbf{x} \in S_K}{argmax} \log{ f(x_1,\cdots,x_K)} &= \underset{\mathbf{x} \in S_K}{argmax}  \left[ \sum_{i=1}^{K}{(\alpha_1 -1) \log{x_i} }  + \log{\left(\frac{1}{{\rm B}(\alpha_1,\cdots,\alpha_K)}\right)} \right] \\ &=  \underset{\mathbf{x} \in S_K}{argmax} \left[ \sum_{i=1}^{K}{(\alpha_1 -1) \log{x_i} } \right] \quad \text{subject to} \, \sum_{i=1}^{K}{x_i} = 1 \end{align*}$$

 

따라서 라그랑주 승수법을 이용하여 다음과 같은 라그랑주 함수의 일계도함수를 풀면 된다.

 

$$\begin{align*} \mathcal{L}(\mathbf{x},\lambda) &= \sum_{i=1}^{K}{(\alpha_i-1) \log{x_i} } + \lambda(1 - x_1 - \cdots - x_K)   \end{align*}$$

 

$$\begin{align*} \text{FOC :} \; \frac{\partial \mathcal{L}}{\partial x_j} &= \frac{\alpha_j -1}{x_j} - \lambda = 0, \quad ^\forall j = 1, \cdots, K  \\ \frac{\partial \mathcal{L}}{\partial \lambda} &= 1 - x_1 - \cdots - x_K = 0 \end{align*} $$

 

위 일계조건을 풀어주면 $\displaystyle \lambda = \sum_{i=1}^{K}{\alpha_i} - K$ 이므로 $\displaystyle x_j = \frac{\alpha_j -1}{\lambda} = \frac{\alpha_j -1}{\alpha_0 - K} $ 가 mode가 된다.

 

마지막으로 위에서 구한 점에서 극대임을 보이기 위해 Hessian matrix를 계산한다.

 

$\displaystyle H(\log{f(x_1,\dots,x_K)}) = \mathrm{diag}\left(-\frac{\alpha_i-1}{x_i^2}\right)_{i=1,\dots,K}$가 negative definite일 때  log pdf가 극대값을 가지게되고, 그 필요충분조건은 $\alpha_i > 1 \, \text{for} ^\forall i = 1,\dots,K$ 이므로

 

$\displaystyle x_j = \frac{\alpha_j -1}{\alpha_0 - K}$는 $\alpha_i >1, i=1,\dots,K$일 때 최빈값이 된다. 어떤 $i$에 대해 $0<\alpha_i<1$이라면 확률밀도함수가 마치 $\displaystyle y=\frac{1}{x}$의 그래프처럼 무한대로 발산하는 구간이 생기므로 최빈값이 정의되지 않는다.

 

분산

$$\begin{align*} E(X_j^2) &=  \int_{\mathbf{x} \in S_K}{ \frac{1}{{\rm B}(\alpha_1, \dots, \alpha_K)} x_j^{\alpha_j +2 -1} \prod_{1 \le i \ne j \le K}{x_i^{\alpha_i -1}} {\rm d}\mathbf{x}}  \\ &= \frac{{\rm B}(\alpha_1, \dots, \alpha_j +2, \dots, \alpha_K)}{{\rm B}(\alpha_1, \dots, \alpha_K)} \boxed{\int_{\mathbf{x} \in S_K}{ \frac{1}{{\rm B}(\alpha_1,\dots, \alpha_j+2,\dots,\alpha_K)}x_j^{\alpha_j +2 -1} \prod_{1 \le i \ne j \le K}{x_i^{\alpha_i -1}}  {\rm d}\mathbf{x}}} \\ &= \frac{{\rm B}(\alpha_1, \dots, \alpha_j +2, \dots, \alpha_K)}{{\rm B}(\alpha_1, \dots, \alpha_K)} \\ &=  \frac{\Gamma(\alpha_1 + \cdots + \alpha_K)}{\Gamma(\alpha_1)\cdots\Gamma(\alpha_K)} \frac{\Gamma(\alpha_1)\cdots\Gamma(\alpha_j+2)\cdots\Gamma(\alpha_K)}{\Gamma(\alpha_1 + \cdots + \alpha_K + 2)} \\ &= \frac{\Gamma(\alpha_j+2)}{\Gamma(\alpha_j)} \frac{\Gamma(\alpha_1+\cdots+\Gamma_K)}{\Gamma(\alpha_1 + \cdots + \alpha_K + 2)}  \\ &= \frac{\alpha_j(\alpha_j+1)}{\alpha_0(\alpha_0+1)} \end{align*} $$

 

여기서 박스 친 부분은 $Dirichlet(\alpha_1, \dots,\alpha_j+2,\dots,\alpha_K)$의 pdf의 적분이므로 1이 된다. 따라서 $\alpha_0 = \sum_i{\alpha_i}$라 할 때,

 

$$\begin{align*} Var(X) &= E(X^2) - E(X)^2 \\ &= \frac{\alpha_j(\alpha_j+1)}{\alpha_0(\alpha_0+1)}- \frac{\alpha_j^2}{\alpha_0^2} \\ &= \frac{\alpha_j(\alpha_0 - \alpha_j)}{\alpha_0^2 (\alpha_0+1)}  \end{align*}$$

 

공분산

서로 다른 $i, j$에 대해

$$\begin{align*} E(X_i X_j) &= \int_{\mathbf{x} \in S_K}{ \frac{1}{{\rm B}(\alpha_1,\dots, \alpha_K)} x_i^{\alpha_i} x_j^{\alpha_j} \prod_{1 \le k \ne i, j \le K}{x_k^{\alpha_k-1}}  {\rm d}\mathbf{x}}  \\ &= \frac{{\rm B}(\alpha_1,\dots,\alpha_i+1, \dots, \alpha_j+1,\dots,\alpha_K)}{{\rm B}(\alpha_1, \dots, \alpha_K)}\int_{\mathbf{x} \in S_K}{\frac{1}{{\rm B}(\alpha_1,\dots,\alpha_i+1,\dots,\alpha_j+1,\dots,\alpha_K)} x_i^{\alpha_i} x_j^{\alpha_j} \prod_{1 \le k \ne i, j \le K}{x_k^{\alpha_k-1}} {\rm d}\mathbf{x}} \\ &= \frac{{\rm B}(\alpha_1, \dots, \alpha_i+1,\dots, \alpha_j+1,\dots,\alpha_K)}{{\rm B}(\alpha_1,\dots,\alpha_K)} \\ &= \frac{\Gamma(\alpha_1)\cdots\Gamma(\alpha_i+1)\cdots\Gamma(\alpha_j+1)\cdots\Gamma(\alpha_K)}{\Gamma(\alpha_0 + 2)} \frac{\Gamma(\alpha_0)}{\Gamma(\alpha_1) \cdots \Gamma(\alpha_K) } \\ &= \frac{\Gamma(\alpha_i+1)\Gamma(\alpha_j+1)}{\Gamma(\alpha_i)\Gamma(\alpha_j)} \frac{\Gamma(\alpha_0)}{\Gamma(\alpha_0+2)} \\ &= \frac{\alpha_i \alpha_j}{\alpha_0(\alpha_0+1)}\end{align*}$$

공분산 공식에 대입하면 다음과 같다.

$$\begin{align*} Cov(X_i,X_j) &= E(X_i X_j) - E(X_i)E(X_j) \\ &= \frac{\alpha_i \alpha_j}{\alpha_0(\alpha_0+1)} - \frac{\alpha_i \alpha_j}{\alpha_0^2} \\ &= -\frac{\alpha_i \alpha_j}{\alpha_0^2(\alpha_0+1)}  \end{align*}$$

 

사전확률벡터로서의 디리클레분포

 

지난 베타분포에 관한 포스팅에서 베타분포는 정의역이 $(0,1)$인 분포로서, 베이즈통계에서 확률을 모수로하는 확률변수의 사전분포로 쓰임을 설명했다.

 

마찬가지로, 디리클레분포는 베타분포의 다변수 확장 분포임을 생각하면, 이항분포의 multivariate version인 다항분포의 모수 확률벡터 $\mathbf{p} \in \mathbb{R}^K$의 사전분포로 쓰일 수 있음을 직관적으로 유추할 수 있다.

 

다항분포의 랜덤 표본 $\mathbf{X}_n = (X_1, \dots, X_K), X_1 + \cdots + X_K = n$이 주어졌을 때, $\mathbf{p} \in \mathbb{R}^K$의 사전분포를 $Dirichlet(\alpha_1, \dots, \alpha_K)$라 하면 다음과 같이 베이즈정리를 이용하여 사후확률밀도함수 $\pi(\mathbf{p}|x_1,\dots,x_K)$를 계산할 수 있다.

 

$\pi(\mathbf{p}|x_1,\dots,x_K) = pdf(x_1,\dots,x_K|\mathbf{p}) \pi(\mathbf{p})$이고, 사전확률 $\mathbf{p}$가 주어졌을 때 다항분포의 pdf는

$$f(x_1,\dots,x_K | \mathbf{p}) = \frac{n!}{x_1! \cdots x_n !}p_1^{x_1}\cdots p_K^{x_K}, x_1+\cdots+x_K=n$$

이고,

사전분포는 $\displaystyle  \pi(\mathbf{p}) = \frac{1}{{\rm B}(\alpha_1,\dots,\alpha_K)}\prod_{i=1}^{K}{x_i^{\alpha_i-1}}$ 이므로

 

$$\begin{align*} \pi(\mathbf{p}|x_1,\dots,x_n) &\propto \frac{n!}{x_1!\cdots x_K!} \prod_{i=1}^{K}{p_i ^{x_i}} \cdot \frac{1}{{\rm B}(\alpha_1,\dots,\alpha_K)} \prod_{i=1}^K{p_i^{\alpha_i-1}} \\ &\propto \boxed{\frac{1}{{\rm B}(x_1+\alpha_1,\cdots,x_K+\alpha_K)} \prod_{i=1}^{K} p_i^{x_i+\alpha_i-1}} \end{align*}$$

 

한편, $\mathbf{p}$의 support는 다항분포의 모수의 조건에 의해

$$P_K = \{ \mathbf{p} : 0 \le p_i \le 1, \sum_{i=1}^{K}{p_i} = 1 \}$$

로 정해지므로, 위 박스친 부분 사후분포는 $Dirichlet(x_1+\alpha_1, \dots, x_K+\alpha_K)$를 따르게 된다. 특히, prior distribution과 posterior distribution이 디리클레분포로 동일하므로 사전분포는 켤레 사전분포(conjugate prior distribution)에 해당하게 된다.