Gaussian Mixture Model & KMeans

각각 200개의 점들로 이루어진 8개의 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 (int k = 0; k < classes.size(); ++k)
            prob_cls[i][k] /= s;
    }
    //M-step;
    for (int 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 (int 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) ;
};