Dirichlet Process: Stick-Breaking Process

Posted on Edited on In
앞에서, 디리슐레 분포와 다항분포의 켤레관계에 대해서 살펴보았다. 디리슐레 분포는 고정된 k개의 클러스터 (선택지)에 대한 확률을 모델링하는 분포로 사전분포로 쓰였다.,디리슐레 과정 (DP; Dirichlet Process)은 디리슐레 분포의 무한대 확장판으로, 고정되지 않은 클러스터에 대한 확률을 모델링하는 분포이다. 명확한 정의는 다음과 같다.

Definition

α>0\alpha>0, HH를 SS에서 정의된 probability measure라 하자. 이 두 모수 (α,H\alpha, H)를 가지는 디리슐레 과정은 SS에서 정의된 random probability measure GG이다. Random probability measure라는 뜻은 GG가 유한한 집합 A=(A1,…,Ar)\mathbf{A} = (A_1,\ldots, A_r)에 확률 G(A)=(G(A1),…,G(Ar))G(\mathbf{A})= (G(A_1),\ldots, G(A_r))을 정해줌을 의미한다.

G∣α,H∼DP(α,H)⇔(G(A1),⋯ ,G(Ar))∣α,H∼Dirichlet(αH(A1),⋯ ,αH(Ar))\begin{aligned} G|\alpha, H &\sim \text{DP}(\alpha,H)\\ \Leftrightarrow (G(A_1),\cdots, G(A_r))|\alpha, H &\sim \text{Dirichlet}(\alpha H(A_1),\cdots ,\alpha H(A_r)) \end{aligned}

cf. 참고로 DP(α,H)\text{DP}(\alpha,H)대신에 DP(αH)\text{DP}(\alpha H)로 표기하기도 한다.

이를 설명하자면 다음과 같다.

  1. HH는 SS에 정의된 base probability measure로, 연속인 확률분포를 가진다.
  2. α>0\alpha>0는 concentration 모수라 불린다. 이는 디리슐레 과정에서 나온 표본들이 얼마나 base 함수 HH와 얼마나 가까운지를 조절하는 모수이다.

Properties of Dirichlet Process

  • 유한 집합 A=(A1,…,Ar)\mathbf{A}=(A_1,\ldots,A_r)는 exclusive & exhaustive한 성질을 가진 파티션이다. exclusive하고 exhaustive하다는 말은 아래의 사진처럼, A1,…,A6A_1,\ldots, A_6은 겹치지 않고 (A1∩A2∩⋯∩Ar=∅A_1\cap A_2\cap \cdots \cap A_r=\varnothing), 전체의 합집합이 전체집합 (A1∪A2∪⋯∪Ar=AA_1\cup A_2 \cup \cdots \cup A_r = \mathbf{A}))임을 의미한다.

파티션

  • DP에서는 무한대 확률 벡터를 rr개로 나누어 무한개의 이산형의 분포를 가진다. 따라서, DP를 다음과 같이 정의할 수 있다.

