2차함수 추정 및 구현

주어진 데이터가 3개의 값을 가지고 있고 그 'x' 위치의 값이 독립적이면 Lagrange Interpolation 을 이용하여 값을 구할 수 있습니다.

그런데 주어진 값이 4개 이상일 경우 그리고 그 함수가 2차함수 형식이라고 알고 있는 경우 추정을 하려면 다른 방법을 생각해 보아야 합니다.

이 글에서는 그 과정을 한번 생각해 보고자 합니다.
다항식 추정을 위해서는 보다 보편적인 규칙을 찾아야 하지만, 이번에는 구현과정을 가장 기초적인 미분과 연립방정식을 이용하여 추정해 보고자 합니다.
잘 모르기도 하지만, 생각을 코드화 하는 과정을 논의해 보고 싶다는 생각이 들었기 때문입니다.

일단 알 수 있는 것은 2차 함수의 미분으로 1차 함수를 구성할 수 있고 그 함수는 미분한 'x' 접선의 기울기라는 사실입니다.
예를 들어 'y=2x^2 - 4x + 8' 를 미분하면 '4x-4' 가 됩니다.  x가 1일때 기울기가 0 이되니 이 지점이 최저점이 될 것입니다.

좌표는 위 수식에 따라
x = { -8, -6, -4, -2, 0, 2, 4, 6, 8 },
y = {168.0, 104.0, 56.0, 24.0, 8.0, 8.0, 24.0, 56.0, 104.0}
라고 가정해 보겠습니다.
물론 추세선이기 때문에 좌표가 함수에 일치하지는 않겠지만 일단 위와 같은 위치로 가정해 보겠습니다.

먼저 살펴보아야 할 순서를 생각 해 보겠습니다.

1. x 좌표를 오름차순으로 정렬하고 y좌표를 그에 맞춰 정렬하기(정렬되어 있다고 가정)
2. 앞 뒤 두 좌표의 기울기 및 x의 평균값으로 1차함수(미분된 값-기울기)를 설정
3. 2차함수의 계수 및 1차 함수의 계수를 평균값으로 구하기
4. 2차 계수와 1차 계수 값으로 상수 값(평균)구하기
를 생각해 볼 수 있습니다.

먼저 x좌표값에 의한 정렬은 되어 있다고 가정하겠습니다.
위의 예에서 보면 두점 사이의 기울기는 두점의 평균값위치에서 기울기와 같을 것이라는 가정을 해 보았습니다.
(104.0-168.0)/(-6-(-8)) = -32, (-6+(-8))/2 = -7
1차 미분식을 a*x+b 로 가정하고 a 가 1차 계수 b가 상수라고 하면
이는 x값이 -7일 때 기울기가 -32 란 의미니까 결국 2차식의 미분함수인 a*x+ b = -32 란 의미고 이는 -7*a + b = -32 란 의미가 됩니다.

(56.0-104.0)/(-4-(-6)) = -24, (-4+(-6))/2 = -5
이 부분을 위와 같은 방식으로 구성하면,
-5*a + b = -24 라는 식으로 구성할 수 있습니다. 두 방정식을 풀기 위해 b 값으로 하나를 환원하면 b = -24 + 5*a 라는 식으로 만들 수 있고 이를 위에 대입하면
-7*a + (-24+5*a ) = -32 입니다. 풀어 보면 a = 4를 도출 할 수 있습니다. b의 값은 a=4라는 값에 의해서 -5*4 + b = -24, b = -4의 값이 나오게 됩니다.

주어진 좌표에서 이와 같은 식으로 a, b를 계산해서 저장한 다음 a의 평균값, b의 평균값을 선택하게 됩니다.

