LaTeX new commands defined for this notebook: $ \newcommand{\bX}{\boldsymbol{X}} \newcommand{\bZ}{\boldsymbol{Z}} \newcommand{\bx}{\boldsymbol{x}} \newcommand{\bz}{\boldsymbol{z}} \newcommand{\btheta}{\boldsymbol{\theta}} \newcommand{\Q}{\mathcal{Q}} \newcommand{\E}{\mathbb{E}} \newcommand{\N}{\mathcal{N}} \newcommand{\bpi}{\boldsymbol{\pi}} \newcommand{\bmu}{\boldsymbol{\mu}} \newcommand{\bSigma}{\boldsymbol{\Sigma}} \newcommand{\L}{\mathcal{L}} \newcommand{\bI}{\boldsymbol{I}} \newcommand{\Mult}{\operatorname{Mult}} \newcommand{\cov}{\operatorname{cov}} \newcommand{\diag}{\operatorname{diag}} \newcommand{\Bern}{\operatorname{Bern}} $

\bX $\bX$, \bZ $\bZ$, \bx $\bx$, \bz $\bz$, \btheta $\btheta$, \Q $\Q$, \E $\E$, \N $\N$, \bpi $\bpi$, \bmu $\bmu$, \bSigma $\bSigma$, \bI $\bI$, \L $\L$, \Mult $\Mult$, \cov $\cov$, \diag $\diag$, \Bern $\Bern$

In [1]:
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
%matplotlib inline

mpl.style.use('seaborn-white')

9.3 An Alternative View of EM

  • 관측 데이터들의 집합을 $N\times D$ 행렬 $\bX$로, 잠재변수들의 집합을 $N\times K$ 행렬 $\bZ$로, 그리고 모수들의 집합을 벡터 $\btheta$로 표시하자. $$ \bX = \begin{pmatrix} \bx_1^T \\ \vdots \\ \bx_N^T \end{pmatrix}, \quad \bZ = \begin{pmatrix} \bz_1^T \\ \vdots \\ \bz_N^T \end{pmatrix} $$

  • EM 알고리즘의 목적은 잠재변수 $\bZ$를 가지는 모델 $p(\bX|\btheta)$의 모수 $\btheta$를 최대우도추정으로 구하는 것이다. 이때 로그우도함수는 $$ \ln p(\bX|\btheta) = \ln \Bigl\{ \sum_\bZ p(\bX,\bZ|\btheta) \Bigr\} \tag{9.29} $$ 를 사용하고 잠재변수 $\bZ$에 대한 총합(summation)이 로그 안에 들어있다는 사실 때문에 계산이 어렵고 복잡해진다.

  • 각 관측 데이터 $\bx_n$에 대응하는 잠재변수 $\bz_n$의 값이 데이터로 함께 주어지면 이 집합 $\{\bX,\bZ\}$를 완전(complete)데이터라 한다. (그림 9.9)

    Figure 9.9

    완전데이터의 최대우도추정은 로그우도함수 $\ln p(\bX,\bZ|\btheta)$를 사용한다.

  • 관측 데이터 $\bX$만 주어지면 불완전(incomplete)데이터라 한다. (그림 9.6)

    Figure 9.6

    불완전데이터인 경우 $\ln p(\bX,\bZ|\btheta)$를 바로 이용할 수 없으므로 대신 로그우도함수의 기댓값이 최대가 되는 모수 $\btheta^{new}$를 추정한다. (M step) $$ \E_\bZ\bigl[\ln p(\bX,\bZ|\btheta)\bigm| \bX\bigr]=\sum_\bZ p(\bZ|\bX)\ln p(\bX,\bZ|\btheta) $$ 그런데 이 기댓값을 구하려면 먼저 잠재변수의 사후 확률분포 $p(\bZ|\bX)$를 알아야 하므로, 이때 모수 $\btheta^{old}$가 필요하다. (E step)

