EM(expectation maximization) 알고리즘 python 구현

프로그래밍/python

EM(expectation maximization) 알고리즘 python 구현

ksyoon 2024. 9. 19. 14:06

em.py

0.00MB

EM 알고리즘은 다음과 같은 목표 함수를 최소화 하는 것이다.

$$J= {\sum_{n=1}^N \sum_{k=1}^K r_{nk} \Vert x_n - u_k \Vert^2}$$

여기서 $N$은 데이터의 갯수 $K$는 집단의 갯수이다.

$ r_{nk}$는 데이터 포인터 $x_n$대한 이진 표시 변수 $r_{nk} \in {0,1}$을 도입한 것이다.

목표는 $ J$를 최소화하는 $\{r_{nk}\}$와 $\{u_{k}\}$ 를 찾는 것이다.

이를 위하여 반복적 과정을 통해 $r_{nk}$ 와 $u_{k}$의 최적화를 수행한다.

일단 $u_{k}$ 에 대한 초기값을 설정한다.

첫단계 $u_{k}$ 를 고정한 채로 $ J$를 최소화하는 $r_{nk}$ 를 찾는다.

두번째 단계 $r_{nk}$ 를 고정하고 $ J$를 최소화하는 $u_{k}$ 를 찾는다.

이런 방식으로 값이 수렴할때 까지 이 두 단계의 최적화 과정을 반복한다.

$r_{nk}$을 업데이트 하는 과정을 E-step ( 현재 파라미터 값(이 경우 평균과 같은 파라미터)을 바탕으로 잠재 변수의 기대값을 계산 )

$u_{k}$을 업데이트 하는 과정을 M-step ( E-step에서 계산된 기대값을 사용해, 모델의 파라미터(평균, 분산 등)를 업데이트. 즉, 데이터를 설명하는 데 최적의 파라미터를 찾기 위해 평균을 새롭게 추정하는 과정이 M-step에서 발생)이라고 한다.

먼저 각각의 $n$에 대하여 $ \Vert x_n - u_k \Vert^2$ 값을 최소화 하는 $k$값에 대해 $r_{nk}$ 값을 1로 설정한다.

$$ r_{nk} = \begin{cases} 1 \quad if \; k=arg \; min_j \Vert x_n - u_j \Vert^2 \\ 0 \quad otherwise.\end{cases}$$

두번째로 $r_{nk}$ 를 고정한 채로 $u_{k}$ 값을 구할 때 $ J$ 는 $u_{k}$ 에 대한 제곱 함수이며 $u_{k}$ 에 대한 미분 값을 구하면

$$2\sum_{n=1}^N r_{nk} (x_n -u_k)=0$$

이를 $u_{k}$ 에 대하여 풀면

$$u_{k}= \cfrac {\sum_n r_{nk}x_n}{\sum_n r_{nk}}$$

이를 python으로 chatgpt를 이용하여 구현하면 아래와 같다.

아래 코드는 $r_{nk}$ 가 한 집단에 속하는지 안하는지 1과 0의 이진수가 아닌, 각 집단에 속하는 확률 밀도 함수를 이용하여 계산 하였다는 차이점이 있다.

# -*- coding: utf-8 -*-
"""
Created on Thu Sep 12 19:17:29 2024
@author: headway
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal
# 데이터 생성 및 초기화
np.random.seed(0)
data = np.random.randn(20, 2) # 20개의 2차원 데이터 생성
# 초기 평균, 공분산, 가중치 설정
means = np.array([[-1, -1], [1, 1], [0, 0]]) # 초기 클러스터 평균 (3개)
covariances = [np.eye(2), np.eye(2), np.eye(2)] # 초기 클러스터 공분산 (3개)
weights = [1/3, 1/3, 1/3] # 초기 클러스터 가중치 (3개)
# E-스텝: 각 데이터 포인트가 각 클러스터에 속할 확률 계산
def e_step(data, means, covariances, weights):
r = np.zeros((len(data), len(means)))
for k in range(len(means)):
r[:, k] = weights[k] * multivariate_normal.pdf(data, means[k], covariances[k])
r /= r.sum(axis=1, keepdims=True)
return r
# M-스텝: 새로운 평균, 공분산, 가중치를 계산
def m_step(data, r):
N_k = r.sum(axis=0)
weights = N_k / len(data)
means = (r.T @ data) / N_k[:, None]
covariances = []
for k in range(len(means)):
diff = data - means[k]
cov = (r[:, k][:, None] * diff).T @ diff / N_k[k]
covariances.append(cov)
return means, covariances, weights
# EM 알고리즘 실행: 10번 반복
for iteration in range(10):
r = e_step(data, means, covariances, weights)
means, covariances, weights = m_step(data, r)
# 클러스터 할당: 가장 높은 확률을 가지는 클러스터로 할당
cluster_assignments = np.argmax(r, axis=1)
# 데이터 시각화: 각 클러스터를 다른 색상으로 표시
plt.figure(figsize=(8, 6))
colors = ['red', 'blue', 'green']
for k in range(3):
cluster_data = data[cluster_assignments == k]
plt.scatter(cluster_data[:, 0], cluster_data[:, 1], color=colors[k], label=f'Cluster {k+1}')
plt.scatter(means[:, 0], means[:, 1], color='black', marker='X', s=200, label='Centroids')
plt.title('EM Algorithm Clustering with 3 Clusters')
plt.xlabel('Feature 1')
plt.ylabel('Feature 2')
plt.legend()
plt.grid(True)
plt.show()

결과는 위와 같다.