다음은 c를 풀기 위해서 원래 함수인 2차 함수로 환원하여(미분된 함수 이기 때문에 다시 원래 함수의 값으로 변환하기 위해 a/2*x^2+ b*x + c 의 형태로 구성합니다.
결국  a = 4/2 , a = 2의 값이 됩니다.

주어진 좌표에서 예를 들면
x = -8일 때 y = 168 이니까, -8^2*2 - 4*-8 + c = 168입니다. 128 + 32 + c = 168, c = 8 입니다.
이렇게 좌표에서 c를 다 구하고 평균값을 취하면 2*x^2 - 4*x + 8 의 식을 도출 할 수 있습니다.

다음은 구현한 코드 입니다.

public static List<Double> getQuadracticEquations(List<? extends Number> xList
, List<? extends Number> yList ) {

if ( xList == null || yList == null || xList.size() < 3 || xList.size() != yList.size() ){
return null;
}

int size = xList.size();
for ( int i = size-1; i >= 0; i-- ) {
if ( xList.get(i) == null || yList.get(i) == null
|| Double.isNaN(xList.get(i).doubleValue()) || Double.isNaN(yList.get(i).doubleValue())
|| Double.isInfinite(xList.get(i).doubleValue()) || Double.isInfinite(yList.get(i).doubleValue()) ) {
xList.remove(i);
yList.remove(i);
}
}

size = xList.size();
if ( size < 3 ) {
return null;
}

List<Double> result = new ArrayList<Double>();

double a = Double.NaN;
double b = Double.NaN;
double c = Double.NaN;

double preXAvg = Double.NaN;
double preSlope = Double.NaN;
List<Double> aCoefList = new ArrayList<Double>();
List<Double> bCoefList = new ArrayList<Double>();
List<Double> cCoefList = new ArrayList<Double>();

for ( int i = 1; i < size; i++) {
double xGap = xList.get(i).doubleValue()-xList.get(i-1).doubleValue();
double yGap = yList.get(i).doubleValue()-yList.get(i-1).doubleValue();
double slope = (yGap/xGap);
double xa = (xList.get(i).doubleValue()+xList.get(i-1).doubleValue())/2.0;
if ( !Double.isNaN(preXAvg) ) {
double aV = ((slope)- preSlope)/(xa - preXAvg);
double bV = (yGap/xGap)-aV*xa;
// System.out.println ( "\tRESULT : " + aV + " : " + bV);
aCoefList.add(aV);
bCoefList.add(bV);
}
preXAvg = xa;
preSlope = slope;
}

a = StatisticsBaseUtil.getAverageValue(aCoefList)/2.0;
b = StatisticsBaseUtil.getAverageValue(bCoefList);

for ( int i = 0; i < size; i++) {
cCoefList.add(yList.get(i).doubleValue() - (Math.pow(xList.get(i).doubleValue(),2)*a + xList.get(i).doubleValue()*b));
}
c = StatisticsBaseUtil.getAverageValue(cCoefList);
result.add(a);
result.add(b);
result.add(c);

return result;
}

테스트 용도로 두가지를 구성해 보았습니다.

먼저 검증용 데이터 입니다.

double a = 2.0;

double b = -4.0;

double c = 8.0;

List<Double> xList = new ArrayList<Double>();

List<Double> yList = new ArrayList<Double>();

int size = 10;

double xx = -8.0;

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

double x = xx+i*2; 

xList.add(x);

yList.add(Math.pow(x,2)*a+b*x+c);

}

System.out.println ( xList + "\n" + yList );

System.out.println ( getQuadracticEquations(xList,yList) );

결과는 아래와 같습니다.

[-8.0, -6.0, -4.0, -2.0, 0.0, 2.0, 4.0, 6.0, 8.0, 10.0]

[168.0, 104.0, 56.0, 24.0, 8.0, 8.0, 24.0, 56.0, 104.0, 168.0]

[2.0, -4.0, 8.0] -> 계수들

다음은 데이터를 약간 흔들어 보았습니다.

double a = 2.0;

double b = -4.0;

double c = 8.0;

List<Double> xList = new ArrayList<Double>();

List<Double> yList = new ArrayList<Double>();

int size = 10;

double xx = -8.0;

Random rd = new Random();

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

double x = xx+i*2; 

xList.add(x);

yList.add(Math.pow(x,2)*a+b*x+c + rd.nextDouble()*10.0);

}

System.out.println ( xList + "\n" + yList );

System.out.println ( getQuadracticEquations(xList,yList) );

결과는 아래와 같습니다.

[-8.0, -6.0, -4.0, -2.0, 0.0, 2.0, 4.0, 6.0, 8.0, 10.0]

[169.37547437750422, 111.00550252621369, 59.85251104741083, 30.723567117837423, 11.018439045994608, 15.543492759715203, 26.299058613035104, 61.63218678109505, 104.34192659530801, 175.75788040181308]

