EM 클러스터링 알고리즘 코드

정확한 코드는 이곳에서 제공하고 있습니다.
—————————————————–
얼마만의 코드 관련 포스팅인지 모르겠다. 휴우~~~

어제 boost/math 라이브러리에 확률관련 모듈이 있는걸 보다가 문득 생각이 들어 EM 클러스터링 알고리즘을 구현해 봤다.

EM 알고리즘은 k-means알고리즘과 접근 방식이 100% 동일하나 단 확률과 통계 관련 지식이 조금 필요하다. 특히나 likelihood 관련 개념 이해는 구현을 위해 필수적이다.

이 알고리즘의 구동 과정은 아래와 같다.

1. 각 클러스터 매개 변수의 초기 값을 할당한다. 랜덤으로 하든 다른 방법으로 하든 마음이지만 역시 초기치가 나중 결과 정확도에 큰 영향을 미치긴 하더라.

2. repeat
3. expectation step : 각 데이터들에 대해서 클러스터에 속할 확률을 구한다. P(클러스터|data) 를 구하는데 이 부분은 Bayes rule을 이용하면 계산이 가능하다.
4. maximization step : 모든 데이터에 대해서 어느 클러스터에 포함되는지 확률값을 계산하고 나서 각 클러스터에 속한 데이터 점들을 이용해 최대 우도(likelihood)의 값을 가지는 매개 변수들을 재 계산한다. 요건 간단하게 3번 스탭에서 계산된 확률 값을 이용해 가중평균과 가중표준편차를 구하면 된다.
5. until 각 클러스터의 매개 변수들의 값들이 변하지 않는 시점까지

[CODE C++]
/*
EM clustering algorithm
Created by gogamza 2009.06.10
You can see full description in http://freesearch.pe.kr/1243
*/
#include <iostream>
#include <boost/random.hpp>
#include <vector>
#include <boost/math/distributions/normal.hpp>
#include <cmath>
#include <cassert>

using boost::math::normal;
using namespace std;

//data structure of data
typedef struct{
double datapoint;
int clusterid;
double prob;
} t_data;

//data structure of each cluster
typedef struct{
double mean;
double stddev;
int clusterid;
} t_cluster;

//cluster vector
vector<t_cluster> cluster;

//clacualte bayes probability
double getProb(const int x, const int clusterid){
if(cluster.size() == 0)
assert(0);
double priorP = 1.0/cluster.size();

vector<t_cluster>::const_iterator iter;
double numer = 0.0;
double denorm = 0.0;
for(iter = cluster.begin();iter != cluster.end();iter++){
if(clusterid == iter->clusterid)
numer = priorP * pdf(normal(iter->mean, iter->stddev), x);
denorm += priorP * pdf(normal(iter->mean, iter->stddev), x);
}
return numer/denorm;
}

int expectation(t_data& data){
vector<t_cluster>::iterator iter;
double x = data.datapoint;
double max = 0.0;
int clusterid = -1;
for(iter = cluster.begin();iter != cluster.end();iter++){
double prob = getProb(x, iter->clusterid);
if(max < prob){
max = prob;
clusterid = iter->clusterid;
}
}

data.clusterid = clusterid;
data.prob = max;
if(clusterid == -1)
return 0;
else
return 1;
}

void maximization(std::vector<t_data>& datas){
vector<t_cluster>::iterator iter;
for(iter = cluster.begin();iter != cluster.end();iter++){
vector<t_data>::iterator iter2;
double numer = 0.0;
double denorm = 0.0;
double mean = 0.0;
double stddev = 0.0;
double numer1 = 0.0;
double denorm1 = 0.0;
vector<t_data> temp;
for(iter2 = datas.begin();iter2 != datas.end();iter2++){
if(iter->clusterid == iter2->clusterid){
numer += iter2->datapoint * iter2->prob;
denorm += iter2->prob;
temp.push_back(*iter2);
}
mean = numer/denorm;

vector<t_data>::iterator iter3 = temp.begin();
while(iter3 != temp.end()){
numer1 += iter3->prob * pow(mean – iter3->datapoint, 2);
denorm1 += iter3->prob;
iter3++;
}

stddev = sqrt(numer1/denorm1);
iter->mean = mean;
iter->stddev = stddev;

}
}
}