The General EM Algorithm

Given a joint distribution $p(\bX,\bZ|\btheta)$ over observed variables $\bX$ and latent variables $\bZ$, governed by parameters $\btheta$, the goal is to maximize the likelihood function $p(\bX|\btheta)$ with respect to $\btheta$.

  1. Choose an initial setting for the parameters $\btheta^{old}$.

  2. (E step) Evaluate the posterior distribution of the latent variables given by $p(\bZ|\bX,\btheta^{old})$.

  3. (M step) Evaluate $\btheta^{new}$ given by $$ \btheta^{new} = \arg_\btheta\max\Q(\btheta,\btheta^{old}) \tag{9.32} $$ where the expectation of the complete-data log likelihood evaluated for some parameter $\btheta$ is denoted by $$ \Q(\btheta,\btheta^{old}) = \E_\bZ\bigl[\ln p(\bX,\bZ|\btheta)\bigm| \bX, \btheta^{old} \bigr] = \sum_{\bZ} p(\bZ|\bX,\btheta^{old}) \ln p(\bX,\bZ|\btheta). \tag{9.33} $$

  4. Check for convergence of either the log likelihood or the parameter values. If the convergence criterion is not satisfied, then let $$ \btheta^{old} \leftarrow \btheta^{new} \tag{9.34} $$ and return to Step 2.

9.3.1 Gaussian mixtures revisited

Complete data set

  • 가우시안 혼합모델에서 완전데이터 $\{\bX,\bZ\}$가 주어진 경우, 각 데이터 $\bx_n$에 해당하는 잠재변수 데이터 $\bz_n=(z_{n1},\dotsc,z_{nK})^T$은 $\bx_n$이 어떤 성분의 가우시안 확률분포 $\N(\bx|\bmu_k,\bSigma_k)$를 따르는지 알려주는 이진(binary) 값 $z_{nk}$들로 구성되어 있다. (1-of-$K$ coding scheme)

  • 성분 $k$에 속한 (즉 $z_{nk}=1$인) 데이터 $\bx_n$들로 이루어진 부분집합을 $\bX_k\subseteq\bX$라 하면, 가우시안 혼합모델의 우도함수는 $$ p(\bX,\bZ|\bpi,\bmu,\bSigma) = \prod_k \prod_{\bx\in\bX_k} \pi_k \N(\bx|\bmu_k,\bSigma_k) = \prod_n \prod_k \pi_k^{z_{nk}} \N(\bx_n|\bmu_k,\bSigma_k)^{z_{nk}} \tag{9.35} $$ 이고, 로그우도함수는 다음과 같다. $$ \ln p(\bX,\bZ|\bpi,\bmu,\bSigma) = \sum_k \Bigl[ |X_k|\ln\pi_k + \sum_{\bx\in\bX_k} \ln\N(\bx|\bmu_k,\bSigma_k) \Bigr] $$

  • (means $\bmu_k$ and covariances $\bSigma_k$) 평균과 공분산의 최대우도추정은 데이터가 부분집합 $\bX_k$로 주어졌을 때 가우시안 분포 $\N(\bx|\bmu_k,\bSigma_k)$에 대한 최대우도추정과 동일하다. 따라서 $N_k \equiv |\bX_k| = \sum_n z_{nk}$로 두면 \begin{align*} \bmu_k^{ML} &= \frac{1}{N_k} \sum_{\bx\in\bX_k} \bx = \frac{1}{N_k} \sum_n z_{nk} \bx_n \\ \bSigma_k^{ML} &= \frac{1}{N_k} \sum_{\bx\in\bX_k} (\bx-\bmu_k^{ML})(\bx-\bmu_k^{ML})^T = \frac{1}{N_k} \sum_n z_{nk} (\bx_n-\bmu_k^{ML})(\bx_n-\bmu_k^{ML})^T \end{align*}

  • (mixing coefficients $\pi_k$) 라그랑주 승수법을 이용하면 혼합계수의 최대우도추정은 $$ \pi_k^{ML} = \frac{N_k}{N} = \frac{1}{N} \sum_n z_{nk} \tag{9.37} $$

    Proof. 제한조건 $g(\bpi)=\sum_k \pi_k-1=0$을 만족시키는 로그우도함수 $f(\bpi)=\ln p(\bX,\bZ|\bpi,\bmu,\bSigma)=\sum_k N_k\ln\pi_k + \text{const}$의 최대를 구한다. 라그랑주 함수 $\L(\bpi,\lambda)=f(\bpi)+\lambda g(\bpi)$에 대해 $$ \nabla_{\bpi,\lambda}\L(\bpi,\lambda)=0 \iff \nabla_\bpi f = \lambda\nabla_\bpi g \quad\text{and}\quad g(\bpi)=0 $$ 이므로, 오른쪽 두 식을 만족시키는 $\lambda$를 먼저 찾으면 \begin{align*} \nabla_\bpi f = \lambda\nabla_\bpi &\:\Rightarrow\: \Bigl(\frac{N_1}{\pi_1},\dotsc,\frac{N_K}{\pi_K}\Bigr)^T = \lambda (1,\dotsc,1)^T \\ &\:\Rightarrow\: \lambda = \frac{N_1}{\pi_1} = \dotsb = \frac{N_K}{\pi_K} \end{align*} 이고, 제한조건 $\sum_k \pi_k=1$을 이용하면 $$ \lambda = \lambda \sum_k \pi_k = \sum_k (\lambda \pi_k) = \sum_k N_k = N \:\Rightarrow\: \pi_k=\frac{N_k}{\lambda}=\frac{N_k}{N} $$

