비선형(Non-linear) 함수의 근사(Fitting) (Part 2)

(지난 회 게시물에서 바로 이어집니다)

그러면 오늘은 지난 회 게시물에서 제시했던 예제 데이터를 그대로 MPFITFUN 함수에 적용하여 근사 결과를 얻는 과정으로 바로 들어가겠습니다. 물론 이미 언급했듯이 MPFITFUN을 사용할 수 있으려면, MPFIT 라이브러리를 받아서 특정 디렉토리에 설치한 후 IDL 경로상에 그 디렉토리를 추가하는 작업이 반드시 선행되어야 합니다. 그 방법에 관해서는 지난 회 게시물의 서두에서 언급하였으므로 혹시라도 못보셨다면 그 내용을 참조하시기 바랍니다. MPFITFUN 함수의 소스코드인 mpfitfun.pro 파일의 내용에 의하면 기본적인 사용 문법은 다음과 같습니다.

result = MPFITFUN(MYFUNCT, X, Y, ERR, start_params, ...)

이와 같이 총 5개의 필수인자들이 반드시 필요합니다. 물론 지금 5개의 필수인자로 투입되어야 할 항목들은 이전 게시물의 내용에서 우리가 이미 만들었기 때문에 그것들을 그대로 사용하면 됩니다. 일단 이전에 LMFIT 함수가 사용되었던 내용을 다음과 같이 MPFITFUN 함수가 사용되는 내용으로 대체해야 합니다. 이 내용만 일단 보면 다음과 같습니다.

coefs = MPFITFUN('my_nl_funct', X, Y, errors, start, NITER=n_iter)

그래서 이렇게만 하면 다 끝날 것 같은데, 아직은 끝이 아닙니다. 중요한 것이 하나 있는데 바로 my_nl_funct, 즉 근사할 수식에 해당되는 사용자 정의 함수의 내용을 바꿔줘야 합니다. 물론 지난 회 게시물에서처럼 편미분을 구하는 등의 번거로운 작업은 다행히 전혀 필요가 없습니다. 그냥 다음과 같이 수식의 형태를 그대로 따라서 단순하게 정의해주면 됩니다.

FUNCTION my_nl_funct, X, P

Y = P[0]*SIN(P[1]*X)

RETURN, Y

END

즉 이와 같이 X와 P를 받아서 Y를 계산하여 돌려주는 형태의 함수로 정의해주면 됩니다. 이와 같이 MPFITFUN을 사용할 경우에는 편미분을 일일이 산출하여 명시하는 것이 전혀 필요가 없습니다. 이것만으로도 부담이 크게 줄어드는 셈입니다. 어쨌든 이러한 변동사항들로 인하여 주프로그램과 부프로그램의 내용이 전반적으로 이전 게시물의 그것과는 달라져야 합니다. 전체 내용은 다음과 같습니다.

FUNCTION my_nl_funct, X, P

Y = P[0]*SIN(P[1]*X)

RETURN, Y

END

PRO test_nonlinear_fitting

X = FINDGEN(101)/50

var1 = RANDOMU(-1, 101)*0.4-0.2

var2 = RANDOMU(-2, 101)*0.6-0.3

Y = (2.3+var1)*SIN((5.7+var2)*X)

errors = MAKE_ARRAY(101, VALUE=0.5)

start = [2.0, 5.0]

yrn = [-3, 3]

win = WINDOW(DIMENSIONS=[600, 500], /NO_TOOLBAR)

pl = PLOT(X, Y, YRANGE=yrn, SYMBOL='circle', /SYM_FILLED, LINESTYLE=6, /CURRENT)

coefs = MPFITFUN('my_nl_funct', X, Y, errors, start, NITER=n_iter)

chisq = TOTAL((Y-my_nl_funct(X, coefs))^2*ABS(1D/errors^2) )

HELP, coefs

PRINT, coefs

PRINT, start

PRINT, chisq

