'Line Fitting'에 해당되는 글 2건

  1. 2008.07.08 Robust Line Fitting
  2. 2008.06.29 EM Algorithm : Line Fitting 예

이미지에서 관찰된 점들이 {(xi, yi) | i = 1, 2,..., N}이 있다. 이 점들을 직선 y = a + b*x 로 피팅을 하고 싶을 때, 보통 최소 자승법을 이용하는 데,  원리는 직선의 방정식이 예측한 y값과 실제 관찰한 y값의 차이의 제곱(=square deviation)을 최소화시키는 직선의 파라미터 a, b를 찾는 것과 동일하다:


확률적으로는, 데이터를 얻는 과정인 측정에서 측정값 yi 는 랜덤에러를 가지게 되고, 이 에러는 참값 y(x) 근방에서 정규분포의 모양을 가지게 된다고 가정을 한다. 만약에 모든 측정의 에러가 동일한 표준편차 σ를 갖게 된다면, N개의 관측데이터가 나타날 확률(likelihood)는 (개별측정은 IID조건을 만족한다고 전제)
     

의 형태가 된다. 따라서 최소자승방법은 이 likelihood를 최대화시키는 파라미터를 찾는 방법이 된다. 이 최소 자승법은 피팅 파라미터를 주어진 데이터를 이용해서 표현할 수 있는 장점이 있지만, outliers에 대해서 매우 취약하다 (아래의 결과 그림을 참조). 이것은 적은 수의 outlier도 χ2 에 큰 기여를 할 수 있기 때문이다. 따라서 피팅을 좀더 robust하게 만들기 위해서는 이들 outlier들이 likelihood에 기여하는 정도를 줄여야 한다. 그러기 위해서는 likelihood의 지수인자가 큰 에러에서 덜 민감하게 반응하는 꼴로 바뀌어야 한다. 이런 형태 중에서 가장 간단한 것 중의 하나가 square deviation 대신에 absolute deviation을 이용하는 것이다:


그러나 이 식을 이용하는 경우에는 최소자승법과는 달리, 파라미터 a, b를 주어진 데이터 {(xi, yi)}로 표현할 수 없고, 반복적인 방법을 이용해서 찾는 수 밖에 없다. 


주어진 수열 {c
i}와 a 와의 차이의 절대값의 합
i |ci - a| 는 a가 이 수열의 median 값인 경우에 최소된다는 사실을 이용하면 (증명: 극값을 구하기 위해서 a에 대해서 미분하면, 0 = (i 1) - (i 1) : 합은 a가 ci 보다 큰 경우와 작은 경우로 분리=> 따라서 0이 되기 위해서는 작은 경우와 큰 경우의 수가 같아야 한다 => 즉, a = median{ci}), 고정된 b값에 대해서 위 absolute deviation을 최소화 시키는 a값은

 


임을 알 수 있다. 그리고 absolute deviation 식을 b에 대해서 미분을 하여서


을 얻는데 위의 a 값을 대입해서 위의 식을 만족시키는 b값을 bracketing and bisection 방법을 이용해서 얻을 수 있다(불연속함수여서, 일반적으로 근을 구하는 방법을 사용하는 것은 위험하다). 아래의 코드는 이를 구현한 것이다.


참고 : Numerical Recipes. 

// qsort-comparator; 

static int comp(const void *a, const void * b) { 

    double d = *((double *)a) - *((double *)b);

    return d > 0 ? 1 : d < 0 ? -1 : 0; 

// 기울기(bb)가 주어진 경우에 y-절편(median = aa)값을 추정하고, 이 때 AD값을 얻는다. 

double RhoFunc(double *x, double* y, int n, 

               double bb, double *aa, double *abdev) { 

    double *h = new double [n];

    for (int i = 0; i < n; i++)  h[i] = y[i] - bb * x[i];

    qsort(h, n, sizeof(double), comp); 

    // median; 

    *aa = (n & 1) ? h[n / 2] : (h[n / 2] + h[n / 2 - 1]) / 2;

    

    double sum = 0;  

    *abdev = 0; 

    for (i = 0; i < n; i++) { 

        double d = y[i] - (bb * x[i] + (*aa)); 

        *abdev += fabs(d); 

        // sign-함수의 원점에서 모호함을 없애기 위해서 증폭을 시킴; 

        if (y[i] != 0.) d /= fabs(y[i]);

        if (fabs(d) > DBL_EPSILON) // sign 함수의 모호함에서 벗어나면,

            sum += (d >= 0 ? x[i] : -x[i]); 

    } 

    delete [] h;

    return sum; // return sum{xi * sign(yi - (b * xi + a))} 

}; 

// 최소자승법을 이용한 라인추정: 

// return (sigma[dy] / sigma[x]);

double FitLine_LS(double *x, double *y, int n, double *a, double *b) { 

    double sx = 0, sy = 0, sxx = 0, sxy = 0; 

    for (int i = 0; i < n; i++) { 

        sx += x[i]; 

        sy += y[i]; 

        sxx += x[i] * x[i] ; 

        sxy += x[i] * y[i] ; 

    }; 

    double det = n * sxx - sx * sx; 

    if (det == 0.) return -1;                   // vertical line;

    *a = (sxx * sy - sx * sxy) / det; 

    *b = (n * sxy - sx * sy) / det;  

    double chi2 = 0;  

    for (i = 0; i < n; i++) { 

        double t = y[i] - (*a + *b * x[i]); 

        chi2 += t * t; 

    }  

    det /= n; // det -> var(x) * n; 

    // chi2 = var(dy) * n;

    // (dy vs x의 편차비)

    return  sqrt(chi2 / det);  ; 

}

// y = a + b * x ; 

// Least Absolute Deviation: 

void FitLine_MAD (double *x, double *y, int n, 

                  double *a, double *b, double *mad) 

    // least square esimates for (aa, bb);      

    double aa, bb; 

    double sigb = FitLine_LS(x, y, n, &aa, &bb); 

    double b1 = bb; 

    double abdev; 

    double f1 = RhoFunc(x, y, n, b1, &aa, &abdev); 

    /* bracket 3-sigma away in the downhill direction; */ 

    double b2 = fabs(3 * sigb); 

    b2 = bb + (f1 < 0 ? -b2 : b2); 

    double f2 = RhoFunc(x, y, n, b2, &aa, &abdev); 

    

    // if conditional added to take care of the case of a

    // line input into this function. It is needed to avoid an

    // infinite loop when (b1 == b2) (within floating point error)

    if (fabs(b2 - b1) > (sigb + 0.005)) { 

        // bracketing;

        while ((f1 * f2) > 0) { 

            bb = 2 * b2 - b1; 

            b1 = b2;

            b2 = bb;

            f1 = f2;  

            f2 = RhoFunc(x, y, n, b2, &aa, &abdev) ; 

        } 

    } 

    // refine until the error is a negligible number of std; 

    sigb *= 0.01; 

    while (fabs(b2 - b1)> sigb) { 

        // bisection; 

        bb = (b1 + b2) / 2.; 

        if ((bb == b1) || (bb == b2)) break ; 

        double f = RhoFunc(x, y, n, bb, &aa, &abdev) ; 

        if ((f * f1) >= 0) { 

            f1 = f; 

            b1 = bb; 

        } 

        else { 

            f2 = f; 

            b2 = bb; 

        } 

    } 

    *a = aa; 

    *b = bb; 

    *mad = abdev / n; 

}

// 붉은 선--> 최소자승법에 의한 피팅.: outlier에 매우 취약한 구조.
// 파란 선--> least absolute deviation을 이용한 피팅: outlier에 매우 robust하다.

"사용자

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

RANSAC : Circle Fit  (0) 2008.07.21
KMeans Algorithm  (0) 2008.07.19
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
Posted by helloktk

댓글을 달아 주세요

아래의 그림을 보면 우리는 데이터를 연결하는 두 개의 직선을 생각할 수 있을 것이다.  그럼 두개의 직선을 어떻게 얻을 것인다? 물론, 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

댓글을 달아 주세요