Incomplete data set

가우시안 혼합모델에서 불완전데이터로 $\bX$만 주어진 경우, EM 알고리즘의 일반화를 적용해 보자.

  • 베이즈 법칙에 의해 $p(\bZ|\bX) p(\bX) = p(\bX|\bZ) p(\bZ)$이고, 가우시안 혼합모델이 \begin{align*} p(\bZ|\bpi) &= \prod_n \prod_k \pi_k^{z_{nk}} \tag{9.10} \\ p(\bX|\bZ,\bmu,\bSigma) &= \prod_n \prod_k \N(\bx_n|\bmu_k,\bSigma_k)^{z_{nk}} \tag{9.11} \end{align*} 이므로, 사후확률분포는 $$ p(\bZ|\bX,\bpi,\bmu,\bSigma) \propto p(\bX|\bZ,\bmu,\bSigma) p(\bZ|\bpi) = \prod_n \prod_k \bigl[ \pi_k \N(\bx_n|\bmu_k,\bSigma_k) \bigr]^{z_{nk}} \tag{9.38} $$ 또한 (by inspection of the directed graph in Figure 9.6 and making use of the d-separation criterion in 8.2.2) $$ p(\bZ|\bX,\bpi,\bmu,\bSigma) = \prod_l p(\bz_l|\bx_l,\bpi,\bmu,\bSigma) \tag{9.80} $$

  • (E step) 잠재변수의 사후확률 $p(\bZ|\bX)$에 대한 $z_{nk}$의 조건부 기댓값은 데이터 $\bx_n$에 대한 $k$번째 성분의 신뢰도 $\gamma(z_{nk})$와 같다. \begin{align*} \E_\bZ[z_{nk}|\bX] &= \sum_\bZ p(\bZ|\bX) z_{nk} = \sum_{\bz_1}\dotsb\sum_{\bz_N} z_{nk} p(\bz_1|\bx_1) \dotsb p(\bz_N|\bx_N) \\ &= \Bigl( \sum_{\bz_n} z_{nk} p(\bz_n|\bx_n) \Bigr) \cdot \prod_{l\neq n} \Bigl( \sum_{\bz_l} p(\bz_l|\bx_l) \Bigr) = \sum_{\bz_n} z_{nk} p(\bz_n|\bx_n) \\ &= \sum_{\bz_n} z_{nk} \frac{p(\bx_n|\bz_n)p(\bz_n)}{p(\bx_n)} = \frac{\sum_{\bz_n} z_{nk} \prod_i \bigl[ \pi_i \N(\bx_n|\bmu_i,\bSigma_i) \bigr]^{z_{ni}}}{\sum_{\bz_n} \prod_j \bigl[ \pi_j \N(\bx_n|\bmu_j,\bSigma_j) \bigr]^{z_{nj}}} \\ &= \frac{\pi_k \N(\bx_n|\bmu_k,\bSigma_k)}{\sum_j \pi_j \N(\bx_n|\bmu_j,\bSigma_j)} = \gamma(z_{nk}) \tag{9.39} \end{align*} 따라서 신뢰도 $\gamma(z_{nk})$를 계산하는데 모수 $\bpi^{old}$, $\bmu^{old}$, $\bSigma^{old}$가 필요하다.

  • (M step) 이제 완전데이터 로그우도함수의 기댓값을 구하면 다음과 같다. \begin{align*} \E_\bZ\bigl[ \ln p(\bX,\bZ|\bpi,\bmu,\bSigma) \bigm| \bX \bigr] &= \sum_\bZ p(\bZ|\bX) \ln p(\bX,\bZ|\bpi,\bmu,\bSigma) \\ &= \sum_\bZ p(\bZ|\bX) \sum_n \sum_k z_{nk} \bigl[ \ln\pi_k + \ln\N(\bx_n|\bmu_k,\bSigma_k) \bigr] \\ &= \sum_n \sum_k \Bigl\{ \sum_\bZ p(\bZ|\bX) z_{nk} \Bigr\} \cdot \bigl[ \ln\pi_k + \ln\N(\bx_n|\bmu_k,\bSigma_k) \bigr] \\ &= \sum_n \sum_k \E_\bZ[z_{nk}|\bX] \cdot \bigl[ \ln\pi_k + \ln\N(\bx_n|\bmu_k,\bSigma_k) \bigr] \\ &= \sum_n \sum_k \gamma(z_{nk}) \bigl[ \ln\pi_k + \ln\N(\bx_n|\bmu_k,\bSigma_k) \bigr] \tag{9.40} \end{align*} 이 기댓값을 최대로 만드는 모수는 9.2.2에서 구했던 최대우도추정과 동일하다. \begin{align*} \bmu_k^{new} &= \frac{1}{N_k} \sum_n \gamma(z_{nk}) \bx_n \tag{9.17} \\ \bSigma_k^{new} &= \frac{1}{N_k} \sum_n \gamma(z_{nk}) (\bx_n-\bmu_k^{new})(\bx_n-\bmu_k^{new})^T \tag{9.19} \\ \pi_k^{new} &= \frac{N_k}{N} \quad\text{where}\quad N_k=\sum_n\gamma(z_{nk}) \tag{9.22} \end{align*}

    Proof. (means $\bmu_k$) (9.40)을 $\bmu_k$에 대한 함수로 다시 기술하면 $$ \E_\bZ\bigl[ \ln p(\bX,\bZ|\bpi,\bmu,\bSigma) \bigm| \bX \bigr] = -\frac{1}{2} \sum_n \sum_j \gamma(z_{nj}) (\bx_n-\bmu_j)^T \bSigma_j^{-1} (\bx_n-\bmu_j) + \text{const} $$ 이므로 $\bmu_k$에 대해 미분한 후, 양변에 $\bSigma_k$를 곱한 다음 0이 되는 값을 구한다. $$ \sum_n \gamma(z_{nk}) \bSigma_k^{-1} (\bx_n-\bmu_k) = 0 \:\Rightarrow\: \sum_n \gamma(z_{nk}) (\bx_n-\bmu_k) = 0 \:\Rightarrow\: \bmu_k = \frac{\sum_n \gamma(z_{nk}) \bx_n}{\sum_n \gamma(z_{nk})} = \frac{\sum_n \gamma(z_{nk}) \bx_n}{N_k} $$

9.3.2 Relation to $K$-means

가우시안 혼합모델에서 $K$-평균 알고리즘은 클러스터 센터(평균) $\bmu_k$를 찾는데 반해 EM 알고리즘은 여기에 더해 혼합계수 $\pi_k$와 공분산 $\bSigma_k$까지 찾는다. 또한 $K$-평균 알고리즘이 각 데이터가 어느 클러스터에 속하는지 라벨로 (견고하게) 알려주는데 반해 EM 알고리즘은 신뢰도를 이용해 (느슨하게) 클러스터를 지정한다. 두 알고리즘은 유사한 모습을 가지고 있는데, 사실 $K$-평균 알고리즘은 EM 알고리즘의 특별한 (극한) 형태로 볼 수 있다.

  • 가우시안 혼합모델에서 모든 공분산을 $\bSigma_k=\epsilon\bI$로 두고 EM 알고리즘의 일반화를 적용해 보자.

  • (E step) 데이터 $\bx_n$에 대한 $k$번째 성분의 신뢰도 $\gamma(z_{nk})$의 극한은 클러스터 라벨 $r_{nk}$가 된다. \begin{align*} \gamma(z_{nk}) &= \frac{\pi_k \N(\bx_n|\bmu_k,\bSigma_k)}{\sum_j \pi_j \N(\bx_n|\bmu_j,\bSigma_j)} = \frac{\pi_k \exp\{-\|\bx_n-\bmu_k\|^2/2\epsilon\}}{\sum_j \pi_j \exp\{-\|\bx_n-\bmu_j\|^2/2\epsilon\}} \tag{9.42} \\ &\to r_{nk} = \begin{cases} 1 & \text{if $k=\arg\min_j\|\bx_n-\bmu_j\|^2$} \\ 0 & \text{otherwise} \end{cases} \tag{9.2} \quad\text{as}\quad \epsilon\to 0 \end{align*}

    Proof. $k_0=\arg\min_j\|\bx_n-\bmu_j\|^2$라 두면 \begin{align*} \frac{\exp\{-\|\bx_n-\bmu_j\|^2/2\epsilon\}}{\exp\{-\|\bx_n-\bmu_{k_0}\|^2/2\epsilon\}} &= \exp\Bigl\{\frac{\|\bx_n-\bmu_{k_0}\|^2-\|\bx_n-\bmu_j\|^2}{2\epsilon}\Bigr\} \\ &\to \begin{cases} 1 & \text{if $j=k_0$} \\ 0 & \text{otherwise} \end{cases} \quad\text{as}\quad \epsilon\to0 \end{align*} 따라서 (9.42)의 분자와 분모 모두를 $\exp\{-\|\bx_n-\bmu_{k_0}\|^2/2\epsilon\}$로 나누고 극한을 구하면 된다.

  • (M step) (9.17)에서 최대우도추정으로 구한 $\bmu_k^{new}$의 극한은 클러스터 센터가 된다. $$ \bmu_k^{new} \to \frac{1}{N_k}\sum_n r_{nk}\bx_n \quad\text{as}\quad \epsilon\to 0 \tag{9.4} $$

  • (distortion measure; inertia) 완전데이터 로그우도함수의 기댓값 (9.40)에 $\epsilon>0$을 곱하면 \begin{align*} \E_\bZ\bigl[\ln p(\bX,\bZ|\bpi,\bmu,\bSigma) \bigm| \bX \bigr] &= \sum_n \sum_k \gamma(z_{nk}) \ln \N(\bx_n|\bmu_k,\bSigma_k) + \text{const} \\ \epsilon \cdot \E_\bZ\bigl[\ln p(\bX,\bZ|\bpi,\bmu,\bSigma) \bigm| \bX \bigr] &= -\frac{1}{2} \sum_n \sum_k \gamma(z_{nk}) \| \bx_n-\bmu_k \|^2 + \epsilon\cdot\text{const} + \epsilon\ln\epsilon\cdot\text{const} \\ &\to -\frac{1}{2} \sum_n \sum_k r_{nk} \| \bx_n-\bmu_k \|^2 \quad\text{as}\quad \epsilon\to0 \tag{9.43} \end{align*} EM 알고리즘은 기댓값의 최대를 찾으므로 $\epsilon$을 곱해도 알고리즘이 작동하는데 전혀 영향을 받지 않는다. 따라서 $\epsilon\to0$이 되면 이 기댓값을 최대로 만드는 것은 $K$-평균 알고리즘에서 $J$의 최소를 구하는 것과 동일함을 알 수 있다.

