'Gaussian Mixture Model'에 해당되는 글 4건

  1. 2017.01.02 Expectation Maximization Algorithm for Two-Component Gaussian Mixture
  2. 2010.01.30 Gaussian Mixture Model & KMeans (4)
  3. 2008.07.01 EM : Binarization
  4. 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

댓글을 달아 주세요

각각 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 감사합니다.

이미지의 히스토그램을 이용하여서 전경과 배경을 분리하는 이진화 과정은 가우시안 mixture model을 이용하여서 EM 알고리즘을 적용하기에 좋은 예다. 전경에 해당하는 픽셀값의 분포와 배경에 해당하는 픽셀값의 분포는 히스토그램상에 섞여서 나타나는데. 이것을 두 가우시안의 혼합으로 생각할 때  em 알고리즘은 각각의 믹싱정도와 각 가우시안의 평균 및 표준편차를 추정한다. mixing parameter πa (a=1, 2,..., nclass)로 표시하는 경우에 특정한 픽셀값에 대한 posterior는


로 쓸수 있다. 여기서 f(x; θ)는 정규분포를 의미한다. 

posterior 정보를 이용하면 mixing parameter πa와 평균, 분산은 다음식으로 갱신이 된다. H[i]는 이미지의 히스토그램을 나타낸다.


   


log-likelihood:

여기서 a는 클래스 인덱스이고, i는 히스토그램 bin 인덱스이다.



// mixing 클래스를 기술하는 클래스;

struct mixclass {
    double prob ;               //mixing parameter;
    double mean ;               // mean
    double var ;                //variance;
} ;
// N(mean, var);

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

// posterior; Pr(Zi = c | xi, Theta);
// 주어진 관측값 x이 클래스 cid에 속할 posterior;
double classprob(double x, int nclass, mixclass*  mclass, int cid) {

// posterior (class_prob[i][c]) table 만들기;
void update_class_prob(int nbins, double * hist, int nclass, mixclass* mclass, double ** class_prob) {
// E-step;  pi[c] = mixture parameter for class c;
// posterior를 이용해서 특정클래스의 mixing 정도를 계산;==> next prior;
void update_prob(int nbins, double * hist, int nclass, mixclass* mclass, double ** class_prob) {
// mu[c]; 클래스의 평균;
void update_mean(int nbins, double * hist, int nclass, mixclass* mclass,  double ** class_prob) {
// var[c]; 클래스의 분산;
void update_var(int nbins, double * hist, int nclass, mixclass* mclass, double ** class_prob) {
// M-step; 
void update_parameters(int nbins, double * hist, int nclass, mixclass* mclass, double ** class_prob) {
// initialization;
void init_em(int nbins, double * hist, int nclass, mixclass* mclass) {
// calculate log-likelihood;
double mixLLK(int nclass, mixclass* mclass) { 
// check termination condition;
bool check_tol(double llk, double llk_p, double  eps) {
// 입력은 이미지의 히스토그램;
double em(int nbins/*=256*/, double hist[/*256*/],
    int nclass/*=2*/, mixclass mclass[/*=2*/], double eps/*=1.e-10*/) {
    double llk = 0, prev_llk = 0;
    // allocate memory buffers for the posterior information;
    double ** class_prob = (double**)malloc(sizeof(double*) * nbins);
    class_prob[0] = (double*)malloc(sizeof(double) * nbins * nclass) ;
    for (int i = 1; i < nbins; i++) class_prob[i] = class_prob[i - 1] + nclass;

    // initialization of algorithm;
    init_em(nbins, hist, nclass, mclass);
    //
    do {
        prev_llk = llk;
        // E-step ;
        update_class_prob(nbins, hist, nclass, mclass, class_prob);
        // M-step;
        update_parameters(nbins, hist, nclass, mclass, class_prob);
        llk = mixLLK(nclass, mclass);
        // TRACE("mean1=%f, mean2=%f\n", mclass[0].mean, mclass[1].mean);
        TRACE("log-likelihood=%e\n", llk);
    } while (!check_tol(llk, prev_llk, eps));
    // clean ;
    free(class_prob[0]);
    free(class_prob) ;
    return llk;
};
  • 적색 : 히스토그램 
  • 청색, 녹색 : posterior(membership); 
  • Otsu 알고리즘을 쓰는 경우에 100에서 threshold 값이 결정되고 EM은 110 정도임.

사용자 삽입 이미지

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

Robust Line Fitting  (0) 2008.07.08
Bayesian Spam Filtering  (0) 2008.07.03
EM : Binarization  (0) 2008.07.01
EM Algorithm : Line Fitting 예  (0) 2008.06.29
Shuffling  (0) 2008.06.21
Bayesian Decision Theory  (1) 2008.06.17
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

댓글을 달아 주세요

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

    좋은 자료 감사합니다!

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

    좋은 글 잘 읽었습니다 ^^

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

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