int main( int argc, char **argv ){
//init
boost::mt19937 Random4;

vector<t_data> datas;

boost::normal_distribution<> norm_gen1(7,3);
boost::normal_distribution<> norm_gen2(3,3);
boost::normal_distribution<> norm_gen3(20,3);

boost::variate_generator< boost::mt19937, boost::normal_distribution<> > genG1( Random4, norm_gen1 );

boost::variate_generator< boost::mt19937, boost::normal_distribution<> > genG2(Random4, norm_gen2);

boost::variate_generator< boost::mt19937, boost::normal_distribution<> > genG3(Random4, norm_gen3);

t_cluster init_cluster1 = {1,3,0};
cluster.push_back(init_cluster1);
t_cluster init_cluster2 = {4,3,1};
cluster.push_back(init_cluster2);
t_cluster init_cluster3 = {18,3,2};
cluster.push_back(init_cluster3);
//random number generation
for( int i = 0; i < 1000; ++i ){
t_data data1 = {0, 0.0, 0.0};
data1.datapoint = genG1();
datas.push_back(data1);

t_data data2 = {0, 0.0, 0.0};
data2.datapoint = genG2();
datas.push_back(data2);

t_data data3 = {0, 0.0, 0.0};
data3.datapoint = genG3();
datas.push_back(data3);
}

//run
int runcnt = 0;
while(runcnt < 50){
vector<t_data>::iterator iter = datas.begin();
//expectation
for(;iter != datas.end(); iter++){
if(0 == expectation(*iter))
cerr << “error occur!\n” << endl;
}
//maximization
maximization(datas);
cout << “iter : ” << runcnt << endl;
vector<t_cluster>::iterator iter2 = cluster.begin();
for(; iter2 != cluster.end(); iter2++){
cout << “cluster id : ” << iter2->clusterid << ‘\n’
<< “mean : ” << iter2->mean << ‘\n’
<< “stddev :” << iter2->stddev << “\n\n” << endl;
}

runcnt++;
}

return 0;
}

[/CODE]

코드에서는 반복 종료 조건을 가장 간단하게 50회 반복으로 책정해 두었다.
main 함수의 init 부분 코드는 정규분포에 근거한 데이터 포인트 generation을 위한 코드가 대부분이다.

초기 generation 데이터들의 평균, 표준편차가

boost::normal_distribution<> norm_gen1(7,3);
boost::normal_distribution<> norm_gen2(3,3);
boost::normal_distribution<> norm_gen3(20,3);

위와 같다.

이미 정답이 나온 상태인데…

그럼 이를 모른다고 가정하고 프로그램을 돌려 나온 결과를 보면…

iter : 15
cluster id : 0
mean : 2.61313
stddev :2.33209

cluster id : 1
mean : 8.94453
stddev :1.84048

cluster id : 2
mean : 20.1655
stddev :2.9029

클러스터링 결과로
(2.6, 2.3), (8,9, 1.8), (20.2, 2.9) 가 나왔는데,

초기 생성 값인
(3,3), (7,3), (20,3) 과 상당히 유사함을 알 수 있다.

점심 시간을 이용해 간단하게 만들어 보려 했는데 핵심 코드보다 데이터 생성 및 구조체 할당 등등의 잡스런 과정이 거의 소스코드의 반 이상을 먹는것을 볼 수 있다.

역시나 데이터 처리를 많이 하는 코드는 코드가 상당히 볼품 없다는 것을 또 느끼는 찰라이다.
이걸 R로 구현하신 분이 있던데… 훔… R로도 보기 힘든 코드가 나오는군.. 쩝

from future import dream

당신의 나의 뜨거운 감자!

관련