9.3.3 Mixtures of Bernoulli distributions

Multivariate Bernoulli distribution

  • Consider a set of $D$ binary variables $x_i$, each of which is governed by a Bernoulli distribution with parameter $\mu_i$, so that $$ \Bern(\bx|\bmu) = \prod_i \mu_i^{x_i}(1-\mu_i)^{(1-x_i)} \tag{9.44} $$

  • The mean and covariance of this distribution are \begin{align*} \E[\bx] &= \bmu \tag{9.45} \\ \cov[\bx] &= \diag\{ \mu_i(1-\mu_i) \} \tag{9.46} \end{align*}

Mixtures of multivariate Bernoulli distributions

  • 베르누이 혼합모델의 확률질량함수는 $$ p(\bx|\bpi,\bmu) = \sum_k \pi_k \Bern(\bx|\bmu_k) \tag{9.47} $$ 이고, 평균과 공분산은 \begin{align*} \E[\bx] &= \sum_k \pi_k \bmu_k \tag{9.49} \\ \cov[\bx] &= \sum_k \pi_k \bigl\{ \bSigma_k + \bmu_k\bmu_k^T \bigr\} - \E[\bx]\E[\bx]^T \tag{9.50} \end{align*} where $\bSigma_k=\diag\{ \mu_{ki}(1-\mu_{ki}) \}$.

    공분산 $\cov[\bx]$는 대각행렬이 아니므로, 베르누이 혼합모델이 확률변수들 사이의 상관관계를 가지고 있음을 알 수 있다.

  • 관측 데이터 $\{\bx_1,\dotsc,\bx_N\}$들로 이루어진 행렬 $\bX$가 주어졌을 때 베르누이 혼합모델의 로그우도함수 $$ \ln p(\bX|\bpi,\bmu) = \sum_n \ln \Bigl\{ \sum_k \pi_k \Bern(\bx_n|\bmu_k) \Bigr\} \tag{9.51} $$ 는 로그 안에 총합을 포함하고 있으므로 단순한(closed) 형태의 최대우도추정 식을 얻을 수 없다.

