[데이터처리와분석] GMM EM Algorithm

 

 

 

[데이터처리와분석] GMM EM Algorithm

요약¶

1. unsupervised classification 문제를 풀 때, 각 포인트를 특정 가우시안 분포에 속하도록 만들고 싶다. 이것이 GMM
2. 처음에는 특정 가우시안 분포(평균, 분산)을 랜덤으로 잡고 각 포인트에 레이블을 준다 (초기화)
3. 레이블 된 포인트으로 다시 가우신안 분포의 평균, 분산을 구해서 다시 레이블링한다. (반복, E step, M step)
4. 그러면 결국 레이블을 잘 구분할 수 있는 것이다. -> 수렴 영상은 참고자료에서 볼 수 있다.

========

이 코드는 기대 최대화(Expectation-Maximization, EM) 알고리즘을 사용하여 가우시안 혼합 모델(Gaussian Mixture Model, GMM)의 매개변수를 추정하는 함수입니다. 각 단계별로 코드를 분석하겠습니다.

함수 서명¶

def em_gmm_orig(xs, pis, mus, sigmas, tol=0.01, max_iter=100):

• xs: 입력 데이터 포인트, n개의 샘플과 p개의 특성을 가진 numpy 배열.
• pis: 초기 혼합 계수 (mixture coefficients), 각 가우시안 분포의 가중치.
• mus: 초기 가우시안 분포의 평균값.
• sigmas: 초기 가우시안 분포의 공분산 행렬.
• tol: 수렴 허용 오차, 로그 가능도 함수의 변화가 이 값보다 작아지면 수렴으로 간주.
• max_iter: 최대 반복 횟수.

초기 설정¶

n, p = xs.shape
    k = len(pis)

    ll_old = 0

• n, p: 데이터 포인트 수와 특성 수.
• k: 가우시안 분포의 수.
• ll_old: 이전 반복의 로그 가능도.

반복 루프¶

for i in range(max_iter):
        exp_A = []
        exp_B = []
        ll_new = 0

• 최대 max_iter만큼 반복.
• exp_A, exp_B: 예상(E-step) 계산에 사용될 리스트(현재 사용되지 않음).
• ll_new: 현재 반복의 로그 가능도.

E-step¶

ws = np.zeros((k, n))
        for j in range(len(mus)):
            for i in range(n):
                ws[j, i] = pis[j] * mvn(mus[j], sigmas[j]).pdf(xs[i])
        ws /= ws.sum(0)

• j 는 GMM의 갯수를, i는 각 포인트를 의미하며, GMM에 따라 모든 포인트의 가중치에 따라 확률을 계산한다.
• mvn(mus[j], sigmas[j]).pdf(xs[i])는 다변량 가우시안 분포(Multivariate Gaussian Distribution)의 확률 밀도 함수를 계산하는 함수 호출입니다. 이를 통해 특정 데이터 포인트 xs[i]가 주어진 가우시안 분포에서 나올 확률을 계산합니다.
• ws: 책임도 행렬, 각 데이터 포인트가 각 가우시안 분포에 속할 확률.
• 각 데이터 포인트에 대해 가우시안 분포의 확률 밀도 함수를 계산하고, 혼합 계수로 가중치를 곱한 값을 계산.
• 각 데이터 포인트에 대한 책임도를 정규화하여 합이 1이 되도록 함.

M-step¶

혼합 계수 갱신¶

pis = np.zeros(k)
        for j in range(len(mus)):
            for i in range(n):
                pis[j] += ws[j, i]
        pis /= n

• 각 가우시안 분포의 혼합 계수를 갱신.
• 책임도의 합을 계산하여 새로운 혼합 계수로 설정.

평균 갱신¶

mus = np.zeros((k, p))
        for j in range(k):
            for i in range(n):
                mus[j] += ws[j, i] * xs[i]
            mus[j] /= ws[j, :].sum()

• 각 가우시안 분포의 평균을 갱신.
• 책임도로 가중평균을 계산하여 새로운 평균으로 설정.

공분산 행렬 갱신¶

sigmas = np.zeros((k, p, p))
        for j in range(k):
            for i in range(n):
                ys = np.reshape(xs[i] - mus[j], (2,1))
                sigmas[j] += ws[j, i] * np.dot(ys, ys.T)
            sigmas[j] /= ws[j,:].sum()

• 각 가우시안 분포의 공분산 행렬을 갱신.
• 책임도로 가중된 공분산을 계산하여 새로운 공분산 행렬로 설정.

로그 가능도 갱신 및 수렴 조건 확인¶

ll_new = 0.0
        for i in range(n):
            s = 0
            for j in range(k):
                s += pis[j] * mvn(mus[j], sigmas[j]).pdf(xs[i])
            ll_new += np.log(s)

        if np.abs(ll_new - ll_old) < tol:
            break
        ll_old = ll_new
        print(ll_new)

• 새로운 로그 가능도를 계산.
• 이전 반복의 로그 가능도와의 차이가 tol보다 작으면 수렴으로 간주하고 반복을 종료.
• 로그 가능도의 변화를 출력하여 수렴 과정을 확인.

