Main Content

이 번역 페이지는 최신 내용을 담고 있지 않습니다. 최신 내용을 영문으로 보려면 여기를 클릭하십시오.

최소제곱 피팅

소개

Curve Fitting Toolbox™는 데이터를 피팅할 때 최소제곱 방법을 사용합니다. 피팅에는 하나 이상의 계수로 응답 변수 데이터를 예측 변수 데이터에 연관 짓는 모수적 모델이 필요합니다. 피팅 과정의 결과는 모델 계수의 추정값입니다.

계수 추정값을 얻기 위해, 최소제곱 방법은 잔차의 제곱합을 최소화합니다. i번째 데이터 점에 대한 잔차 ri는 관측된 응답 변수 값 yi와 피팅된 응답 변수 값 ŷi 사이의 차이로 정의되며, 데이터와 관련 있는 오차입니다.

ri=yiy^iresidual=datafit

잔차의 제곱합은 다음과 같이 지정됩니다.

S=i=1nri2=i=1n(yiy^i)2

여기서 n은 피팅에 포함된 데이터 점의 개수이고 S는 오차 제곱합 추정값입니다. 다음과 같은 유형의 최소제곱 피팅이 지원됩니다.

  • 선형 최소제곱

  • 가중 선형 최소제곱

  • 로버스트 최소제곱

  • 비선형 최소제곱

오차 분포

불규칙 변동을 포함하는 데이터를 피팅할 때는 일반적으로 오차에 대해 다음과 같은 두 가지 중요한 가정이 이루어집니다.

  • 오차는 응답 변수 데이터에만 존재하고 예측 변수 데이터에는 존재하지 않는다.

  • 오차는 무작위적이며 평균이 0이고 분산이 상수 σ2인 정규(가우스) 분포를 따른다.

두 번째 가정은 다음과 같이 표현되기도 합니다.

errorN(0,σ2)

오차가 정규분포를 따른다고 가정하는 이유는 여러 측정량의 분포에서 정규분포가 적절한 근사를 제공하는 경우가 많기 때문입니다. 최소제곱 피팅 방법은 파라미터 추정값을 계산할 때 정규분포된 오차를 가정하지는 않지만, 그래도 이 방법은 극값을 갖는 다량의 랜덤 오차를 포함하지 않는 데이터에서 가장 좋은 결과를 나타냅니다. 정규분포는 극단적인 랜덤 오차가 흔치 않은 확률 분포 중 하나입니다. 그러나 신뢰한계 및 예측한계와 같은 통계적 결과가 타당성을 가지려면 오차가 반드시 정규분포를 이루어야 합니다.

오차의 평균이 0이면 오차는 순수하게 무작위적입니다. 평균이 0이 아니면 모델이 데이터에 대해 적절한 선택이 아니거나 오차가 순수하게 무작위적이지 않고 계통 오차가 포함된 것일 수 있습니다.

데이터의 분산이 상수라는 것은 오차의 “산포”가 일정함을 의미합니다. 동일한 분산을 갖는 데이터는 동일한 질을 갖는다고도 합니다.

가중 최소제곱 회귀에서는 랜덤 오차가 상수 분산을 갖는다는 가정이 내포되어 있지 않습니다. 대신, 피팅 절차에서 제공되는 가중치가 데이터에 존재하는 질의 서로 다른 수준을 올바르게 나타낸다고 가정합니다. 가중치는 각 데이터 점이 피팅된 계수의 추정값에 미치는 영향도를 적절한 수준으로 조정하는 데 사용됩니다.

선형 최소제곱

Curve Fitting Toolbox는 선형 최소제곱 방법을 사용하여 선형 모델을 데이터에 피팅합니다. 선형 모델은 계수가 선형인 방정식으로 정의됩니다. 예를 들어, 다항식은 선형이지만 가우스 분포 함수는 선형이 아닙니다. 선형 최소제곱 피팅 과정을 살펴보기 위해 1차 다항식으로 모델링할 수 있는 n개의 데이터 점이 있다고 가정하겠습니다.

y=p1x+p2

알려지지 않은 계수 p1과 p2에 대해 이 방정식을 풀기 위해 S를 2개의 알려지지 않은 계수로 이루어진 n개의 선형 연립방정식으로 표현합니다. n이 알려지지 않은 계수의 개수보다 크면 연립방정식이 과결정인 것입니다.