EM algorithm of Bernoulli mixtures

  • 가우시안 혼합모델과 동일하게 $\sum_k z_k=1$인 $K$차원 이진변수 $\bz=(z_1,\dotsc,z_K)^T$를 잠재변수로 도입하면 사전확률분포와 조건부분포는 \begin{align*} p(\bz|\bpi) &= \prod_k \pi_k^{z_k} \tag{9.53} \\ p(\bx|\bz,\bmu) &= \prod_k \Bern(\bx|\bmu_k)^{z_k} \tag{9.52} \end{align*}

  • 결합분포를 주변화(marginalization) 하면 (9.47)과 동일한 결과를 얻는다. \begin{align*} p(\bx|\bpi,\bmu) &= \sum_\bz p(\bx,\bz|\bpi,\bmu) = \sum_\bz p(\bz|\bpi) p(\bx|\bz,\bmu) \\ &= \sum_\bz \prod_k \bigl[ \pi_k \Bern(\bx|\bmu_k) \bigr]^{z_k} = \sum_k \pi_k \Bern(\bx|\bmu_k) \end{align*}

  • (E step) 잠재변수의 사후확률 $p(\bZ|\bX)$에 대한 $z_{nk}$의 조건부 기댓값은 데이터 $\bx_n$에 대한 $k$번째 성분의 신뢰도 $\gamma(z_{nk})$와 같다. $$ \gamma(z_{nk}) \equiv \E_\bZ[z_{nk}|\bX] = \sum_{\bz_n} z_{nk} \frac{p(\bx_n|\bz_n)p(\bz_n)}{p(\bx_n)} = \frac{\sum_{\bz_n} z_{nk} \prod_i \bigl[ \pi_i \Bern(\bx_n|\bmu_i) \bigr]^{z_{ni}}}{\sum_{\bz_n} \prod_j \bigl[ \pi_j \Bern(\bx_n|\bmu_j) \bigr]^{z_{nj}}} = \frac{\pi_k \Bern(\bx_n|\bmu_k)}{\sum_j \pi_j \Bern(\bx_n|\bmu_j)} \tag{9.39} $$ 신뢰도 $\gamma(z_{nk})$를 계산하는데 모수 $\pi_k^{old}$와 $\bmu_k^{old}$가 필요하다.

  • (M step) 이제 완전데이터 로그우도함수의 기댓값을 구하면 다음과 같다. \begin{align*} \ln p(\bX,\bZ|\bpi,\bmu) = \sum_n \sum_k z_{nk} \Bigl\{ \ln\pi_k + \sum_i \bigl[ x_{ni}\ln\mu_{ki} + (1-x_{ni})\ln(1-\mu_{ki}) \bigr] \Bigr\} \tag{9.54} \\ \E_\bZ\bigl[ \ln p(\bX,\bZ|\bpi,\bmu) \bigm| \bX \bigr] = \sum_n \sum_k \gamma(z_{nk}) \Bigl\{ \ln\pi_k + \sum_i \bigl[ x_{ni}\ln\mu_{ki} + (1-x_{ni})\ln(1-\mu_{ki}) \bigr] \Bigr\} \tag{9.55} \end{align*} 이 기댓값을 최대로 만드는 모수는 가우시안 혼합모델과 동일하다. \begin{align*} \bmu_k^{new} &= \bar{\bx}_k = \frac{1}{N_k} \sum_n \gamma(z_{nk}) \bx_n \tag{9.58} \\ \pi_k^{new} &= \frac{N_k}{N} \quad\text{where}\quad N_k=\sum_n\gamma(z_{nk}) \tag{9.60} \end{align*}

  • 가우시안 혼합모델과 달리 베르누이 혼합모델에서는 우도함수가 무한대로 가는 특이점 문제가 발생하지 않는다. 베르누이 분포에서 확률질량은 $0\leq \Bern(\bx_n|\bmu_k)\leq 1$이므로 혼합모델의 우도함수 $p(\bx_n)=\sum_j \pi_j \Bern(\bx_n|\bmu_j)$도 상계가 정해지기 때문이다. 또한 우도함수가 0이 되는 특이점도 초기값에서 이러한 문제가 발생하지 않도록 설정하면 다시 나타나지 않는다. EM 알고리즘은 우도함수가 증가하는 방향으로 모수를 찾기 때문이다.

