Geometry & Recognition

주어진 N개의 관측값을 k개의 가우시안 mixture로 모델링한 후 EM 알고리즘을 적용하여 모델을 결정하는 파라미터를 추정하고자 한다. 추정해야 할 모델 파라미터는 k개의 가우시안을 기술하는 평균값(μ)과 공분산 행렬(Covariance Matrix: Σ)과 이들 가우시안 분포의 mixing(α) 정도를 나타내는 파라미터이다.

다음은 이 과정을 수행하는 C++ class와 사용 예제, 및 결과를 보여준다.

//가우시안 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 * dxy) * 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};
    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|θ);