2차 다항식 회귀 (polynomial regression of order 2)

   회귀 (regression) 이란, 주어진 데이터 (x_k,y_k) 가 함수 y=f(x) 로부터 측정된 것이란 가정 하에, 주어진 데이터에 가장 잘 맞는 함수의 변수를 찾는 것을 의미한다. 다시 말해, 주어진 데이터 집합 

{ (x,y) }k = { (x1,y1), (x2, y2), (x3, y3), ..., (xN,yN) }

이 있고, 이 데이터가 y = ax + b  의 선형식으로부터 측정된 것이라 가정하였을 때, 데이터와 실제 함수와의 오차를 가장 작게 하는 a 와 b 는 무엇일까? 이 값을 찾아 내는 것이 regression 이고, 일반적으로 함수  f 를 다항식으로 가정하였을 경우 polynomial regression 이라 한다. 예를 들면 다음과 같다.

이와 같은 것을 수학적으로 나타내면 다음과 같다.

위에서 보는 바와 같이 regression model 을 y = a0 + a1x + ei 로 한다. 이 때 ei  는 실험 측정치에 포함되어 있는 오차이다. 이 모델에서 오차의 합을 가장 작게 하는 a0과 a1 을 찾아 내면 된다. solution 에 있는 바와 같이 주로 오차의 제곱의 합을 a0와 a1 의 함수로 가정하여 ai 로 편미분하여 0 이 되게 하는 a0, a1 을 찾아 낸다. 1차 다항식의 계산은 간단하니 여기서는 하지 않고(이 글에 있는 코드 중 get_correlation 함수가 1차 다항식의 회귀를 구현한 것이므로 참고할 것), 2차 다항식을 다루어 본다. 2차 다항식도 결국은 1차 다항식과 똑같은데, regression model 만 2차로 확장한다. 그 후 오차의 제곱의 합을 ai 로 편미분하여 0 이 되게 놓고 선형방정식을 풀면 된다.

위의 식을 maple 을 이용하여 표현하면 다음과 같다.

위의 선형방정식을 다음과 같은 형태로 변환하여,

풀면 다음과 같은 해를 얻는다.

지금까지의 과정을 C++ 코드로 작성하면 다음과 같다.

double square(double init, double x)

{

    return init + x*x;

}

double cubic(double init, double x)

{

    return init + x*x*x;

}

double forth_power(double init, double x)

{

    return init+x*x*x*x;

}

 

// Yi = b0 + b1*Xi + b2*Xi^2 + ei 으로, sum(e_i)를 최소화시키게끔

// regression 된 변수를 찾도록 정규방정식을 편미분해서 0 이 되는

// 각각의 bi를 구하면 다음과 같다.

bool get2ndOrderRegression(std::vector<double>* srcX, std::vector<double>* srcY, double *b0, double *b1, double* b2)

{

    double Y = std::accumulate(srcY->begin(), srcY->end(), 0.0);

    double X = std::accumulate(srcX->begin(), srcX->end(), 0.0);

    double X2 = std::accumulate(srcX->begin(), srcX->end(), 0.0, square);

    double X3 = std::accumulate(srcX->begin(), srcX->end(), 0.0, cubic);

    double X4 = std::accumulate(srcX->begin(), srcX->end(), 0.0, forth_power);

    double K = 0.0;

    double L = 0.0;

    int i = 0;

    int n = (int)srcX->size();

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

        K += ((*srcY)[i]*(*srcX)[i]*(*srcX)[i]);

        L += ((*srcY)[i]*(*srcX)[i]);

    }

 

    double denominator = -n*X4*X2 + X4*X*X + X2*X2*X2 + X3*X3*n - 2*X3*X*X2;

    double b0p = -(Y*X4*X2 - Y*X3*X3 - X*L*X4 + X*X3*K - X2*X2*K + X2*X3*L);

    double b1p = X*Y*X4 - X*K*X2 - L*n*X4 + X3*n*K - Y*X2*X3 + X2*X2*L;

    double b2p = -(K*n*X2 - K*X*X - X2*X2*Y - X3*n*L + X3*X*Y + X*X2*L);

    *b0 = b0p/denominator;

    *b1 = b1p/denominator;

    *b2 = b2p/denominator;

 

    return true;

}