In [2]:
from sklearn.datasets import fetch_mldata
mnist = fetch_mldata('MNIST original')

# We now fit a data set of N = 600 such digits, comprising
# the digits ‘2’, ‘3’, and ‘4’.
N = 600
indices = np.random.choice(
    np.where((mnist.target >= 2) & (mnist.target <= 4))[0], size=N)
target = mnist.target[indices]
data = mnist.data[indices] / 255
# Here the digit images have been turned into binary vectors
# by setting all elements whose values exceed 0.5 to 1. and
# setting the remaining elements to 0.
mask = data > 0.5
data[mask] = 1.0
data[~mask] = 0.0

def draw_a_digit(ax, data, title):
    ax.imshow(data.reshape((28, 28)))
    ax.set_title(title)
    ax.tick_params(labelbottom=False, labelleft=False)

print('Number of samples = {}'.format(
    [(k, np.count_nonzero(target == k)) for k in (2, 3, 4)] ))

nrows = 3
ncols = 10
fig, ax = plt.subplots(nrows, ncols, figsize=(10,nrows*1.2))
ax = ax.ravel()
for i in range(nrows*ncols):
    n = np.random.randint(N)
    draw_a_digit(ax[i], data[n], str(int(target[n])))

Number of samples = [(2, 199), (3, 222), (4, 179)]
In [3]:
from scipy.stats import uniform

