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 https://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;
  &
nbsp; 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로도 보기 힘든 코드가 나오는군.. 쩝