Gaussian Mixture Model

Gaussian Mixture Model(GMM)은 주어진 데이터가 하나의 가우시안 분포가 아니라 여러 개의 가우시안 분포가 혼합되어 생성되었다고 가정하는 확률 모델이다. 각 가우시안은 하나의 데이터 군집(cluster)을 나타내며, 전체 확률밀도함수는 여러 개의 가우시안 분포 $\mathcal N(y\, ;\mu,\Sigma)$의 가중합으로 표현된다.

$$p(y)=\sum_{j=1}^{k}\alpha_j  \mathcal N(y\, ;\mu_j,\Sigma_j)$$

여기서 \(k\)는 가우시안의 개수, \(\alpha_j\)는 각 가우시안의 혼합비(mixture weight), \(\mu_j\)는 평균, \(\Sigma_j\)는 공분산 행렬을 나타낸다. 혼합비는 \(\sum_j\alpha_j=1\)을 만족한다.

GMM은 각 데이터가 하나의 군집에 완전히 속한다고 가정하는 K-means와 달리, 각 데이터가 여러 가우시안에 속할 확률을 계산하는 확률적(soft) 군집화 방법이다. 따라서 서로 겹치는 군집도 자연스럽게 표현할 수 있으며, 군집의 크기나 방향, 분산이 서로 다른 경우에도 효과적으로 모델링할 수 있다.

GMM의 모수(평균, 공분산, 혼합비)는 일반적으로 EM 알고리즘을 이용하여 추정한다. E-step에서는 각 데이터가 각 가우시안에 속할 확률을 계산하고, M-step에서는 이 확률을 이용하여 평균, 공분산, 혼합비를 다시 추정한다. 이 두 과정을 반복하면 likelihood가 증가하는 방향으로 모수가 갱신되어 최종 모델이 얻어진다.

영상 처리에서는 컬러 분할(color segmentation), 배경 제거(background subtraction), 객체 추적, 영상 압축 등 다양한 분야에서 GMM이 널리 사용된다. 특히 배경 제거에서는 각 픽셀의 시간에 따른 밝기 변화를 여러 개의 가우시안으로 모델링하여, 반복적으로 나타나는 배경은 높은 확률의 가우시안으로 표현하고, 새롭게 등장한 물체는 전경으로 검출하는 데 많이 활용된다.

 

다음은 이 과정을 수행하는 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(std::vector<POINT> &pts, int sx, int sy, int w) {
        posterior.resize(pts.size());
        wgauss.resize(pts.size());
        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(std::vector<POINT> &pts, double prior = 0) {
        double sx = 0, sy = 0, sxx = 0, syy = 0, sxy = 0;
        for (int j = pts.size(); j-->0;) {
            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(std::vector<double> &px) {
        weight = 0;
        for (int j = px.size(); j-->0;) {
            posterior[j] = wgauss[j] / px[j];    
            weight += posterior[j];
        }
        weight /= px.size(); 
    } 
    //P(x|thetal) * prior;
    void setProb(std::vector<POINT> &pts) {
        for (int i = pts.size(); i-->0;)
            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(std::vector<POINT> &pts) {
    GaussKernel2D kernel[NKERNELS];
    int nclass = NKERNELS;
    double weights[20] = {1};
    std::vector<double> px(pts.size());
    double wsum = 0 ;

    for (int i = 0; i < nclass; i++) {
        kernel[i].init(pts, 400, 400, weights[i]);
        wsum += weights[i];
    };    
#define MAX_ITER 50
    for (int iter = 0; iter < MAX_ITER; ++iter) {   
        for (i = px.size(); i-->0;) px[i] = 0;
        for (int k = 0; k < nclass; k++){
            GaussKernel2D &gker = kernel[k];
            gker.setProb(pts); 
            for (int i = px.size(); i-->0;)
                px[i] += gker.wgauss[i] ;
        }        
        for (k = 0; k < nclass; k++) {
            kernel[k].estimate(px);
            kernel[k].getParams(pts);
        }
        //또는 log-likelihood를 계산하여서 그 변화가 적으면 loop-끝내면 된다..
    }
}

//참고 : 아래의 데이터는 사전에 라벨링이 된 것이 아니다. 컬러링은 한번 계산한 후에 분포에 맞게 컬러를 조절하여서 다시 계산한 것이다.

 

f(y|θ);