K = 3
D = data.shape[1] # 28 by 28

# The mixing coefficients were initialized to $\pi_k = 1/K$, and the parameters
# $\mu_{kj}$ were set to random values chosen uniformly in the range (0.25, 0.75)
# and then normalized to satisfy the constraint that $\sum_j \mu_{kj}=1$.
pi_init = np.repeat(1/K, K) # shape (K,)
x = uniform.rvs(loc=0.25, scale=0.5, size=(D, K))
mu_init = (x / np.sum(x, axis=0)).T # shape (K, D)

def draw_components(x, text=None):
    fig, ax = plt.subplots(ncols=K, figsize=(K*2,2))
    for k in range(K):
        draw_a_digit(ax[k], x[k], '$k={}$'.format(k))
    ax[-1].text(1.05, 0.9, text, fontsize=14, transform=ax[-1].transAxes)

print('Initial pi = {}'.format(pi_init))
draw_components(mu_init, text='Initialization')

Initial pi = [0.33333333 0.33333333 0.33333333]
In [4]:
from scipy.special import xlogy, xlog1py, logsumexp

# x=(N, D), pi=(K,), mu=(K, D), resp=(N, K)

# much faster than scipy.stats.bernoulli
def log_multi_bernoulli(x, mu):
    assert x.shape[-1] == mu.shape[-1]
    return np.sum(xlogy(x[:,np.newaxis], mu) +
                  xlog1py(1 - x[:,np.newaxis], -mu), axis=-1) # shape (N, K)

def E_step(x, pi, mu):
    assert pi.shape[0] == mu.shape[0] 
    z = np.log(pi) + log_multi_bernoulli(x, mu) # shape (N, K)
    return np.exp(z - logsumexp(z, axis=1)[:,np.newaxis]) # shape (N, K)

def M_step(x, resp):
    assert x.shape[0] == resp.shape[0]
    N_k = np.sum(resp, axis=0)
    pi = N_k / x.shape[0]
    mu = np.sum(resp[:,:,np.newaxis] * x[:,np.newaxis,:], axis=0) / N_k[:,np.newaxis]
    return pi, mu

mu = mu_init
pi = pi_init

n_times = 10
for i in range(n_times):
    resp = E_step(data, pi, mu)
    pi, mu = M_step(data, resp)
    draw_components(mu, text='Iteration {}'.format(i+1))

fig, ax = plt.subplots()
draw_a_digit(ax, mu.mean(axis=0), 'mean after {} interations'.format(n_times))

Figure 9.10 Illustration of the Bernoulli mixture model in which the top row shows examples from the digits data set after converting the pixel values from grey scale to binary using a threshold of 0.5. On the bottom row the first three images show the parameters $\mu_{ki}$ for each of the three components in the mixture model. As a comparison, we also fit the same data set using a single multivariate Bernoulli distribution, again using maximum likelihood. This amounts to simply averaging the counts in each pixel and is shown by the right-most image on the bottom row.