S=i=1n(yi(p1xi+p2))2

최소제곱 피팅 과정에서는 잔차의 제곱합을 최소화하므로, 계수는 S를 각 파라미터에 대해 미분하고 그 결과를 0과 같다고 놓음으로써 결정됩니다.

Sp1=2i=1nxi(yi(p1xi+p2))=0Sp2=2i=1n(yi(p1xi+p2))=0

실제 파라미터의 추정값은 보통 b로 표현됩니다. p1과 p2를 b1과 b2로 바꾸면 위의 방정식이 다음과 같이 됩니다.

xi(yi(b1xi+b2))=0    (yi(b1xi+b2))=0

여기서 합은 i = 1부터 n까지 진행됩니다. 정규 방정식은 다음과 같이 정의됩니다.

b1xi2+b2xi=xiyib1xi+nb2=yi

b1에 대해 풀면 다음과 같습니다.

b1=nxiyixiyinxi2(xi)2

b1 값을 사용하여 b2에 대해 풀면 다음과 같습니다.

b2=1n(yib1xi)

앞에서 보았듯이 계수 p1과 p2를 추정하는 데는 몇 개의 간단한 계산만 필요합니다. 이 예제를 보다 고차의 다항식으로 확장하는 것은 조금 지루하지만 간단한 작업입니다. 각 일차항에 대해 정규 방정식 하나씩을 모델에 추가하기만 하면 됩니다.

선형 모델을 행렬 형식으로 나타내면 다음 수식이 됩니다.

y = Xβ + ε

여기서

  • y는 응답 변수로 구성된 n×1 벡터입니다.

  • β는 계수로 구성된 m×1 벡터입니다.

  • X는 모델에 대한 n×m 설계 행렬입니다.

  • ε은 오차로 구성된 n×1 벡터입니다.

1차 다항식의 경우, 2개의 알려지지 않은 계수로 이루어진 n개의 방정식은 y, X 및 β로 다음과 같이 표현됩니다.

[y1y2y3...yn]=[x11x21x31...xn1]×[p1p2]

문제에 대한 최소제곱해는 알려지지 않은 계수 벡터 β를 추정하는 벡터 b입니다. 정규 방정식은 다음과 같이 지정됩니다.

(XTX)b = XTy

여기서 XT는 설계 행렬 X의 전치입니다. b에 대해 풀면 다음과 같습니다.

b = (XTX)–1 XTy

알려지지 않은 계수에 대해 선형 연립방정식을 풀려면 MATLAB® 백슬래시 연산자(mldivide)를 사용하십시오. XTX의 역행렬을 구하면 허용되지 않는 반올림 오차가 발생할 수 있으므로 백슬래시 연산자는 피벗 연산을 사용한 QR 분해를 사용하는데, 이는 수치적으로 매우 안정적인 알고리즘입니다. 백슬래시 연산자와 QR 분해에 대한 자세한 내용은 산술 연산 항목을 참조하십시오.

b를 다시 모델식에 적용하여 예측된 응답 변수 값 ŷ를 얻을 수 있습니다.

ŷ = Xb = Hy

H = X(XTX)–1 XT

문자 위의 해트(곡절) 기호는 파라미터에 대한 추정값 또는 모델에서 도출한 예측값을 나타냅니다. 투영 행렬 H는 y 위에 해트를 적용하므로 해트 행렬이라고 합니다.

잔차는 다음과 같이 지정됩니다.

r = y – ŷ = (1–H)y

가중 최소제곱

일반적으로 응답 변수 데이터는 동일한 질을 가지며 따라서 상수 분산을 갖는다고 가정됩니다. 이 가정이 위배될 경우 피팅이 질 낮은 데이터에 의해 과하게 영향을 받을 수 있습니다. 피팅을 개선하려면 피팅 과정에 추가적인 스케일링 인자(가중치)가 포함되는 가중 최소제곱 회귀를 사용하면 됩니다. 가중 최소제곱 회귀는 다음과 같은 오차 추정값을 최소화합니다.

s=i=1nwi(yiy^i)2

여기서 wi는 가중치입니다. 가중치는 각 응답 변수 값이 최종 파라미터 추정값에 얼마만큼의 영향을 주는지 결정합니다. 질 높은 데이터 점은 질 낮은 데이터 점보다 피팅에 더 큰 영향을 줍니다. 가중치가 알려져 있거나 가중치가 특정한 형태를 따른다는 근거가 있는 경우 데이터에 가중치를 적용하는 것이 좋습니다.

가중치는 파라미터 추정값 b에 대한 표현식을 다음과 같은 방식으로 수정합니다.

b=β^=(XTWX)1XTWy

여기서 W는 가중치 행렬 w의 대각선 요소로 지정됩니다.

데이터를 피팅하고 잔차를 플로팅하면 많은 경우 분산이 상수가 아닌지 여부를 확인할 수 있습니다. 아래 플롯에서 데이터는 다양한 질의 반복 실험 데이터를 포함하며 피팅은 올바르다고 가정됩니다. 잔차의 플롯을 보면, 작은 예측 변수 값이 큰 예측 변수 값보다 응답 변수 값의 더 큰 산포를 생성하는 “깔때기” 모양에서 데이터의 질이 낮은 것을 알 수 있습니다.

사용자가 제공하는 가중치는 응답 변수 분산을 상수 값으로 변환해야 합니다. 데이터의 측정 오차의 분산을 알고 있다면 가중치는 다음과 같이 지정됩니다.

wi=1/σi2

또는 각 데이터 점에 대해 오차 변수 추정값만 알고 있다면 실제 분산 대신 추정값을 사용하는 것으로도 충분한 경우가 많습니다. 분산을 모른다면 상대적인 척도로 가중치를 지정하는 것으로 충분합니다. 가중치가 지정된 경우에도 전체 분산 항이 추정된다는 점에 주의하십시오. 여기서 가중치는 피팅의 각 점에 대한 상대 가중치를 정의하지만, 각 점의 정확한 분산을 지정하는 데 사용되지는 않습니다.

예를 들어, 각 데이터 점이 여러 독립 측정값의 평균이라면 측정값의 수치를 가중치로 사용하는 것이 적절할 수 있습니다.

로버스트 최소제곱

일반적으로 응답 변수 오차는 정규분포를 따르고 극값은 드물다고 가정됩니다. 그렇기는 해도 이상값이라 부르는 극값이 발생하기는 합니다.

최소제곱 피팅의 주요 장점은 이상값에 대한 민감도입니다. 잔차를 제곱하면 극단에 있는 데이터 점의 영향이 확대되기 때문에 이상값은 피팅에 큰 영향을 줍니다. 이상값의 영향을 최소화하려면 로버스트 최소제곱 회귀를 사용하여 데이터를 피팅하면 됩니다. 이 툴박스는 다음과 같은 두 가지 로버스트 회귀 방법을 제공합니다.

  • 최소 절대 잔차(LAR) — LAR 방법은 잔차의 차이의 제곱이 아니라 잔차의 차이의 절대값을 최소화하는 곡선을 찾습니다. 따라서 극값이 피팅에 보다 적은 영향을 줍니다.

  • 겹제곱(bisquare) 가중치 — 이 방법은 가중 제곱합을 최소화합니다. 여기서 각 데이터 점에 부여되는 가중치는 피팅된 선으로부터 그 점이 얼마나 멀리 떨어져 있는지에 따라 결정됩니다. 선에 가까이 있는 점에는 온전한 가중치가 부여됩니다. 선에서 보다 멀리 있는 점에는 줄어든 가중치가 부여됩니다. 우연에 의해 기대되는 것보다 더 선에서 멀리 있는 점에는 가중치 0이 부여됩니다.

    대부분의 경우, 겹제곱 가중치 방법이 일반적인 최소제곱 방법을 써서 데이터 대부분을 피팅하는 곡선을 동시에 찾으면서 이상값의 영향을 최소화하기 때문에 LAR보다 선호됩니다.

겹제곱 가중치를 사용하는 로버스트 피팅은 반복 재가중 최소제곱 알고리즘을 사용하며, 다음 절차를 따릅니다.

  1. 가중 최소제곱으로 모델을 피팅합니다.

  2. 조정된 잔차를 계산하고 표준화합니다. 조정된 잔차는 다음과 같이 지정됩니다.

    radj=ri1hi

    ri는 일반적인 최소제곱 잔차이고 hi는 최소제곱 피팅에 큰 영향을 주는, 높은 지렛대값을 갖는 데이터 점의 가중치를 줄여서 잔차를 조정하는 지렛대값입니다. 표준화된 조정된 잔차는 다음과 같이 지정됩니다.

    u=radjKs

    K는 값이 4.685인 조율 상수이고, s는 MAD/0.6745로 지정되는 로버스트 표준편차입니다. 여기서 MAD는 잔차의 중앙값 절대편차(Median Absolute Deviation)입니다.

  3. 로버스트 가중치를 u의 함수로 계산합니다. 겹제곱 가중치는 다음과 같이 지정됩니다.

    wi={(1(ui)2)2|ui|<10|ui|1

    사용자 고유의 회귀 가중 벡터를 제공할 경우 최종 가중치는 로버스트 가중치와 회귀 가중치의 곱이 됩니다.

  4. 피팅이 수렴하면 작업이 완료됩니다. 수렴하지 않을 경우, 첫 번째 단계로 돌아가서 피팅 절차의 다음 반복을 수행합니다.

아래에 표시된 플롯에서는 일반적인 선형 피팅을 겹제곱 가중치를 사용한 로버스트 피팅과 비교합니다. 로버스트 피팅이 데이터의 대부분을 따르면서 이상값의 영향을 강하게 받지 않음을 알 수 있습니다.

로버스트 회귀를 사용하여 이상값의 영향을 최소화하는 대신, 피팅에서 제외할 데이터 점을 지정해 주는 방법도 가능합니다. 자세한 내용은 이상값 제거하기 항목을 참조하십시오.

비선형 최소제곱

Curve Fitting Toolbox는 비선형 최소제곱 공식을 사용하여 비선형 모델을 데이터에 피팅합니다. 비선형 모델은 계수가 비선형이거나 선형과 비선형이 혼합된 방정식으로 정의됩니다. 예를 들어, 가우스 분포 함수, 다항식의 비, 멱함수는 모두 비선형입니다.

비선형 모델을 행렬 형식으로 나타내면 다음 수식이 됩니다.

y = f(X,β) + ε

여기서

  • y는 응답 변수로 구성된 n×1 벡터입니다.

  • f는 β와 X의 함수입니다.

  • β는 계수로 구성된 m×1 벡터입니다.

  • X는 모델에 대한 n×m 설계 행렬입니다.

  • ε은 오차로 구성된 n×1 벡터입니다.

비선형 모델은 간단한 행렬 기법을 사용하여 계수를 추정할 수 없기 때문에 선형 모델보다 피팅하기가 더 어렵습니다. 비선형 모델에는 다음 단계를 따르는 반복법이 필요합니다.

  1. 각 계수에 대한 초기 추정값으로 시작합니다. 일부 비선형 모델의 경우에는 합리적인 시작값을 생성하는 발견법적 접근 방식이 제공됩니다. 그 밖의 모델의 경우에는 구간 [0,1] 내의 임의 값이 제공됩니다.

  2. 현재 계수 집합에 대해 피팅된 곡선을 생성합니다. 피팅된 응답 변수 값 ŷ는 다음과 같이 지정됩니다.

    ŷ = f(X,b)

    이 값을 구하려면 f(X,b)의 야코비 행렬을 계산해야 합니다. 야코비 행렬은 계수에 대한 편도함수로 구성된 행렬로 정의됩니다.

  3. 계수를 조정하고 피팅이 개선되는지 확인합니다. 조정의 방향과 크기는 피팅 알고리즘에 따라 달라집니다. 이 툴박스는 다음과 같은 알고리즘을 제공합니다.

    • Trust-Region — 디폴트 알고리즘으로, 계수 제약 조건을 지정한 경우 반드시 사용해야 합니다. 이 알고리즘은 널리 사용되는 Levenberg-Marquardt 알고리즘을 개선한 것으로, 여타 알고리즘보다 까다로운 비선형 문제를 더 효율적으로 풉니다.

    • Levenberg-Marquardt — 이 알고리즘은 오랜 기간 사용되어 왔으며 다양한 비선형 모델과 시작값에서 효과가 있다는 사실이 입증되었습니다. Trust-Region 알고리즘이 적절한 피팅을 생성하지 못하고 계수 제약 조건이 없는 경우 Levenberg-Marquardt 알고리즘을 시도해 봐야 합니다.

  4. 2단계로 돌아가서 피팅이 지정된 수렴 조건에 도달할 때까지 이 과정을 반복합니다.

비선형 모델에 대해 가중치와 로버스트 피팅을 사용할 수 있습니다. 그에 따라 피팅 과정이 수정됩니다.

근사 과정의 특성 때문에 모든 비선형 모델, 데이터 세트 및 시작점에 대해 실패 없이 사용할 수 있는 알고리즘은 없습니다. 따라서 디폴트 시작점, 알고리즘 및 수렴 조건을 사용하여 합리적인 피팅을 달성하지 못했다면 다른 옵션으로 시험해 보아야 합니다. 디폴트 옵션을 수정하는 방법에 대한 설명은 피팅 옵션 및 최적화된 시작점 지정하기 항목을 참조하십시오. 비선형 모델은 시작점에 특히 민감할 수 있으므로 시작점이 가장 먼저 수정하는 피팅 옵션이 되어야 합니다.

로버스트 피팅

이 예제에서는 이상값을 제외하는 것의 효과와 로버스트 피팅의 효과를 비교하는 방법을 보여줍니다. 이 예제에서는 모델에서 1.5 표준편차보다 큰 임의의 거리에 있는 이상값을 제외하는 방법을 보여줍니다. 그런 다음 이상값을 제거하는 것과 이상값에 보다 적은 가중치를 부여하는 로버스트 피팅을 지정하는 것을 비교합니다.

기준 정현파 신호를 만듭니다.

xdata = (0:0.1:2*pi)'; 
y0 = sin(xdata);

신호에 일정하지 않은 분산을 갖는 잡음을 추가합니다.

% Response-dependent Gaussian noise
gnoise = y0.*randn(size(y0));

% Salt-and-pepper noise
spnoise = zeros(size(y0)); 
p = randperm(length(y0));
sppoints = p(1:round(length(p)/5));
spnoise(sppoints) = 5*sign(y0(sppoints));

ydata = y0 + gnoise + spnoise;

잡음이 있는 데이터를 기준 정현파 모델로 피팅하고, 3개의 출력 인수를 지정하여 잔차를 포함한 피팅 정보를 얻습니다.

f = fittype('a*sin(b*x)'); 
[fit1,gof,fitinfo] = fit(xdata,ydata,f,'StartPoint',[1 1]);

fitinfo 구조체의 정보를 검토합니다.

fitinfo
fitinfo = struct with fields:
           numobs: 63
         numparam: 2
        residuals: [63x1 double]
         Jacobian: [63x2 double]
         exitflag: 3
    firstorderopt: 0.0883
       iterations: 5
        funcCount: 18
     cgiterations: 0
        algorithm: 'trust-region-reflective'
         stepsize: 4.1539e-04
          message: 'Success, but fitting stopped because change in residuals less than tolerance (TolFun).'

fitinfo 구조체에서 잔차를 가져옵니다.

residuals = fitinfo.residuals;

"이상값"을 기준 모델에서 1.5 표준편차보다 큰 임의의 거리에 있는 점으로 식별하고, 이상값을 제외한 상태에서 데이터를 다시 피팅합니다.

I = abs( residuals) > 1.5 * std( residuals );
outliers = excludedata(xdata,ydata,'indices',I);

fit2 = fit(xdata,ydata,f,'StartPoint',[1 1],...
           'Exclude',outliers);

이상값을 제외하는 효과와 로버스트 피팅에서 이상값에 보다 적은 겹제곱 가중치를 부여하는 효과를 비교합니다.

fit3 = fit(xdata,ydata,f,'StartPoint',[1 1],'Robust','on');

데이터, 이상값 및 피팅의 결과를 플로팅합니다. 정보를 표시하는 범례를 지정합니다.

plot(fit1,'r-',xdata,ydata,'k.',outliers,'m*') 
hold on
plot(fit2,'c--')
plot(fit3,'b:')
xlim([0 2*pi])
legend( 'Data', 'Data excluded from second fit', 'Original fit',...
    'Fit with points excluded', 'Robust fit' )
hold off

이상값을 고려하여 두 피팅에 대한 잔차를 플로팅합니다.

figure 
plot(fit2,xdata,ydata,'co','residuals') 
hold on
plot(fit3,xdata,ydata,'bx','residuals')
hold off