[2.027905088403056, -7.334878784772051, 14.941109705660864]

유사한 추세를 보여 주긴 하지만 무엇인가 조금 어긋나 보이는 것 같습니다.

그래서 조금 각도를 틀어 생각해 보기로 하겠습니다.

각 좌표값의 변화가 평균적으로만 적용되기 때문에 전체적으로는 튀는 값이 보다 많이 적용될 수 있다는 생각을 해 볼 수 있습니다.

미분이 1차 방정식이기 때문에 위에서 구한 x의 위치 'xa'와 slope 을 이용해서 선형 최소제곱법을 적용해 보겠습니다.

public static List<Double> getQuadracticEquations(List<? extends Number> xList

, List<? extends Number> yList ) {

if ( xList == null || yList == null || xList.size() < 3 || xList.size() != yList.size() ){

return null;

}

int size = xList.size();

for ( int i = size-1; i >= 0; i-- ) {

if ( xList.get(i) == null || yList.get(i) == null 

|| Double.isNaN(xList.get(i).doubleValue()) || Double.isNaN(yList.get(i).doubleValue()) 

|| Double.isInfinite(xList.get(i).doubleValue()) || Double.isInfinite(yList.get(i).doubleValue()) ) {

xList.remove(i);

yList.remove(i);

}

}

size = xList.size();

if ( size < 3 ) {

return null;

}

List<Double> result = new ArrayList<Double>();

double a = Double.NaN;

double b = Double.NaN;

double c = Double.NaN;

double preXAvg = Double.NaN;

double preSlope = Double.NaN;

List<Double> cCoefList = new ArrayList<Double>();

List<Double> xxList = new ArrayList<Double>();

List<Double> yyList = new ArrayList<Double>();

for ( int i = 1; i < size; i++) {

double xGap = xList.get(i).doubleValue()-xList.get(i-1).doubleValue();

double yGap = yList.get(i).doubleValue()-yList.get(i-1).doubleValue();

double slope = (yGap/xGap);

double xa = (xList.get(i).doubleValue()+xList.get(i-1).doubleValue())/2.0;

xxList.add(xa);

yyList.add(slope);

}

List<Double> leastList = StatisticsBaseUtil.getLeastSquaresValues(xxList,yyList);

a = leastList.get(0)/2.0;

b = leastList.get(1);

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

cCoefList.add(yList.get(i).doubleValue() - (Math.pow(xList.get(i).doubleValue(),2)*a + xList.get(i).doubleValue()*b));

}

c = StatisticsBaseUtil.getAverageValue(cCoefList);

result.add(a);

result.add(b);

result.add(c);

return result;

}

결과는 다음과 같습니다.

[-8.0, -6.0, -4.0, -2.0, 0.0, 2.0, 4.0, 6.0, 8.0, 10.0]

[173.84227469091417, 106.62900546618599, 57.65556936199552, 32.05383516544042, 8.677565603651919, 8.091549432210913, 28.18148589027254, 62.75161780452815, 107.05732207887412, 173.9211479980657]

[2.0415744582349165, -4.078767066072524, 10.551372835299308]

또다른 방식으로 구한 계수는 다음과 같습니다.

[10.811528071220785, -4.012475708175465, 2.031973087828489]

그 결과를 R Graph 로 표현해 보겠습니다.

xGap 에 따른 오차등을 평균값으로 구성한 것에서 오류가 발생할 가능성이 있고, 그래서 그 오차를 줄이기 위해 계산한 x좌표와 slope을 x,y좌표로 재 구성하여 선형 최소제곱법을 수행해 보았습니다. 결과는 다른 방식으로 구한 계수값과 상당히 유사함을 보여 주고 있습니다.   x 사이값이 (평균값)이 해당 위치의 slope 이라는 가정은 위험하기는 하지만 직관적으로 보이기는 합니다.

3차, 4차 등으로 고차 함수로 추정해야 할 때 한계점도 있습니다.
다만, 이 글에서는 2차 함수로 추정하면서 코드화 하는 과정을 함께 살펴 보고 싶어서 구성한 글입니다.


다시 말씀 드리지만, 다른 추정방법과 그 결과값이 약간 다른 방법입니다.
코드의 StatisticsBaseUtil는 이전에 구성한 평균,분산,표준편차, 추세선에서 구현한 메소드 들입니다.