선형 최소 제곱법

참고

번역중

이 단원에서는 선형 최소 제곱법을 이용해 실험 값을 여러 함수들의 조합으로 표현하는 기능을 서술합니다. 각 계산은 알려진 오차를 포함하거나 포함하지 않을 수도 있습니다. 가중치가 있는 값들에 대해, 각 함수들은 최적의 계수를 계산하고 그들에 관한 공분산 행렬를 찾습니다. 가중치가 없는 값들은 각 지점의 분포를 이용해 공분산 행렬을 찾고, 분산-공분산 행렬을 제공합니다.

이 단원의 함수들은 여러 형태로 나누어져 있습니다. 단순한 1개, 2개의 계수를 받는 분석 함수와 여러 계수를 사용 가능한 분석 함수로 나누어져 있습니다.

개요

최소 제곱법은 \(\chi^2\) 을 최소화 시켜 찾을 수 있습니다. 이 값은 모델 \(Y(c,x)\) 에 대해, \(n\) 개의 실험 값들의 가중 잔차 제곱합을 의미합니다.

\[\chi^2 = \sum_i w_i (y_i - Y(c,x_i))^2\]

모델의 \(p\) 계수는 \(c=\{ c_0, c_1, \dots \}\) 입니다. 가중 계수 \(w_i\) 는 \(w_i = 1/\sigma_i^2\) 로 주어집니다. \(\sigma_i\) 는 지점 \(y_i\) 의 실험 오차입니다. 이 오차들은 정규 분포(가우스 분포)를 따르고 비상관(uncorrelated) 관계에 있다고 가정합니다. 가중치가 없는 값들에 대해, \(\chi^2\) 은 가중치 계수 없이 계산됩니다.

피팅 함수들은 가장 잘 맞는 계수 \(c\) 와 \(p\times p\) 차원의 공분산 행렬을 반환합니다. 공분산 행렬은 값들에 있는 오차( \(\sigma_i\) )로 인해 생기는 계수 \(c\) 의 통계적 오차를 측정합니다. 이 행렬은 다음과 같이 정의됩니다.

\[C_{ab} = \braket{\delta{c_a} \delta{c_b}}\]

\(\braket{}\) 은 해당 값 아래, 가우스 오차 분포 함수의 평균값을 나타냅니다.

공분산 행렬은 각 지점의 오차 \(\sigma_i\) 의 오차 전파로 계산할 수 있습니다. 값 \(\sigma_{y_i}\) 의 변화로 생기는 피팅 계수의 변화량 \(\delta_{c_a}\) 은 다음과 같이 주어집니다.

\[\sigma{c_a} = \sum_i \frac{\partial c_a}{\partial y_i} \delta y_i\]

비상관 값들은, 다음을 만족합니다.

\[\braket{\delta y_i \delta y_j} = \sigma_i^2 \delta_{ij}\]

이때, 해당하는 공분산 행렬은 다음과 같이 계산할 수 있습니다.

\[C_{ab} = \sum_i \frac{1}{w_i} \frac{\partial c_a}{\partial y_i} \frac{\partial c_b}{\partial y_i}\]

가중치가 없는 값들은 공분산 행렬을 계산할때, 급수에서 가중 계수 \(w_i\) 가 단일값 \(w = \frac{1}{\sigma^2}\) 로 대체됩니다. \(\sigma^2\) 는 최적 추정 계수에 대한, 잔차 분포로 \(\sigma^2 = \sum (y_i - Y(c,x_i))^2/(n-p)\) 와 같이 계산됩니다. 이를 분산-공분산 행렬이라 합니다.

최적 추정 계수들에 대한 표준 편차는 해당하는 공분산 함수의 대각 성분의 제곱근으로 계산됩니다. \(\sigma_{c_a} = \sqrt{C_{aa}}\) . 계수 \(c_a\) 와 \(c_b\) 의 상관 계수는 \(\rho_{ab} = C_{ab}/\sqrt{C_{aa}C_{bb}}\) 로 주어집니다.

선형 회귀(Linear regression)

이 단원의 함수들은 계수가 1, 2개인 간단한 선형 회귀 모델을 계산합니다. 이 단원의 함수들은 헤더 파일 gsl_fit.h 에 정의되어 있습니다.

상수항을 가지는 선형 회귀 모델

이 단원의 함수들은 일직선 형태의 모델 \(Y(c,x) = c_o +c_1 x\) 에 대한 선형 회귀 분석을 계산합니다.

상수항이 없는 선형 회귀 모델

여러 계수를 가지는 선형 모델

정규화된 회귀 분석(Regularized regression)

로버스트 회귀 분석(Robust regression)