PRINT, n_iter

plo = PLOT(X, my_nl_funct(X, coefs), COLOR='red', THICK=2, $

  TITLE='MPFITFUN Result', FONT_SIZE=11, /OVERPLOT)

END

이 전체 내용을 보면 전반적인 틀은 그대로이지만 my_nl_funct 함수의 내용이 변경되었고, LMFIT 대신 MPFITFUN 함수가 사용되었습니다. 그리고 MPFITFUN 함수의 결과로 돌려받게 되는 coefs의 내용이 LMFIT 함수를 사용할 경우와는 다릅니다. LMFIT 함수가 돌려주는 coefs는 X의 각 포인트에 해당되는 근사 계산값들로 구성된 배열이었지만, MPFITFUN 함수가 돌려주는 coefs는 근사된 함수의 계수값들로 구성된 배열이 됩니다. 이것은 coefs의 값들을 출력해보면 알 수 있습니다. 출력된 coefs의 값들은 다음과 같습니다.

      2.25254      5.70780

이전 게시물의 내용을 떠올려보면 LMFIT을 사용할 경우 계수 초기 추정값 배열인 start의 값들이 나중에 변경되는 방식이었던 것과는 다릅니다. 어쨌든 이 계수값들을 보면 LMFIT으로 얻은 결과와 거의 유사하다는 것을 알 수 있습니다. 근사된 함수의 패턴을 그림으로 표출하려면 맨 하단의 PLOT 함수에서 Y에 해당되는 부분을 위와 같이 my_nl_funct(X, coefs)로 처리하면 됩니다. 이렇게 그려진 그림은 다음과 같습니다.

보시다시피 이전 게시물에서 봤던 LMFIT을 사용한 결과와 매우 유사합니다. 어차피 계수값들이 거의 유사하므로 함수의 패턴도 유사할 수 밖에 없습니다. 그리고 계산 과정에서 iteration이 실제로 반복된 횟수를 담은 변수 n_iter의 값은 5인데, LMFIT을 사용했을 경우에는 이 횟수값이 7이었습니다. LMFIT와 MPFITFUN 함수의 내부 알고리즘이 똑같지는 않기 때문에 이런 정도의 차이는 충분히 발생 가능합니다. Chi-Square 값의 경우는 LMFIT 함수의 경우 이 값을 계산하여 되돌려주는 전담 키워드(CHISQ)가 있었던 반면, MPFITFUN에는 제가 찾아본 바로는 그런 역할을 하는 키워드가 따로 없습니다. 그래서 지금은 Chi-Square 값을 직접 계산해야 하는데, 위의 프로그램에 있는 다음과 같은 내용이 그 역할을 합니다.

chisq = TOTAL((Y-my_nl_funct_2a(X, coefs))^2*ABS(1D/errors^2) )

이렇게 계산된 Chi-Square 값은 39.27로 LMFIT에서 얻은 값과 역시 유사합니다. 그래서 LMFIT과 MPFITFUN은 사용하는 방식은 다르지만 최종적인 계산 결과는 거의 유사한 것으로 보입니다. 물론 iteration 횟수나 소수점 하위 자리의 숫자상의 차이는 약간 있을 수 있습니다. 다만 이러한 두 함수 사이의 결과 차이가 근사할 함수의 형태나 데이터 포인트들의 분포 등에 따라 더 커질 수도 있는가에 대해서는 저도 아주 다양한 경우들을 테스트해본 것은 아니어서 정확히는 모르겠습니다. 어쨌든 한번 정리해보는 의미에서 LMFIT을 사용한 경우와 MPFITFUN을 사용한 경우의 산출값들을 다시 한번 비교해봅시다. 먼저 이전 게시물에서 봤던 LMFIT을 사용한 경우입니다.

계수들 :       2.2525569       5.7077310

Chi-Square :       39.274361

Iteration 횟수 :           7

그리고 MPFITFUN을 사용한 경우입니다.

