파이썬 EM 알고리즘 구현 예제(Gaussian Mixture Model)

Python 가우시안 혼합 모형(GMM) 구현(EM 알고리즘 예시)

파이썬에서 EM 알고리즘을 직접 구현한 예시를
대표적인 케이스인 가우시안 혼합 모형(GMM)의 경우를 통하여 살펴보도록 하겠습니다.

 

Step 1. 데이터 가정

두 개의 독립적인 정규분포에서 샘플을 10개씩 추출되어 총 20개의 데이터가 혼재하는
상황을 가정해보도록 하겠습니다.

여기서는 단순한 1차원 정규분포의 상황으로 예시를 살펴보겠습니다.

실제 데이터가 추출된 두 가우시안 분포는 다음과 같습니다.
분포 1 : 평균 = 0, 표준편차 = 2
분포 2 : 평균 = 4, 표준편차 = 1

Numpy 모듈을 통하여 위에서 가정한 분포를 따라 데이터들을 추출한 코드는 아래와 같습니다.

import numpy as np

data1 = np.random.randn(10) * 2 # mu = 0, sigma = 2
data2 = np.random.randn(10) + 4 # mu = 4, sigma = 1
data = np.append(data1, data2)

print(data)

# 데이터 추출 결과
[-0.665264 -0.53535   0.409343  4.297033 -0.748832  0.984836 -2.806239
  0.579069 -0.542261 -0.210497  4.007117  4.232208  1.873942  3.294335
  4.298381  3.880806  3.685262  4.185403  4.175795  3.923003]

이제 20개의 데이터가 어느 분포에서 추출되었는지와 각 분포의 평균과 표준편차를
EM 알고리즘을 통하여 예측해보는 과정을 구현해보겠습니다.

 

Step 2. E Step

EM 알고리즘의 E step은 초기 가정한 파라미터에 대하여 각 데이터가 해당 분포에서 추출될
likelihood를 계산하는 과정으로 구성됩니다.

초기 파라미터인 각 분포의 평균, 표준편차 및 가중치(alpha)를 아래와 같이 설정해보겠습니다.
(어느 정도 범위 내에서 다른 값으로 설정해도 무관합니다.)

# 초기 mu, sigma, alpha 가정
mu = [-1, 3]
sigma = [1, 1]
alpha = [0.5, 0.5]

 

반응형

 

이제 각 20개의 데이터가 분포 1, 분포 2에 포함될
상대적인 likelihood 값을 아래와 같이 계산해볼 수 있습니다.
정규분포의 pdf 값을 가져오는 과정은 scipy의 함수를 사용하였습니다.

import scipy as sp
import scipy.stats

# E Step

# 정규분포 객체 생성
rv1 = sp.stats.norm(loc = mu[0], scale = sigma[0])
rv2 = sp.stats.norm(loc = mu[1], scale = sigma[1])
rvs = [rv1, rv2]

W = np.zeros((2, len(data))) # likelihood를 저장할 행렬

# likelihood 행렬 계산
for i in range(2):
    for j in range(len(data)):
        W[i, j] = (rvs[i].pdf(data[j]) * alpha[i]) / ((rv1.pdf(data[j]) * alpha[0]) + (rv2.pdf(data[j]) * alpha[1]))

np.set_printoptions(precision=6, suppress=True) # 출력 양식 지정

print(W)

# likelihood 행렬 출력 결과(각 20개 데이터가 분포1, 분포2에서 추출되었을 likelihood)
[[0.999989 0.999977 0.995299 0.000002 0.999993 0.905729 1.       0.988236
  0.999978 0.999852 0.000006 0.000002 0.092602 0.000129 0.000002 0.00001
  0.000023 0.000003 0.000003 0.000008]
 [0.000011 0.000023 0.004701 0.999998 0.000007 0.094271 0.       0.011764
  0.000022 0.000148 0.999994 0.999998 0.907398 0.999871 0.999998 0.99999
  0.999977 0.999997 0.999997 0.999992]]

 

 

Step 3. M Step

EM 알고리즘의 후속 단계로 E Step에서 계산된 likelihood를 이용하여
파라미터들(각 분포의 평균, 표준편차 및 가중치)을 최적화하는 과정인 M Step이 진행됩니다.

평균, 표준편차 및 가중치 값들을 최적화하는 코드는 다음과 같습니다.

# M Step
for i in range(2):
    mu[i] = (W[i] @ data) / W.sum(axis = 1)[i]
    sigma[i] = np.sqrt((W[i] * ((data - mu[i]) ** 2)).sum() / W.sum(axis = 1)[i])
alpha = W.sum(axis = 1) / len(data)

업데이트된 파라미터 값들의 출력 결과는 아래와 같습니다.

print("mu : \n", mu)
print("sigma : \n", sigma)
print("alpha : \n", alpha)

# 출력 결과
mu : 
 [-1.397374924944736, 3.8720418289165672]
sigma : 
 [1.1152866047436278, 0.9275519014237806]
alpha : 
 [0.470922 0.529078]

이제, EM 알고리즘의 다음 iteration을 진행하고 싶다면 업데이트된 파라미터 값들을 가지고
위의 Step 2의 E 스텝과 Step 3의 M 스텝 과정을 반복해주시면 됩니다.