아래의 그림을 보면 우리는 데이터를 연결하는 두 개의 직선을 생각할 수 있을 것이다.  그럼 두개의 직선을 어떻게 얻을 것인다? 물론, ICA 를 이용하는것이 한가지 방법이 될 것이다. 여기서는 EM 알고리즘을 이용하여서 두개의 직선을 기술하는 기울기와 y-절편의 값을 구하는 방법을 알아본다.

사용자 삽입 이미지


두 개의 직선을 각각 y=a1 * x + b1 , y = a2 * x + b2,라고 하면, 문제는 (a1, b1), (a2, b2)를 찾는 구하는 문제이다. 만약에 각각의 data 점의 측정에 의한 에러분포가 정규분포를 따른다면(여기서는 두 개의 모델 모두 갖은 표준편차를 갖는다고 가정했다) i-번째의 데이터 점이 각각의 직선모델1 과 2에 속할 확률은 (posterior with equal prior) 베이즈 공식에 의해서 

로 주어진다. 여기서 r1(i)와 r2(i)는 residual error이다:

(*이 값 대신에 직선에서의 거리=로 대체하여도 된다)

이제 각각의 데이터에 대해서 posterior를 구하였으므로(E-step) 이 값을 가중치로 하여서 직선의 방정식을 다시 갱신한다. 즉, 각각의 data 점들에 대한 w1(i)를 가중치로 하여서 다시 직선모델1의 파라미터를 재계산하고, 동일한 과정을 w2(i)를 가지고 직선모델2를 계산하여서 (a1,b1), (a2, b2)를 재계산한다(M-step). 이 갱신된 파라미터를 이용하여서 다시 가중치를 구하는 작업과, 직선의 파라미터를 구하는 작업을 일정하게 수렴할 때까지 반복을 하는 과정을 취하면 된다.

아래의 그림은 3번만에 원하는 결과를 얻는 것을 보여준다. 직선의 파라미터에 대한 초기값과 residual error의 표준편차 파리미터에 대한 적절한 값의 선택은 중요하다.

사용자 삽입 이미지

//코드의 일부...
std::vector<CPoint> data ;                                         // data,
std::vector<double> w1(data.size()), w2(data.size());  // weights.
double a1, b1, a2, b2 ;                                              // line params;
double sigma ;
// E-step;
void calcWeights() {
    for(int i=0; i<data.size(); i++){
        double  x = data[i].x,
                y = data[i].y ;
        double r1 = a1 * x + b1 - y ;
        double r2 = a2 * x + b2 - y ;
        double n1 = SQR(r1) / SQR(sigma);
        double n2 = SQR(r2) / SQR(sigma);
        double p1 = exp( - n1);
        double p2 = exp( - n2);
        w1[i] = p1 / (p1 + p2);
        w2[i] = p2 / (p1 + p2);
    }
};
//  M-step
void estimModels() {
    double s1xx=0, s1x=0, s1=0, s1y=0, s1xy=0;
    double s2xx=0, s2x=0, s2=0, s2y=0, s2xy=0;
    for(int i=0; i<data.size(); i++){
        double  x = data[i].x,
                y = data[i].y ;
            s1xx += w1[i] * SQR(x);
            s1xy += w1[i] * x * y;
            s1x  += w1[i] * x ;
            s1y  += w1[i] * y ;
            s1   += w1[i];
            //
            s2xx += w2[i] * SQR(x) ;
            s2xy += w2[i] * x * y ;
            s2x  += w2[i] * x ;
            s2y  += w2[i] * y ;
            s2   += w2[i];
    };
    double det1=s1xx*s1 - SQR(s1x);
    double det2=s2xx*s2 - SQR(s2x);
    a1 =(  s1 * s1xy  - s1x * s1y )/det1 ;
    b1 =(  s1xx * s1y - s1x * s1xy)/det1 ;
    a2 =(  s2 * s2xy  - s2x * s2y )/det2 ;
    b2 =(  s2xx * s2y - s2x * s2xy)/det2 ;
}

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

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
Gaussian Mixture Model  (2) 2008.06.07
Posted by helloktk