각각 200개의 점들로 이루어진 8개의 2차원 가우시안 군집을 랜덤하게 만들고, 이 2차원 데이터를 kmeans알고리즘을 써서 8개로 분할하였다. 아래의 시뮬레이션은 이 정보를 초기 조건으로 하여서 Gaussian Mixture Model(GMM)에 적용한 결과이다. 두 개의 군집에 대해서 kmeans결과와 GMM의 결과가 서로 많이 차이가 남을 보여준다.


코드 추가: 2010.02.23
struct Data2d {
    double x, y ;
    int id ;
    Data2d(){};
    Data2d(double x, double y) : x(x), y(y), id(-1){}
};
struct Gauss2d {   
    double cov[4];
    double mx, my ; //mean ;
    double mix;
    //
    double nfactor; //1/(2*Pi*sqrt(det));
    double det ;    //det(cov)
    double icov[4];
    void   prepare() ;
    //
    double pdf(double x, double y);
} ;
void Gauss2d::prepare() {
    //det(cov);
    det = cov[0]*cov[3] - cov[1]*cov[2];
    if(det < 1e-10) {
        AfxMessageBox("not converging");
        return ;
    };
    nfactor = 1./(2.*MPI*sqrt(det)) ;
    //inv(cov);
    icov[0] = cov[3] / det ;
    icov[1] = -cov[1] / det;
    icov[2] = -cov[2] / det;
    icov[3] = cov[0] / det;
}
double Gauss2d::pdf(double x, double y) {
    x -= mx ;
    y -= my ;
    double a = x*(icov[0]*x + icov[1]*y) +
               y*(icov[2]*x + icov[3]*y);
    return nfactor * exp(-0.5*a) ;
};
void init_classes(std::vector<Data2d>& data, std::vector<Gauss2d>& classes) {
    /*
    for(int i = 0; i < classes.size(); i++) {
        Gauss2d& cls = classes[i] ;
        cls.cov[0] = 10 + 50 * rand() / double(RAND_MAX);
        cls.cov[1] = 0;
        cls.cov[2] = 0;
        cls.cov[3] = 10 + 50 * rand() / double(RAND_MAX);
        cls.mx = 100 + 300 * rand() / double(RAND_MAX);  
        cls.my = 100 + 300 * rand() / double(RAND_MAX);  
        cls.mix = 1;
    }
    */
    KMeans(data, classes);
    //use kmeans to locate initial positions;
}
void test_step(std::vector<Data2d>& data,
               std::vector<Gauss2d>& classes,
               std::vector<std::vector<double> >& prob_cls)
{
    //E-step ;
    for(int k = 0; k < classes.size(); k++) {
        Gauss2d& cls = classes[k];
        cls.prepare();
        //
        for(int i = 0; i < data.size(); i++) {
            prob_cls[i][k] = cls.mix * cls.pdf(data[i].x, data[i].y);
        };
    }
    //normalize-->임의의 데이터는 각 어떤 클레스에 속할 활률의 합=1;
    for(int i = 0; i < data.size(); i++) {
        double s = 0;
        int bc = 0; double bp = 0;  // to determine membership(debug);
        for(int k = 0; k < classes.size(); ++k) {
            s += prob_cls[i][k];
            // find maximum posterior for each data;
            if(bp < prob_cls[i][k]) {
                bp = prob_cls[i][k] ;
                bc = k ;
            };
        }
        data[i].id = bc;
        // normalize to 1;
        for(k = 0; k < classes.size(); ++k)
            prob_cls[i][k] /= s;
    }
    //M-step;
    for(k = 0; k < classes.size(); k++) {
        Gauss2d & cls = classes[k];
        //get mean;
        double meanx = 0;
        double meany = 0;
        double marginal = 0;
        for(int i = 0; i < data.size(); i++) {
            meanx += prob_cls[i][k] * data[i].x ;
            meany += prob_cls[i][k] * data[i].y ;
            marginal += prob_cls[i][k];
        };
        cls.mx = meanx = meanx / marginal ;
        cls.my = meany = meany / marginal ;
        //get mixing;
        cls.mix = marginal / classes.size();
        //get stdev;
        double sxx = 0, syy = 0, sxy = 0;
        for(i = 0; i < data.size(); i++) {
            double dx = data[i].x - meanx ;
            double dy = data[i].y - meany ;
            sxx += prob_cls[i][k] * dx * dx ;
            syy += prob_cls[i][k] * dy * dy ;
            sxy += prob_cls[i][k] * dx * dy ;
        };
        //set covariance;
        cls.cov[0] = sxx / marginal;
        cls.cov[1] = sxy / marginal;
        cls.cov[3] = syy / marginal;
        cls.cov[2] = cls.cov[1]; //symmetric;
    }  
}
void test() {
    int max_iter = 100;
    int nclass = 8;
    int ndata = 500;
    std::vector<Gauss2d> classes(nclass);
    std::vector<Data2d> data(ndata);
    // prepare posterior space;
    std::vector<std::vector<double> > prob_cls;
    for(int i=0; i < data.size(); ++i) {
        prob_cls.push_back(std::vector<double>(classes.size()));
    } ;
   // generate data...
   ..................................
    //init_classes
    init_classes(data, classes) ;

    int iter = 0;
    do {
        iter++;
        test_step(data, classes, prob_cls);
    } while(iter < max_iter) ;
};

'Image Recognition' 카테고리의 다른 글

Watershed Algorithm Code  (0) 2010.03.19
Retinex 알고리즘  (11) 2010.02.03
Gaussian Mixture Model & KMeans  (4) 2010.01.30
Image Morphing  (0) 2010.01.24
Fant's Algorithm  (0) 2010.01.22
Affine Transformation  (0) 2010.01.20
Posted by helloktk

댓글을 달아 주세요

  1. 마틴 2010.06.09 15:33  댓글주소  수정/삭제  댓글쓰기

    좋은 정보 잘 보았습니다.
    그림으로 보니 이해가 빠르네요 ^^

  2. 다솜 2016.09.03 16:34  댓글주소  수정/삭제  댓글쓰기

    이상합니다.
    소스를 모두 작성해보았는데
    물론 에러 없이 ㅠㅠ
    근데 그래픽의 움직임이 없네요
    혹 위와 같이 움직이는 소스를 받아 볼수 있을까요?

    dasom_01@hanmail.net
    감사합니다. 좋은 하루 되세요

    • helloktk 2016.09.10 12:11 신고  댓글주소  수정/삭제

      블로그의 코드는 Gaussian Mixture 알고리즘만 구현된 것입니다. 시간에 따라 데이터가 클러스터링되는 과정을 보이려면 코드의 중간에 일정한 시간마다 결과를 화면에 갱신하게 하게 만들어 주는 timer routine(또는 event handler)을 넣어야 합니다. 아니면 일정한 iteration 횟수간격으로 중간결과를 화일로 저장한 한 후 GIF를 만들어주는 프로그램이나 GIF 클래스를 이용해서 움직이는 결과를 만들수도 있습니다. 이과정은 스스로 만들어 보는 것이 더 좋을 것입니다.

  3. 멀더 2016.10.12 06:37  댓글주소  수정/삭제  댓글쓰기

    좋은 정보 감사합니다. 위와 같은 그림으로 출력할 수 있는 예제를 전달 받을 수 있을지요?
    설명을 주셨지만 어떻게 처리해야 할지 감이 잘 오지 않아서요.. ㅠㅠ

    hicys76@gmail.com 감사합니다.