계수들 :      2.25254      5.70780

Chi-Square :       39.274379

Iteration 횟수 :           5

자 그러면 이번에는 전혀 다른 예제 데이터들을 역시 전혀 다른 형태의 비선형 함수에 근사하는 경우를 보겠습니다. 다만 원본 데이터 및 함수 패턴이 이전과 달라진 것을 제외하면 전반적인 과정 자체는 유사합니다. 따라서 이 내용은 전체 코드를 바로 제시하는 정도로만 소개하겠습니다. 역시 LMFIT를 사용한 경우와 MPFITFUN을 사용한 경우로 나눠서 소개합니다. 먼저 LMFIT을 사용한 경우부터 봅시다.

FUNCTION my_nl_funct, X, P

result = [P[0]*EXP(P[1]*X)+P[2]+P[3]*SIN(X), EXP(P[1]*X), $

  P[0]*EXP(P[1]*X)*X, 1.0, SIN(X)]

RETURN, result

END

PRO test_nonlinear_fitting

X = FINDGEN(40)/20.0

Y = 8.7*EXP(-12.9*X)+2.6+5.3*SIN(X)+RANDOMU(-1, 40)-0.5

errors = MAKE_ARRAY(40, VALUE=1.0)

start = [10.0, -5.1, 6.0, 4.0]

yrn = [0, 15]

win = WINDOW(DIMENSIONS=[600, 500], /NO_TOOLBAR)

pl = PLOT(X, Y, YRANGE=yrn, SYMBOL='circle', /SYM_FILLED, LINESTYLE=6, /CURRENT)

coefs = LMFIT(X, Y, start, MEASURE_ERRORS=errors, FUNCTION_NAME='my_nl_funct', $

  ITER=n_iter, CHISQ=chisq, /DOUBLE)

HELP, coefs

PRINT, coefs

PRINT, start

PRINT, chisq

PRINT, n_iter

plo = PLOT(X, coefs, COLOR='red', THICK=2, TITLE='LMFIT Result', $

  FONT_SIZE=11, /OVERPLOT)

END

이번 예제에서는 SIN과 Exponential이 공존하는 함수의 형태를 대략 따라가는 가상의 예제 데이터를 먼저 생성하고, 이 데이터를 동일한 형태의 함수로 근사하고 있습니다. 주프로그램의 내용에서 알 수 있듯이 X는 0.0~2.0의 범위내에 있고 Y는 기본적으로 다음과 같은 형태를 갖지만 역시 RANDOMU 함수를 함께 사용하여 어느 정도의 퍼짐(Spread)가 있도록 인위적으로 손질을 좀 했습니다.

Y = 8.7*EXP(-12.9*X)+2.6+5.3*SIN(X) + (random spread)

그래서 이 식에서 8.7, -12.9, 2.6, 5.3의 자리에 미지의 계수 4개가 있다고 가정하고, 이 계수들을 찾는 과정이 됩니다. 바로 최종 결과 그림을 보면 다음과 같습니다.

그리고 출력된 계수들, Chi-Square값, iteration 횟수는 다음과 같습니다.

계수들 :       8.1715688      -11.944544       2.6251145       5.3703969

Chi-Square :       2.8507093

Iteration 횟수 :           7

이번에는 동일한 예제를 MPFITFUN을 사용하여 근사하는 경우입니다.

FUNCTION my_nl_funct, X, P

Y = P[0]*EXP(P[1]*X)+P[2]+P[3]*SIN(X)

RETURN, Y

END

PRO test_nonlinear_fitting

X = FINDGEN(40)/20.0

Y = 8.7*EXP(-12.9*X)+2.6+5.3*SIN(X)+RANDOMU(-1, 40)-0.5

errors = MAKE_ARRAY(40, VALUE=1.0)

start = [10.0, -5.1, 6.0, 4.0]

yrn = [0, 15]