최종 결과 반환¶

return ll_new, pis, mus, sigmas

• 최종 로그 가능도, 혼합 계수, 평균, 공분산 행렬을 반환.

종합 분석¶

이 함수는 GMM의 매개변수를 추정하기 위해 EM 알고리즘을 사용합니다. 각 반복에서 책임도를 계산하고, 이를 바탕으로 혼합 계수, 평균, 공분산 행렬을 갱신합니다. 수렴 조건이 만족되면 반복을 종료하고 최종 매개변수를 반환합니다. EM 알고리즘을 통해 데이터를 가우시안 혼합 모델로 잘 설명할 수 있게 됩니다.

======== 참고할만한 자료

1. 공돌이의 수학노트 - https://angeloyeo.github.io/2021/02/08/GMM_and_EM.html
   • GMM와 EM에 대해 너무 잘 설명한 내용이다 + 영상

 

 

GMM EM Algorithm¶

 

 

We learn how to perform Expectation Maximization to perform inference and learning the GMM.

EM algorithm consists of E-step and M-step

In E-step, we compute expectation based on the current guess of parameters.

In M-step, we use the expectation and optimize the parameter.

 

 

이번 시간에는 EM 알고리즘을 직접 구현해보도록 한다.

 

In [ ]:

from scipy.stats import multivariate_normal as mvn

 

 

위 함수를 통해 gaussian pdf 를 공식으로 직접 작성하는 것을 생략한다.

 

In [ ]:

def em_gmm_orig(xs, pis, mus, sigmas, tol=0.01, max_iter=100):

    n, p = xs.shape
    k = len(pis)

    ll_old = 0
    for i in range(max_iter):
        exp_A = []
        exp_B = []
        ll_new = 0

        # E-step
        ws = np.zeros((k, n))
        for j in range(len(mus)):
            for i in range(n):
                ws[j, i] = pis[j] * mvn(mus[j], sigmas[j]).pdf(xs[i])
        ws /= ws.sum(0)

        # M-step
        pis = np.zeros(k)
        for j in range(len(mus)):
            for i in range(n):
                pis[j] += ws[j, i]
        pis /= n

        mus = np.zeros((k, p))
        for j in range(k):
            for i in range(n):
                mus[j] += ws[j, i] * xs[i]
            mus[j] /= ws[j, :].sum()

        sigmas = np.zeros((k, p, p))
        for j in range(k):
            for i in range(n):
                ys = np.reshape(xs[i]- mus[j], (2,1))
                sigmas[j] += ws[j, i] * np.dot(ys, ys.T)
            sigmas[j] /= ws[j,:].sum()

        # update complete log likelihoood
        ll_new = 0.0
        for i in range(n):
            s = 0
            for j in range(k):
                s += pis[j] * mvn(mus[j], sigmas[j]).pdf(xs[i])
            ll_new += np.log(s)

        if np.abs(ll_new - ll_old) < tol:
            break
        ll_old = ll_new
        print(ll_new)

    return ll_new, pis, mus, sigmas

 

 

Let's generate dataset based on GMM

np.random.multivariate_normla 함수를 사용해보자.

 

In [ ]:

import numpy as np
np.random.seed(123)


# create data set
n = 1000  # number of dataset 
_mus = np.array([[0,4], [-2,0]])
_sigmas = np.array([[[3, 0], [0, 0.5]], [[1,0],[0,2]]])

ci = np.random.binomial(1, 0.6, n)    # 0.6의 확률로 1, 0.4의 확률로 0   1000 개를 샘플링 
 
xs = []
for i in range(n):
    cluster_index = ci[i]
    mu = _mus[cluster_index]
    sigma = _sigmas[cluster_index]
    xs.append( np.random.multivariate_normal(mu, sigma ))
    
xs = np.array(xs)

 

 

그래프로 먼저 확인해보자.

 

In [ ]:

%matplotlib inline
import matplotlib.pyplot as plt
plt.scatter(xs[:,0],xs[:,1])

 

Out[ ]:

<matplotlib.collections.PathCollection at 0x20541b22c50>

 

 

 

Assume we don't know the parameters

 

 

Initial guessing

 

In [ ]:

# initial guesses for parameters
pis = np.random.random(2)
pis /= pis.sum()
mus = np.random.random((2,2))
sigmas = np.array([np.eye(2)] * 2)

 

In [ ]:

ll2, pis2, mus2, sigmas2 = em_gmm_orig(xs, pis, mus, sigmas)

 

 

-3961.2210365915125
-3946.892674616119
-3928.560775439659
-3909.4275385395467
-3890.101351406793
-3869.06655723916
-3843.4506957557232
-3809.528520183506
-3766.390706837679
-3725.545861794156
-3701.4657368664407
-3690.0882743474376
-3684.312733498831
-3681.33687615515
-3679.8858958123888
-3679.2136261388578
-3678.9115798572
-3678.778057158327
-3678.719523216521
-3678.693975748356
-3678.6828523308864

 

 

Compare the parameters

 

In [ ]: