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

오늘은 데이터를 비선형 함수(Non-Linear Function)의 형태로 근사(Fitting)하는 작업을 IDL에서 어떻게 하는가에 관하여 다뤄보고자 합니다. 데이터 포인트들이 주어졌을 때 그 패턴을 특정 함수에 맞추는 작업은 과학기술 분야에서 매우 자주 하게 되는 일입니다. 물론 그러한 작업들 중에서는 선형 함수(Linear Function)에 맞추는 경우가 더 일반적이긴 합니다. 1차, 2차 또는 다차 다항식에 맞추는 것으로도 원하는 결과를 어느 정도는 산출할 수 있습니다. 하지만 맞춰야 하는 대상 모델이 비선형 함수의 형태를 갖는 경우들도 종종 있습니다. 이러한 작업을 IDL에서 어떤 방식으로 처리하느냐가 바로 오늘의 주제입니다.

일단 큰 틀에서 보면 두가지 방법이 있습니다. 첫번째는 IDL에서 기본적으로 지원되는 내장함수를 사용하는 방법이고, 두번째는 외부 라이브러리를 이용하는 방법입니다. 첫번째 방법에 해당되는 내장함수는 바로 LMFIT 함수입니다. 그리고 두번째 방법에서 언급한 외부 라이브러리는 바로 MPFIT 라이브러리인데, 이 라이브러리 안에서 지원되는 여러 루틴들 중에서도 MPFITFUN이라는 함수가 바로 그러한 역할을 합니다. 그래서 오늘은 LMFIT와 MPFITFUN 두 함수를 사용하여 데이터 포인트들을 비선형 함수에 근사하는 실습을 해보고자 합니다. 다만 아까 언급했듯이 LMFIT 함수는 IDL에서 기본 지원되기 때문에 그냥 바로 사용할 수 있습니다. 그러나 MPFITFUN 함수는 IDL에 기본 내장된 기능이 아니기 때문에, 사용을 위해서는 별도로 MPFIT 라이브러리를 받아서 설치하는 작업이 반드시 선행되어야 합니다. 이 라이브러리는 아래 웹페이지에서 받을 수 있습니다.

http://www.physics.wisc.edu/~craigm/idl/fitting.html

이 웹페이지의 중간쯤을 보면 mpfit이라는 이름의 tar 및 zip 파일을 받을 수 있습니다. 이 파일을 받아서 압축을 푼 다음 그 디렉토리를 IDL의 경로에 추가해야 MPFIT 라이브러리의 여러 기능들을 바로바로 사용할 수 있습니다. IDL에서 외부 라이브러리를 설정하는 방법에 관해서는 제가 예전에 올렸던 이 게시물의 내용을 참조하시면 됩니다. 그러면 이러한 준비가 된 상태로 가정하고 내용을 이어나가겠습니다. 먼저 예제로 사용할 데이터는 그냥 가상으로 다음과 같이 만들어보았습니다. 다만 가상 데이터를 그냥 아무렇게나 만든 것은 아니고, 대략적으로는 함수의 패턴을 가지면서 흩어짐(Spread)이 어느 정도 있도록 하였습니다.

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)

이 가상 데이터는 전반적으로는 Y=2.3*SIN(5.7*X)의 형태를 따라가지만 난수 생성 함수인 RANDOMU 함수를 사용하여 인위적으로 오차를 추가하는 방식으로 생성하였습니다. 데이터 포인트들의 총 갯수는 101개입니다. 그래서 다음과 같이 데이터 포인트들만 플롯으로 먼저 표출해서 그 형태를 봅시다.

yrn = [-3, 3]

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

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

이렇게 하면 그 모습은 다음 그림과 같습니다.

그래서 이런 형태를 갖는 데이터 포인트들이 주어진 상태에서 그 패턴을 비선형 함수에 맞춰보는 작업을 해봅시다. 이를 위해서는 어떤 형태의 함수가 되어야 할 것인지를 먼저 결정해야 합니다. 그래서 Y=P0*SIN(P1*X)의 형태를 갖는 함수로 맞추기로 하겠습니다. 결국 가장 적절한 P0, P1 계수값들을 알아내는 과정이 될 것입니다. 사실은 이 작업이 제대로 수행되었을 경우의 결과는 우리는 이미 충분히 예상 가능합니다. 당연히 P0는 2.3에 가까운 값이 될 것이고 P1은 5.7에 가까운 값이 될 것입니다. 하지만 그러한 결론을 얻기 위하여 어떤 방식으로 작업을 해야 할 것인지를 LMFIT을 사용하는 방식과 MPFITFUN을 사용하는 방식으로 나눠서 지금부터 살펴보겠습니다. 사실 제가 살짝 스포일러를 말씀드리면, LMFIT을 사용하는 방법보다는 MPFITFUN을 사용하는 방법이 약간 더 쉽습니다. 좀 더 정확히 얘기한다면 LMFIT을 사용하는 방법은 뭔가 좀 더 번거로운 측면이 있습니다. 그 이유를 먼저 IDL의 내장함수인 LMFIT을 사용하는 방법부터 보면서 알아보겠습니다.

IDL 도움말에서 LMFIT 함수에 관한 내용을 보면 사용 문법은 다음과 같습니다.

Result = LMFIT( X, Y, A [, ALPHA=variable] [, CHISQ=variable] [, CONVERGENCE=variable] [, COVAR=variable] [, /DOUBLE] [, FITA=vector] [, FUNCTION_NAME=string] [, ITER=variable] [, ITMAX=value] [, ITMIN=value] [, MEASURE_ERRORS=vector] [, SIGMA=variable] [, TOL=value] )

보시다시피 일단 필수인자가 세 개가 무조건 필요합니다. 처음에 나오는 X, Y는 당연히 데이터 포인트들에 해당되기 때문에 위의 예제코드에서 정의된 X, Y와 동일합니다. 그런데 A라고 되어 있는 세번째 인자는 "initial estimate for each coefficient"라고 나와 있습니다. 즉 계수들에 대한 초기 추정값들로 구성된 배열이란 얘기이고, 결국 계수들에 대한 초기 추정값들을 우리가 직접 넣어줘야 한다는 얘기입니다. 이렇게 계수 추정값들을 처음에 주는 것이 이러한 작업에서는 반드시 필요합니다. 이것은 MPFITFUN에서도 마찬가지입니다. 사실 우리는 이미 정답을 알고 있습니다. 그래서 아마도 초기 추정값들을 [2.3, 5.7]로 부여하면 당연히 작업은 한방에 끝날 것입니다. 하지만 추정값이므로 그냥 적당히 유사한 값들을 주기로 합시다. 사실 나중에 확인해볼 수 있겠지만, 이 초기 추정값들을 어떻게 설정하느냐가 최종 결과의산출에 꽤 민감하게 작용을 합니다.

그리고 키워드들 중에 주목해야 할 것이 두가지 정도가 있는데 바로 MEASURE_ERRORS와 FUNCTION_NAME 키워드들입니다. 먼저 MEASURE_ERRORS 키워드에는 각 데이터 포인트마다의 에러 분포의 표준편차값들로 구성된 배열을 부여해야 합니다. 이 값이 크면 에러 폭을 여유있게 주는 것이고 작으면 에러 폭을 깐깐하게 주는 것입니다. 그 적정한 값은 데이터마다 다릅니다. 이 값은 근사 결과과 원본 데이터와 어느 정도 유사한가를 정량화하는 Chi-Square라는 값을 산출하는 근거가 됩니다. 그리고 FUNCTION_NAME 키워드가 있는데, 이 키워드에는 어떤 사용자 정의 함수의 이름을 문자값으로 부여해야 합니다. 그런데 그 "사용자 정의 함수"라는 것은 도대체 무엇이냐 이것이 관건입니다. 하지만 일단 LMFIT을 사용한 내용을 먼저 본 다음에 얘기를 이어나가겠습니다. 그 내용은 다음과 같습니다.

errors = MAKE_ARRAY(101, VALUE=0.5)

start = [2.0, 5.0]

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

  ITER=n_iter, CHISQ=chisq, /DOUBLE)

여기서는 계수들의 초기 추정값들을 start라는 배열에 담아서 LMFIT 함수의 세번째 인수로 부여하였습니다. 그리고 에러분포의 표준편차는 101개의 데이터 포인트들 모두에 대하여 0.5로 가정을 하고, 0.5라는 값 101개로 구성된 errors라는 배열을 만들어서 사용하였습니다. 그런데 FUNCTION_NAME 키워드를 보면 'my_nl_funct'라는 이름이 보입니다. 사실 이 시점에서는 my_nl_funct라는 이름의 함수 부프로그램을 다음과 같이 따로 만들어야 합니다. 이것은 사용자 정의 함수로서 현재의 작업에서 반드시 필요합니다.

FUNCTION my_nl_funct, X, P

result = [P[0]*SIN(P[1]*X), SIN(P[1]*X), P[0]*X*COS(P[1]*X)]

RETURN, result

END

이 내용에 관해서는 좀 있다 설명을 하겠습니다. 일단 지금까지의 내용을 IDL 프로그램으로 저장해두는 것이 좋겠습니다. 전체적으로는 주프로그램인 test_nonlinear_fitting과 부프로그램인 my_nl_funct이 공존하는 다음과 같은 내용이 되어야 합니다. 그리고 이 전체 내용을 그대로 test_nonlinear_fitting.pro라는 프로그램 파일로 저장해두시기 바랍니다.

FUNCTION my_nl_funct, X, P

result = [P[0]*SIN(P[1]*X), SIN(P[1]*X), P[0]*X*COS(P[1]*X)]

RETURN, result

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 = 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

그러면 my_nl_funct라는 함수 부프로그램의 내용으로 다시 돌아오겠습니다. 이 내용을 보면 X, P 두 개의 인자를 받는 사용자 정의 함수임을 알 수 있습니다. 이 함수는 X 인자를 통해서 데이터 포인트의 X값 배열인 X를 그대로 받게 됩니다. 그리고 P 인자는 현재 우리가 근사하고자 하는 비선형 함수의 계수값들로 구성된 배열에 대응됩니다. 앞서 언급했듯이 두 개의 계수들로 정의되기 때문에 이 값들로 구성된 배열이 P 인자에 대응된다고 보면 됩니다. 이와 같이 my_nl_funct라는 사용자 정의 함수는 근사하고자 하는 함수의 형태를 가지면서 X와 P 인자들을 통하여 데이터 X값 배열 및 계수값 배열을 인자로 받을 수 있도록 정의되어야 합니다. 그러면 자연스럽게 생각해보면이 함수는 내부적으로 Y에 해당되는 값들을 계산하는 형태로 정의되어야 할 것으로 생각할 수 있습니다. 그렇게 따지면 사실 다음과 같은 형태가 되는 것이 더 자연스럽게 보일 수도 있습니다.

FUNCTION my_nl_funct, X, P

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

RETURN, result

END

즉 이와 같이 함수 내부에서 X와 P를 이용하여 대상 함수의 형태에 맞는 연산을 하도록 정의되는 것이 맞지 않겠냐 하는 생각이 들지 않으신지요? 그런데 실제로는 LMFIT를 사용하려면 근사할 함수가 이런 식으로 정의되면 안됩니다. 원래대로 다음과 같이 정의되어야 합니다.

FUNCTION my_nl_funct, X, P

result = [P[0]*SIN(P[1]*X), SIN(P[1]*X), P[0]*X*COS(P[1]*X)]

RETURN, result

END

왜냐하면 LMFIT 함수를 온전하게 사용하려면 사용자 정의 함수를 반드시 이런 방식으로 정의해야 하기 때문입니다. 그러면 위의 내용에서 result는 왜 저렇게 세 개의 값들로 구성된 배열이 되어야 할까요? 사실 저 세 개의 값들 중 첫번째 것은 Y에 해당되는 값으로 쉽게 알아볼 수 있지만, 두번째 및 세번째는 각 계수(P[0] 및 P[1])에 대한 편미분(Partial Derivative) 결과입니다. 이런 방식으로 정의된 항목들로 구성된 배열이 되어야 합니다. 그 이유는 사실 LMFIT 함수가 내부적으로 사용하는 알고리즘의 특성으로 인한 것이다 정도만 언급드릴 수 있을 것 같습니다. 어쨌든 LMFIT 함수를 사용하려면 이렇게 해야 합니다. 일단 "편미분"이란 용어가 등장하는 것만으로도 벌써 골치가 아파옵니다. 그나마 다시 정리해보면 다음과 같습니다.

result 배열의 첫번째 원소 : P[0]*SIN(P[1]*X)

result 배열의 두번째 원소 : P[0]*SIN(P[1]*X)를 P[0]에 대하여 편미분한 SIN(P[1]*X)

result 배열의 세번째 원소 : P[0]*SIN(P[1]*X)를 P[1]에 대하여 편미분한 P[0]*X*COS(P[1]*X)]

그러면 일단 이렇다고 치고 다음 과정으로 넘어가보겠습니다. LMFIT 함수가 이런 방식으로 실행이 되면 결과가 나옵니다. 먼저 주프로그램에서 LMFIT 함수가 사용된 라인에서 좌변을 보면 coefs라는 항목이 정의되어 있는데요. 이 항목 원본 데이터 포인트의 갯수와 동일한 101개의 근사값들로 구성된 배열로 얻어집니다. 그래서 coefs라는 이름이 좀 무색하긴 하지만 하여간 그렇습니다. 그런데 이러한 결과도 필요하지만 더 중요한 것은 함수의 계수값들입니다. 이건 어디에 있을까요? 사실은 처음에 start라는 배열로 줬던 계수 초기 추정값들이 결국 LMFIT의 작업이 끝나면서 계수 산출값들로 대체됩니다. 그래서 주프로그램의 내용에서 coefs에 대하여 HELP 및 PRINT 명령을 사용하여 그 형태 및 값들을 출력해서 볼 수 있도록 한 것인데, 실제로 101개의 실수형(2배 정밀도) 근사값들이 출력되는 것을 확인할 수 있습니다. 그리고 start의 값들을 출력해본 결과는 다음과 같습니다.

       2.2525569       5.7077310

이것이 바로 P0, P1에 해당되는 계수 산출값들입니다. 애초에 가상 데이터를 만들 때 사용했던 계수값들이 2.3과 5.7이었던 것을 감안하면 거의 근접한 결과를 얻었음을 알 수 있습니다. 그리고 coefs의 101개의 값들을 원본 데이터 포인트들과 중첩하여 그려보는 과정이 마지막에 들어가 있습니다. 그 결과는 다음 그림과 같습니다.

그리고 LMFIT 함수에서 ITER와 CHISQ라는 키워드도 사용된 것을 볼 수 있습니다. ITER 키워드는 실제 계산 과정에서 iteration이 수행된 횟수값을 돌려주는 역할을 합니다. 그리고 CHISQ는 Chi-Square 값을 돌려주는 역할을 합니다. 위 내용에서는 n_iter, chisq 변수를 통해서 이러한 값들을 획득하였는데 출력된 두 값을 보면 다음과 같습니다.

       39.274361

           7

