EM Algorithm: Line Fitting

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

 

 

두 직선을 각각 $y=a_1 x + b_1$, $y = a_2  x + b_2$라고 하면, $(a_1, b_1)$, $(a_2, b_2)$를 구하는 문제이다. 만약 각각의 data가 에러가 수반되는 측정에 의해서 얻어졌다고 하자. 에러 분포가 정규분포를 따른다면(여기서는 두 개의 모델 모두 갖은 표준편차를 갖는다고 가정했다) $i$-번째의 데이터가 각각의 직선 모델 1과 2에 속할 확률은 (posterior with equal prior) Bayes 공식에 의해서 

$$ w_1[i] = \frac{ e^{ - \frac{ r_1^2[i]}{2\sigma^2} } } {e^{ - \frac{ r_1^2[i]}{2\sigma^2} } + e^{ - \frac{ r_2^2[i]}{2\sigma^2} } }, \quad w_1[i] = \frac{ e^{ - \frac{ r_2^2[i]}{2\sigma^2} } } {e^{ - \frac{ r_1^2[i]}{2\sigma^2} } + e^{ - \frac{ r_2^2[i]}{2\sigma^2} } }, \quad i = 0,1,2,... $$

로 주어진다. 여기서 $r_1(i)$와 $r_2(i)$는 residual error이다:

$$r_1[i] = a_1 x[i] + b_1 - y[i],\quad r_2[i] = a_2 x[i] + b_2 - y[i], \quad i=0,1,2,...$$

(*이 값 대신에 직선까지 거리=$\frac{|a_k x + b_k - y|}{\sqrt{1+ a_k^2}},~ k=1,2)$로 대체해도 된다)

이제 각각의 데이터에 대해서 posterior를 구하였으므로(E-step) 이 값을 가중치로 하여서 직선의 방정식을 다시 갱신한다. 즉, 각각의 data 점들에 대한 $w_1(i)$를 가중치로 하여서 다시 직선 모델 1의 파라미터를 재계산하고, 동일한 과정을 $w_2(i)$를 가지고 직선 모델 2를 계산하여서 $(a_1, b_1)$, $(a_2, b_2)$를 재계산한다(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 = data.size(); i-- > 0;) {
        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) / 2;
        double n2 = SQR(r2) / SQR(sigma) / 2;
        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 = data.size(); i-- > 0;) {
        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;
}