N개의 관측값이 주어지고, 이것들이 k개의 가우시안의 mixture로 생각을 할 때, EM알고리즘을 적용하여서 모델을 결정하는 파라미터를 추저하고자 한다. 모델 파라미터는 K개의 가우시안을 기술하는 평균값과 표준편차(Covariance Matrix)와 더블어 mixing정도를 나타내는 파라미터 K개를 추정하여야 한다.

사용자 삽입 이미지


아래는 이 일련의 과정을 수행하는 클래스와 사용예제, 및 결과샘플을 보여준다.

//가우시안 kernel Class;
class GaussKernel2D {
    double mx, my ;
    double sdx, sdy, sdxy;
    double weight;
public:
    std::vector<double> posterior;
    std::vector<double> wgauss;

    void init(POINT* pts, int np, int sx, int sy, int w) {
        posterior.resize(np);
        wgauss.resize(np);
        weight=w ;
        // initialize model parameters(random하게 선택함-->try another method!!::
        // 주어진 데이터로 전체 분포의 범위와 중심을 알수 있고, 이것을 주어진 클래스 수 만큼으로 임의로
        // 분배하는 방식으로 초기 조건을 설정하면 보다 안정적이다)
        mx= sx * double(rand())/RAND_MAX;
        my= sy * double(rand())/RAND_MAX;
        sdx= sx/4 + 1;
        sdy= sy/4 + 1;
        mx += rand() % 100 ;
        my += rand() % 100 ;
        sdxy=0;
    }
    double gauss2d(double x, double y) {
        double varx = sdx * sdx;
        double vary = sdy * sdy;
        double det = varx*vary - sdxy*sdxy;
        double idet = 1.0/det ;
        double dxx = (x - mx)*(x - mx);
        double dyy = (y - my)*(y - my);
        double dxy = (x - mx)*(y - my);
        return (1.0/sqrt(det)/6.28319)*exp(-0.5*(dxx*vary+dyy*varx-2.*dxy*sdxy)*idet);
    }
    void getParams(POINT *pts, int np, double prior=0)
    {
        double sx=0, sy=0, sxx=0, syy=0, sxy=0;
        for(int j = 0; j < np; j++) {
            double x = pts[j].x, y = pts[j].y;
            sx += posterior[j]* x;
            sy += posterior[j]* y;
            sxx += posterior[j]* x * x;
            syy += posterior[j]* y * y;
            sxy += posterior[j]* x * y;
        }
        double denorm = np * weight;
        mx= sx/denorm; my= sy/denorm;
        double devx=sxx/denorm - mx*mx ;
        if(devx<=0) devx=0.001;
        sdx=sqrt(devx);
        double devy=syy/denorm - my*my;
        if(devy<=0) devy=0.001;
        sdy=sqrt(devy);
        sdxy = sxy / denorm - mx*my;
        //if prior non-zero--> weight = weight * (1-alpha) + alpha * prior;   alpha=0.1?
    };
    //weight; //posterior;
    void estimate(double *px, int np) {
        weight = 0;
        for(int j = 0; j < np; j++) {
            posterior[j] = wgauss[j] / px[j];   
            weight += posterior[j];
        }
        weight /= np;
    } 

    //P(x|thetal) * prior;
    void setProb(POINT *pts, int np) {
        for(int i=0; i<np; i++){
            wgauss[i] = weight * gauss2d(pts[i].x, pts[i].y);
        }
    }
    void Draw(CDC* pDC, DWORD color=RGB(0xFF,0,0)) {
        CPen pen0(PS_SOLID,1,color);
        CPen* pOld=pDC->SelectObject(&pen0);
        drawCon(pDC, mx, my, sdx, sdy, sdxy); // draw ellipses;
        pDC->Ellipse(mx-2,my-2,mx+2,my+2);     // draw center of ellipse;
        pDC->SelectObject(pOld);           
    }
} ;
//
void em_main(POINT *pts, int np) {
    GaussKernel2D kernel[NKERNELS];
    int nclass=NKERNELS;
    double weights[20] = {1,1,1,1,1,1};
    std::vector<double> px(np);
    double wsum = 0 ;

    for(int i=0;i<nclass; i++){
        kernel[i].init(pts, np, 400, 400, weights[i]);
        wsum += weights[i];
    };    
#define MAX_ITER 50
    for(int iter = 0; iter < MAX_ITER; ++iter){  
        for(i=0; i<np; i++) px[i]=0;
        for(int k=0; k<nclass; k++){
            GaussKernel2D & gker = kernel[k];
            gker.setProb(pts, np);
            for(int i=0; i<np; i++){
                px[i] += gker.wgauss[i] ;
            }
        }        
        for(k=0; k<nclass; k++) {
            kernel[k].estimate(&px[0], np);
            kernel[k].getParams(pts, np);
        }
        //또는 log-likelihood를 계산하여서 그 변화가 적으면 loop-끝내면 된다..
    }
}
//참고 : 아래의 데이터는 사전에 라벨링이 된것이 아니다. 컬러링은 한번 계산한 후에 분포에 맞게 컬러를 조절하여서 다시 계산한 것이다.

사용자 삽입 이미지

 

f(y|θ);

사용자 삽입 이미지


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

Shuffling  (0) 2008.06.21
Bayesian Decision Theory  (1) 2008.06.17
Gaussian Mixture Model  (2) 2008.06.07
Rasterizing Voronoi Diagram  (0) 2008.05.26
RANSAC Algorithm  (0) 2008.05.24
Contour Tracing  (0) 2008.05.22
Posted by helloktk

댓글을 달아 주세요

  1. NoWhere 2009.04.05 01:14  댓글주소  수정/삭제  댓글쓰기

    좋은 자료 감사합니다!

  2. bbrto21 2011.11.18 01:19  댓글주소  수정/삭제  댓글쓰기

    좋은 글 잘 읽었습니다 ^^

    저는 GMM 을 응용한 프로젝트를 진행 중인 대학생입니다.

    혹시 GMM 관련하여 참고할 만한 자료나 예제 소스 같은 것을 공유 해 주실수 있으십니까?