win = WINDOW(DIMENSIONS=[600, 500], /NO_TOOLBAR)

pl = PLOT(X, Y, YRANGE=yrn, SYMBOL='circle', /SYM_FILLED, LINESTYLE=6, /CURRENT)

coefs = MPFITFUN('my_nl_funct', X, Y, errors, start, NITER=n_iter)

chisq = TOTAL((Y-my_nl_funct(X, coefs))^2*ABS(1D/errors^2) )

HELP, coefs

PRINT, coefs

PRINT, start

PRINT, chisq

PRINT, n_iter

plo = PLOT(X, my_nl_funct(X, coefs), COLOR='red', THICK=2, $

  TITLE='MPFITFUN Result', FONT_SIZE=11, /OVERPLOT)

END

최종 결과 그림은 다음과 같습니다.

그리고 출력된 계수들, Chi-Square값, iteration 횟수는 다음과 같습니다.

계수들 :      8.17166     -11.9435      2.62498      5.37054

Chi-Square :       2.8507078

Iteration 횟수 :           7

이번에도 역시 LMFIT의 결과와 MPFIT의 결과가 서로 거의 유사한 것을 볼 수 있습니다. 그리고 Chi-Square 값을 사용하여 데이터와 함수 사이의 유사성을 다음과 같은 방식으로 계산할 수 있습니다. 여기서는 DoF가 데이터 포인트 갯수 40에서 계수의 갯수 4를 뺀 값이된다는 점만 유의하면 됩니다.

PRINT, 1-CHISQR_PDF(2.85, 40-4)

그런데 이 계산값은 1.0이 나옵니다. 유사성 확률이 100%란 얘기입니다. 역시 너무 과하게 "좋은" 결과인 것 같습니다. 이렇게 좋게 나온 이유는 errors에서 정의된 에러분포의 표준편차가 1.0으로 너무 크게 즉 너그럽게 잡혀있기 때문입니다. 이 값을 0.4나 0.3 정도로 줄여보면, Chi-Square 값은 더 커지면서 유사성 확률이 점차 줄어드는 것을 확인할 수 있습니다. 물론 이전 게시물에서도 언급했듯이 에러분포 표준편차가 실제로 어느 정도인지는 데이터의 특성 및 여러가지를 고려해서 작업자가 직접 결정해야 합니다.

지금까지 데이터 포인트들을 비선형 함수로 근사하는 작업을 IDL 내장함수인 LMFIT 그리고 외부 라이브러리 함수인 MPFITFUN을 사용하여 수행하는 두 종류의 과정들을 살펴보았습니다. 이미 우리가 느꼈듯이 사용 방법의 편의성의 측면에서보면 아무래도 LMFIT 보다는 MPFITFUN이 더 나아보이는 것이 사실입니다. 수식의 계수별 편미분을 따로 구하는 것이 매우 번거로운 작업일 수 밖에 없기 때문입니다. 따라서 이러한 차이점을 감안해서 적절한 방법을 선택하면 될 것 같습니다.

그리고 이러한 근사 작업에 있어서 한가지 더 유념해야 할 것은, 계수들의 초기 추정값이 최종 계산값과 너무 차이가 나면 계산 결과가 제대로 나오지 않을 수도 있다는 것입니다. 예를 들어 최종 계산된 계수값은 2.25와 5.70인데, 초기 추정치를 1.0과 2.0으로 줄 경우에는 수렴을 제대로 하지 못해서 계산 결과가 제대로 나오지 않습니다. 그나마 좀 더 근접한 값을 초기 추정치로 주어야 결과가 제대로 나옵니다. 이런 불편함을 극복하는 방법이 따로 있을 수도 있겠지만, 어쨌든 비선형 함수로의 근사라는 작업 자체가 좀 이와 같은 민감성이 있는 작업이라는 특성 때문일지도 모르겠습니다. 자 그러면 비선형 함수의 근사에 관한 내용은 이 정도로 마무리하겠습니다.