'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 알고리즘  (10) 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

이미지의 히스토그램을 이용하여서 전경과 배경을 분리하는 이진화 과정은 가우시안 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