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

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의 이진수가 아닌, 각 집단에 속하는 확률 밀도 함수를 이용하여 계산 하였다는 차이점이 있다.

1. # -*- coding: utf-8 -*-
2. """
3. Created on Thu Sep 12 19:17:29 2024
4.  
5. @author: headway
6. """
7.  
8. import numpy as np
9. import matplotlib.pyplot as plt
10. from scipy.stats import multivariate_normal
11.  
12. # 데이터 생성 및 초기화
13. np.random.seed(0)
14. data = np.random.randn(20, 2)  # 20개의 2차원 데이터 생성
15.  
16. # 초기 평균, 공분산, 가중치 설정
17. means = np.array([[-1, -1], [1, 1], [0, 0]])  # 초기 클러스터 평균 (3개)
18. covariances = [np.eye(2), np.eye(2), np.eye(2)]  # 초기 클러스터 공분산 (3개)
19. weights = [1/3, 1/3, 1/3]  # 초기 클러스터 가중치 (3개)
20.  
21. # E-스텝: 각 데이터 포인트가 각 클러스터에 속할 확률 계산
22. def e_step(data, means, covariances, weights):
23.     r = np.zeros((len(data), len(means)))
24.     for k in range(len(means)):
25.         r[:, k] = weights[k] * multivariate_normal.pdf(data, means[k], covariances[k])
26.     r /= r.sum(axis=1, keepdims=True)
27.     return r
28.  
29. # M-스텝: 새로운 평균, 공분산, 가중치를 계산
30. def m_step(data, r):
31.     N_k = r.sum(axis=0)
32.     weights = N_k / len(data)
33.     means = (r.T @ data) / N_k[:, None]
34.  
35.     covariances = []
36.     for k in range(len(means)):
37.         diff = data - means[k]
38.         cov = (r[:, k][:, None] * diff).T @ diff / N_k[k]
39.         covariances.append(cov)
40.    
41.     return means, covariances, weights
42.  
43. # EM 알고리즘 실행: 10번 반복
44. for iteration in range(10):
45.     r = e_step(data, means, covariances, weights)
46.     means, covariances, weights = m_step(data, r)
47.  
48. # 클러스터 할당: 가장 높은 확률을 가지는 클러스터로 할당
49. cluster_assignments = np.argmax(r, axis=1)
50.  
51. # 데이터 시각화: 각 클러스터를 다른 색상으로 표시
52. plt.figure(figsize=(8, 6))
53. colors = ['red', 'blue', 'green']
54. for k in range(3):
55.     cluster_data = data[cluster_assignments == k]
56.     plt.scatter(cluster_data[:, 0], cluster_data[:, 1], color=colors[k], label=f'Cluster {k+1}')
57. plt.scatter(means[:, 0], means[:, 1], color='black', marker='X', s=200, label='Centroids')
58. plt.title('EM Algorithm Clustering with 3 Clusters')
59. plt.xlabel('Feature 1')
60. plt.ylabel('Feature 2')
61. plt.legend()
62. plt.grid(True)
63. plt.show()
64.  

 

결과는 위와 같다.