G(⋅)=∑r=1∞πrδθr(⋅)δθr(S)={1,θr∈S0,θr∉S \begin{aligned} G(\cdot) &= \sum^\infty_{r=1} \pi_r \delta_{\theta_r}(\cdot)\\ \delta_{\theta_r}(S)&= \begin{cases} 1, & \theta_r \in S\\ 0, & \theta_r \not\in S \end{cases} \end{aligned}

여기서 πr\pi_r는 rr번째 클러스터로 분류된 확률 질량 (probability mass)이고 θr\theta_r는 rr번째 클러스터의 평균, δ\delta는 디랙 델타 함수이다. 또한 ∑r=1∞πr=1\sum^\infty_{r=1} \pi_r = 1이다.

  • DP의 평균과 분산은 다음과 같다.

  • E[G(A)]=H(A)E[G(A)] = H(A)

  • V[G(A)]=H(A)(1−H(A))α+1V[G(A)] = \frac{H(A)(1-H(A))}{\alpha+1}

Stick-Breaking Process

DP와 짝꿍으로 나오는 설명은 항상 Stick-breaking process이다. 필자는 처음에 왜 이게 갑자기 나오나 하는 생각이 들었지만, KAIST Youtube 강의
에서 그 해답을 찾을 수 있었다.

DP는 무한대에서 정의되는 과정이기 때문에 이 과정에서 표본을 추출하기 위한 알고리즘 (sampling schemes/construction scheme)이 필요하다. 즉, 위의 DP의 정의에서 θr\theta_r과 πr\pi_r을 construct하기 위한 방법론이 필요하다. 우리가 살펴볼 방법은 세 가지이다.

  1. Chinese restaurant process
  2. Polya-Urn Scheme
  3. Stick-breaking process

그 중 Stick-breaking process를 먼저 쓰는 이유는 이해하기 가장 직관적이였기 때문이다.

Stick-Breaking Process

β1,β2,…∣α∼Beta(1,α)\beta_1,\beta_2,\ldots|\alpha \sim \text{Beta}(1,\alpha)를 따르는 무한 개의 β\beta가 있을 때, Stick-breaking process를 통해 πr=βr∏i=1r−1(1−βi)=βr(1−β1)(1−β2)⋯(1−βr−1)\pi_r = \beta_r \prod^{r-1}_{i=1} (1-\beta_i) = \beta_r (1-\beta_1) (1-\beta_2)\cdots (1-\beta_{r-1})를 만든다. βi\beta_i는 베타분포를 따르기 때문에 0과 1사이의 값는 비율로 해석할 수 있다.

방법을 설명하면 다음과 같다.

  1. 11의 길이를 가진 막대가 있다.
  2. 이 막대를 π1=β1\pi_1=\beta_1의 비율로 두 개로 자르면, 두 개는 각각 π1=β1\pi_1=\beta_1과 1−β11-\beta_1의 길이를 갖는다.
  3. β1\beta_1길이의 막대는 냅두고, 1−β11-\beta_1의 길이의 막대만 β2\beta_2의 비율로 나누면, π2=β2(1−β1)\pi_2 = \beta_2(1-\beta_1) 길이의 막대와 (1−β1)(1−β2)(1-\beta_1)(1-\beta_2) 길이의 막대로 나눌 수 있다.
  4. 이를 ∞\infty번 반복하면 다음의 식을 얻는다.

π1=β1π2=β2(1−β1)π3=β3(1−β1)(1−β2)⋮\begin{aligned} \pi_1 &= \beta_1 \\ \pi_2 &= \beta_2 (1-\beta_1)\\ \pi_3 &= \beta_3 (1-\beta_1) (1-\beta_2)\\ &\vdots \end{aligned}

만약 π1,π2,…∼H\pi_1,\pi_2,\ldots \sim H라면,

G=∑r=1∞πrδθr∼DP(α,H)G = \sum^\infty_{r=1} \pi_r \delta_{\theta_r} \sim \text{DP}(\alpha,H)

를 따른다.

Example

그렇다면, stick breaking process를 파이썬으로 구현해보자. 이는 Austin Rochford Blog을 참조하였다.

Step1. 필요한 모듈 정의

1
2
3
4
5
6
7
8
%matplotlib inline
from matplotlib import pyplot as plt
import numpy as np
import pymc3 as pm
import scipy
import seaborn as sns
from statsmodels.datasets import get_rdataset
from theano import tensor as T
  • 파이썬의 모듈이란 R 상의 패키지의 개념이다.
  • %matplotlib inline은 Jupyter notebook 내에서 그림을 출력하기 위한 magic command이다.
  • seaborn 모듈은 Matplotlib을 기반으로 다양한 색상 테마와 통계용 차트 등의 기능을 추가한 시각화 패키지이다.
  • get_rdataset 모듈은 R패키지에 내장되어 있는 데이터셋을 불러올 수 있는 모듈이다.

Step2. Stick-breaking process를 통해 β\beta, ω\omega, π\pi 표본추출

1
2
3
4
5
6
7
8
9
10
np.random.seed(433)
N=20            # n_obs
K=30            # n_sticks
alpha=2.        # parameter in Beta dist 
H = ss.norm     # base dist
beta = ss.beta.rvs(1,alpha, size=(N,K)) # Beta(1,alpha) dist
pi = np.empty_like(beta)
pi[:, 0] = beta[:,0]                    # pi1 = beta1
pi[:, 1:] = beta[:, 1:] * (1-beta[:, :-1]).cumprod(axis=1)
omega = H.rvs(size=(N,K))

이 과정이 Stick-breaking process를 구현하는 핵심 과정이다.

  • beta.rvs: 앞에서 봤듯이 β1,β2,…∼Beta(1,α)\beta_1,\beta_2,\ldots \sim \text{Beta}(1,\alpha)를 따른다. 이때 무한개의 β\beta를 생성할 수는 없으므로 NN개의 관측치를 가진 KK개의 β\beta를 생성한다.
  • np.empty_like로 0으로 채워진 N×KN\times K의 배열을 생성해 π\pi의 초기값으로 저장한다.
  • 이후, cumprod 메쏘드를 이용해 누적곱을 계산해준다. 이때 axis=1은 행의 축을 기준으로 곱해준다는 의미로, r−1r-1개의 β\beta를 곱해주는 역할을 한다.

π1=β1πr=βr⋅∏i=1r−1(1−βi), r=2,⋯ ,K\begin{aligned} \pi_1 &= \beta_1\\ \pi_r &= \beta_r \cdot \prod^{r-1}_{i=1}(1-\beta_i),\ r=2,\cdots,K \end{aligned}

  • omega는 base 분포인 HH에서 생성한 난수이다.

Step3. 플랏 그리기

1
2
3
4
5
6
7
8
9
10
11
x_plot = np.linspace(-3,3,200)            
sample_cdfs = (pi[..., np.newaxis]* np.less.outer(omega, x_plot)).sum(axis=1)

fig, ax = plt.subplots(figsize=(8,6))

ax.plot (x_plot, sample_cdfs[0],c="gray", alpha=0.75, label = "DP sample CDFs")
ax.plot(x_plot, sample_cdfs[1:].T, c="gray", alpha=0.75) #.T: transpose (20,200) => (200,20)
ax.plot(x_plot, H.cdf(x_plot), c= "k", label = "Base CDF")

ax.set_title(r'$\alpha = {}$'.format(alpha))
ax.legend(loc=2)

파이썬으로 그래프를 그릴 때는 matplotlib.pyplot을 자주 사용한다. 이를 위한 과정은 다음과 같다.

  • x_plot: 그래프를 그릴 때 필요한 xx의 grid를 그려준다.
  • pi[...,np.newaxis]: π\pi는 β\beta와 같은 차원인 (20,30)의 배열이었다. 이에 np.newaxis를 사용하면 2차원의 배열이 3차원 (20,30,1)로 차원이 하나 늘어난다. 이런 과정을 거친 이유는 np.less.outer(omega,x_plot)의 3차원의 배열과 곱해주기 위해서이다.
  • np.less.outer(omega,x_plot): 이 부분이 이해하기 가장 어려웠다. 기존의 np.outer(x,y)는 x=(x1,…,xn)x=(x_1, \ldots, x_n)와 y=(y1,…,ym)y = (y_1,\ldots,y_m)의 외적인 (n,m)행렬을 반환한다.

outer(x,y)=(x1y1x1y2⋯x1ymx2y1x2y2⋯x2ym⋮⋮…⋮xny1xny2⋯xnym)\mathtt{outer(x,y)} = \left(\begin{array}{cccc} x_1y_1 & x_1y_2 & \cdots & x_1 y_m\\ x_2y_1 & x_2y_2 & \cdots & x_2 y_m\\ \vdots & \vdots & \ldots & \vdots \\ x_ny_1 & x_ny_2 & \cdots & x_n y_m \end{array}\right)

그러나 np.less.outer(x,y)는 더이상 외적에 대한 내용이 아니다. 이는 xx와 yy의 원소들끼리 비교했을 때 xi<yjx_i < y_j이면 TRUE (1)를 반환하고, xi≥yjx_i \ge y_j이면 FALSE(0)을 반환한다. 여기서 omega는 (20,30)의 배열이고, x_plot은 (200,)의 배열이었기 때문에, omega의 원소 값이 x_plot의 원소 값보다 작으면 1을 반환하는 (20,30,200)의 차원을 계산해준다.

outer(x,y)=(1x1<y11x1<y2⋯1x1<ym1x2<y11x2<y2⋯1x2<ym⋮⋮…⋮1xn<y11xn<y2⋯1xn<ym)\mathtt{outer(x,y)} = \left(\begin{array}{cccc} \mathbf{1}_{x_1< y_1} & \mathbf{1}_{x_1 < y_2} & \cdots & \mathbf{1}_{x_1< y_m}\\ \mathbf{1}_{x_2< y_1} & \mathbf{1}_{x_2 < y_2} & \cdots & \mathbf{1}_{x_2< y_m}\\ \vdots & \vdots & \ldots & \vdots \\ \mathbf{1}_{x_n< y_1} & \mathbf{1}_{x_n< y_2} & \cdots & \mathbf{1}_{x_n< y_m} \end{array}\right)

즉, 이는 디랙 델타 함수 δθr\delta_{\theta_r}를 계산한 것이다. x_plot의 범위 내에 있으면 1이고, 범위 밖에 없으면 0의 값을 계산해준다.

  • sample_cdfs: 마지막으로 이를 axis=1 (행의 축)으로 더해줌으로써 G=∑r=1KπrδθrG = \sum_{r=1}^K\pi_r \delta_{\theta_r}를 계산해준다.

이를 α\alpha값에 따라 플랏을 그려주면,

과 같다. 즉, α\alpha가 커질수록 표본 분포들이 base 분포 HH로 집중된다. 이것이 α\alpha를 concentration 모수로 부르는 이유이다.

α→∞\alpha \rightarrow \infty, G→HG \rightarrow H

References