'Expectation Maximization'에 해당되는 글 2건

  1. 2017.01.02 Expectation Maximization Algorithm for Two-Component Gaussian Mixture
  2. 2008.06.07 Gaussian Mixture Model (2)

주어진 히스토그램 데이터를 두 개의 가우시안 분포의 혼합으로 모델링하여 데이터를 분리하기 위한 decision boundary 값 (threshold value)을 expectation maximization 알고리즘을 적용하여 구한다. 

E-step: compute responsibility (of class 2; for class 1, 1-gamma_i)

M-step: compute the weighted means, variances and mixing probability

log-likelihood: 

void estimGaussParams(double data[], int start, int end, double *mean, double *var) {

void initGuess(double data[], int n, double mean[], double var[], double *mixprob) {

#define PI (4.0 * atan(1.))

double gaussDist(double x, double mean, double var) {   

double responsibility2(double x, double mean[], double var[], double mixprob) {   

double weightedMeanVar(double data[], double gamma[], int n, double mean[], double var[]) {

#define EPSILON  1e-6

// Expectation Maximization algorithm applied to Two component Gaussian Mixture Model;

double emTwoCompGMM(double data[], int n) {

    double mean[2], var[2], mixprob;

    double *gamma = new double [n];     // responsibilities for class 2;

    initGuess(data, n, mean, var, &mixprob);

    

    // begin algorithm;

    while (1) {

        // E-step;

        for (int i = 0; i < n; i++) 

            gamma[i] = responsibility2(i, mean, var, mixprob);


        double old_mixprob mixprob;


        // M-step;

        mixprob = weightedMeanVar(data, gamma, n, mean, var);

        TRACE("mixing probability= %f\n", mixprob);


        // check convergence(usually loglikelihood is tested);

        if (fabs(mixprob - old_mixprob) < EPSILON)

            break;

    }

    // estimate decision boundary;

    for (int k = n - 1; k >= 0; k--) 

        if (gamma[k] < 0.5) break;


    delete [] gamma;

    return (2 * k + 1) / 2.; // = average of {k, k+1};

}



Posted by helloktk

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