대규모 선형계(Large dense liear systems)

1

예제

다음의 프로그램은 선형 최소 제곱법을 사용해 주어진 데이터에 맞는 직선 식을 찾습니다. 그리고 최적의 피팅선과 버금 표준 편차를 출력합니다.

#include <stdio.h>
#include <gsl/gsl_fit.h>

int
main (void)
{
  int i, n = 4;
  double x[4] = { 1970, 1980, 1990, 2000 };
  double y[4] = {   12,   11,   14,   13 };
  double w[4] = {  0.1,  0.2,  0.3,  0.4 };

  double c0, c1, cov00, cov01, cov11, chisq;

  gsl_fit_wlinear (x, 1, w, 1, y, 1, n,
                   &c0, &c1, &cov00, &cov01, &cov11,
                   &chisq);

  printf ("# best fit: Y = %g + %g X\n", c0, c1);
  printf ("# covariance matrix:\n");
  printf ("# [ %g, %g\n#   %g, %g]\n",
          cov00, cov01, cov01, cov11);
  printf ("# chisq = %g\n", chisq);

  for (i = 0; i < n; i++)
    printf ("data: %g %g %g\n",
                   x[i], y[i], 1/sqrt(w[i]));

  printf ("\n");

  for (i = -30; i < 130; i++)
    {
      double xf = x[0] + (i/100.0) * (x[n-1] - x[0]);
      double yf, yf_err;

      gsl_fit_linear_est (xf,
                          c0, c1,
                          cov00, cov01, cov11,
                          &yf, &yf_err);

      printf ("fit: %g %g\n", xf, yf);
      printf ("hi : %g %g\n", xf, yf + yf_err);
      printf ("lo : %g %g\n", xf, yf - yf_err);
    }
  return 0;
}

다음의 명령어들은 프로그램의 출력값으로부터 데이터를 뽑아내고 GNU plotutils “graph” 도구를 이용해 시각 그래프를 만들어줍니다.

$./demo > tmp
$more tmp
# best fit: Y = -106.6 + 0.06 X
# covariance matrix:
# [ 39602, -19.9
#   -19.9, 0.01]
# chisq = 0.8

$for n in data fit hi lo ;
   do
     grep "^ :math:`n" tmp | cut -d: -f2 > ` n ;
   done
$graph -T X -X x -Y y -y 0 20 -m 0 -S 2 -Ie data
     -S 0 -I a -m 1 fit -m 2 hi -m 2 lo

결과는 다음과 같습니다.

참고 문헌과 추가 자료

최소 제곱법과 관련된 수식과 기법들은 Particle Data Group에서 출판한 The Review of Particle Physics의 “Statistics” 단원을 참고할 수 있습니다.

  • Review of Particle Properties, R.M. Barnett et al., Physical Review D54, 1 (1996) http://pdg.lbl.gov

The Review of Particle Physics은 위의 링크에서 볼 수 있습니다.

이 단원에서 구현된 기능들을 검사하는 데 NIST Statistical Reference Datasets을 사용했습니다. 해당 값과 문서들은 NIST 사이트를 참고할 수 있습니다.

http://www.nist.gov/itl/div898/strd/index.html

Tikhonov regularization에 대한 자세한 정보는 다음을 참고할 수 있습니다.

  • Hansen, P. C. (1998), Rank-Deficient and Discrete Ill-Posed Problems: Numerical Aspects of Linear Inversion. SIAM Monogr. on Mathematical Modeling and Computation, Society for Industrial and Applied Mathematics

  • M. Rezghi and S. M. Hosseini (2009), A new variant of L-curve for Tikhonov regularization, Journal of Computational and Applied Mathematics, Volume 231, Issue 2, pages 914-924.

GSL의 로버스트 선형 회귀 구현체는 다음 출판물에 기반해 있습니다.

  • DuMouchel, W. and F. O’Brien (1989), “Integrating a robust option into a multiple regression computing environment,” Computer Science and Statistics: Proceedings of the 21st Symposium on the Interface, American Statistical Association

  • Street, J.O., R.J. Carroll, and D. Ruppert (1988), “A note on computing robust regression estimates via iteratively reweighted least squares”, The American Statistician, v. 42, pp. 152-154.

정규 방정식들과 TSQR을 이용한 대규모 선형 최소 제곱계의 풀이는 다음을 참고할 수 있습니다.

  • Trefethen, L. N. and Bau, D. (1997), “Numerical Linear Algebra”, SIAM.

  • Demmel, J., Grigori, L., Hoemmen, M. F., and Langou, J. “Communication-optimal parallel and sequential QR and LU factorizations”, UCB Technical Report No. UCB/EECS-2008-89, 2008.