즉 총 7회에 걸쳐 iteration이 수행되었고, Chi-Square값은 약 39.27이 나왔습니다. Chi-Square의 값은 근사된 함수와 원본 데이터 사이의 유사성이 어느 정도인지를 나타내는 유사성 확률을 계산하는데 사용될 수 있습니다. 이것은 CHISQR_PDF라는 함수를 사용하여 계산합니다. 즉 CHISQR_PDF 함수에 Chi-Square 값과 DoF(Degree of Freedom) 값을 인자로 주면 됩니다. DoF는 데이터 포인트의 총 갯수에서 계수의 갯수를 뺀 값입니다. 지금의 경우라면 101-2=99가 됩니다. 그래서 유사성 확률은 다음과 같이 계산할 수 있습니다.

PRINT, 1-CHISQR_PDF(39.27, 101-2)

이렇게 계산된 유사성 확률은 1.00이 나옵니다. 즉 100%란 얘기입니다. 너무 "잘" 나온 감이 분명히 있습니다. 그 이유는 errors에서 정의된 에러분포 표준편차 값이 0.5이기 때문입니다. 다시 말하면 에러의 폭을 너무 크게(후하게) 잡았기 때문입니다. 이 값을 좀 더 줄여보면 Chi-Square값이 더 크게 나오면서 유사성 확률은 줄어듭니다. 물론 에러분포의 표준편차를 어느 정도로 잡아야 할 것이냐는 데이터의 특성에 따라 다릅니다.

어쨌든 위의 그림을 보면 결과가 그래도 우리가 기대했던 만큼은 나온 것 같습니다. 그런데 여기까지 오는 과정을 보면 다른 것은 뭐 그렇다고 쳐도, 근사할 함수의 형태에 해당되는 사용자 정의 함수를 만드는데 있어서 각 인자별 편미분을 산출해야 하는 등 과정 자체가 너무 번거롭고 복잡합니다. 그나마 제가 여기서는 나름 간단한 형태의 비선형 함수를 사용하긴 했지만, 실전에서 직면하게 될 것들은 아마도 SIN, COS, Expential, Log 등이 난무하는 더 무시무시한 경우가 될 가능성이 높습니다. 더구나 산출해야 할 계수들이 더 많아질 가능성도 높습니다. 이런 상황에서 대상 비선형 함수에 대한 각 계수별 편미분을 일일이 구한다는 것은 굉장히 번거로운 일입니다. 저는 이 글 작성하면서 SIN의 미분이 COS이란 것도 갑자기 기억이 나지않아서(분명히 고등학교 시절에 배웠는데도) 검색해서 찾았는데, 더 복잡한 경우는 아예 생각하기도 싫어집니다.

사실 IDL에서 비선형이 아닌 선형 함수 근사에 있어서는 LINFIT, POLY_FIT, REGRESS 등 사용하기 편하고 좋은 것들이 많습니다. 하지만 비선형 함수로 가면 얘기가 좀 달라집니다. 앞서 언급한 사용 방법상의 번거로움 때문에 아무래도 IDL의 LMFIT 함수는 사용을 적극적으로 추천해드리기가 약간 애매한 측면이 있습니다. LMFIT 외에도 유사한 느낌의 CURVEFIT이라는 함수도 있는데, 이것도 역시 비슷한 방식으로 사용해야 합니다. 그래서 예전부터 IDL 사용자들 사이에서는 비선형 함수 근사에 있어서 사용법이 좀 더 간편한 다른 대체 수단에 대한 수요가 많았는데, 바로 MPFIT 라이브러리가 그러한 역할을 상당 부분 담당해왔던 것도 사실입니다. 그래서 오늘은 일단 여기까지 하고, 다음 회차에서는 이 MPFIT 라이브러리의 MPFITFUN 함수를 이용하여, 앞서 실습했던 동일한 예제 데이터에 대하여 비선형 근사 결과를 얻어보는 과정을 살펴보고 기타 관련 얘기들도 함께 이어나가도록 하겠습니다.