전산물리학:선형_회귀_분석_linear_regression_analysis

선형 회귀 분석(Linear regression analysis)은 $x$와 $y$의 값을 알고 있을 때, $y = ax +b$의 꼴에서 $a$와 $b$의 값을 추측하는데 사용하는 방법이다. 데이터 포인트 하나만 가지고는 $a$와 $b$값을 추론할 수 없기 때문에 $x$에 대한 $y$ 데이터가 많이 있을 때 함수 형태가 어떻게 되는지 추론할 때 쓰인다. 데이터의 수가 많기 때문에 아래와 같이 행렬 꼴로 많이 쓰인다. $$ \mathbf{Y} = a\mathbf{X} + b $$

이 때, $\mathbf{Y} = (y_1, y_2, y_3, \ldots)^\intercal$, $\mathbf{X} = (x_1, x_2, x_3, \ldots)^\intercal$ 이다.

그 중 널리 쓰이는 방법이 최소제곱법(Least-square fitting)이다. 최소제곱법은 데이터로 얻은 $\mathbf{Y}$와 우리가 추론한 함수의 결과인 $f(\mathbf{X})$ 사이의 오차의 제곱을 최소화 하는 함수 $f(x)$를 찾는 방법이다. 즉, 다음의 식에서 $\chi^2$를 최소화하는 함수 $f(x)$를 찾는 것이다. $$ \chi^2 = \sum_{i=1}^{N}\left(\frac{y_i - f(x_i)}{\sigma_i}\right)^2$$

Python의 경우 scipy의 curve_fitting 함수를 사용할 수 있고, C/C++의 경우 gsl(gnu scientific library)의 gsl_fit_linear 혹은 gsl_fit_wlinear를 사용할 수 있다. 데이터에 표준오차가 포함되어 있는 경우에는 wlinear를, 포함되지 않은 경우에는 linear를 사용한다. wlinear를 쓸 때는 표준오차 값의 제곱의 역수를 w값으로 넣어줘야 한다. gsl의 경우 $$\chi^2 = w(y-f(x))^2;\quad w = \frac{1}{\sigma^2}$$ 으로 계산하기 때문이다.

scipy를 쓸 때는 표준오차가 포함되어 있는 경우에는 absolute_sigma 옵션을 True로 줘야 한다.

실습

다음의 데이터를 예로 각 언어로 적합(fitting)해보자.

x y $\sigma$
0.0 1.5891700001638700 0.5320185951347170
0.25 1.2046979031630900 0.45427229254919400
0.5 2.598484110379440 0.5077027054798450
0.75 2.2043375760987900 0.5205826529502230
1.0 3.002229375836470 0.4698090780151540
1.25 2.928655222350250 0.4773624121076510
1.5 3.4769844321649000 0.5018787752525470
1.75 4.007335261529810 0.4897362325985810
2.0 4.534192485682100 0.49462141854510300
2.25 5.76185759925369 0.5405394898833090
2.5 6.41954394976063 0.4523449367708450
2.75 6.023919923933230 0.48898047501417500
3.0 7.443742735575060 0.5342659697011440
3.25 7.870686775946450 0.4730380798563350
3.5 8.022587710435350 0.4572104512757930
3.75 9.395420889477870 0.49902577926821200
4.0 9.948288024421440 0.530527413902035
4.25 9.986332770533490 0.48969676619036700
4.5 9.694674681091840 0.5443240285483060
4.75 10.519384785632700 0.47616567472190800

먼저 Python 코드는 아래와 같다.

| curve_fit.py
import numpy as np
from scipy.optimize import curve_fit
 
func = lambda x, a, b: a*x + b #f(x) = ax + b
 
x = np.linspace(0, 4.75, 20) #x = [0.00, 0.25, 0.50, ...]
 
y = np.array([1.5891700001638700, 1.2046979031630900, 2.598484110379440, 2.2043375760987900, 3.002229375836470, 2.928655222350250, 3.4769844321649000, 4.007335261529810, 4.534192485682100, 5.76185759925369, 6.41954394976063, 6.023919923933230, 7.443742735575060, 7.870686775946450, 8.022587710435350, 9.395420889477870, 9.948288024421440, 9.986332770533490, 9.694674681091840, 10.519384785632700])
 
sig = np.array([0.5320185951347170, 0.45427229254919400, 0.5077027054798450, 0.5205826529502230, 0.4698090780151540, 0.4773624121076510, 0.5018787752525470, 0.4897362325985810, 0.49462141854510300, 0.5405394898833090, 0.4523449367708450, 0.48898047501417500, 0.5342659697011440, 0.4730380798563350, 0.4572104512757930, 0.49902577926821200, 0.530527413902035, 0.48969676619036700, 0.5443240285483060, 0.47616567472190800])
 
#if saved file is existed
"""
x, y, err = np.loadtxt('line.data', unpack = True)
"""
 
NDF = x.size - 2 # # of Degree of Freedom
popt, pcov = curve_fit(func, x, y, sigma = sig, absolute_sigma = True) #popt: inferenced coefficient, pcov: covariance between coefficients
perr = np.sqrt(np.diag(pcov))# diagonal elements of covariance matrix are error of each coefficient
chi2 = np.sum((y - func(x,*popt))**2 / sig**2)
 
print("y = {0}x + {1}".format(popt[0],popt[1])

다음은 C++에서 gsl을 이용한 코드이다.

| curve_fit.cpp
#include <iostream>
#include <gsl/gsl_fit.h>
#include <cmath>
 
using std::cout;
using std::pow;
using std::sqrt;
 
// #include <fstream>
// #include <sstream>
// #include <string>
// #include <vector>
 
//using std::ifstream;
//using std::istringstream;
//using std::vector;
//using std::string;
//using std::stod;
//using std::copy;
 
int main()
{
    double x[] = {0.0, 0.25, 0.5, 0.75, 1.0, 1.25, 1.5, 1.75, 2.0, 2.25, 2.5, 2.75, 3.0, 3.25, 3.5, 3.75, 4.0, 4.25, 4.5, 4.75};
    double y[] = {1.5891700001638700, 1.2046979031630900, 2.598484110379440, 2.2043375760987900, 3.002229375836470, 2.928655222350250, 3.4769844321649000, 4.007335261529810, 4.534192485682100, 5.76185759925369, 6.41954394976063, 6.023919923933230, 7.443742735575060, 7.870686775946450, 8.022587710435350, 9.395420889477870, 9.948288024421440, 9.986332770533490, 9.694674681091840, 10.519384785632700};
    double sig[] = {0.5320185951347170, 0.45427229254919400, 0.5077027054798450, 0.5205826529502230, 0.4698090780151540, 0.4773624121076510, 0.5018787752525470, 0.4897362325985810, 0.49462141854510300, 0.5405394898833090, 0.4523449367708450, 0.48898047501417500, 0.5342659697011440, 0.4730380798563350, 0.4572104512757930, 0.49902577926821200, 0.530527413902035, 0.48969676619036700, 0.5443240285483060, 0.47616567472190800};
 
    //if saved file is existed
    /*
    string tmp;
    std::ifstream fin;
    fin.open("line.dat");
 
    vector<double> xv;
    vector<double> yv;
    vector<double> sigv;
    vector<string> sv;
    string sbf;
 
    while(!fin.eof())
    {
        getline(fin, tmp);
        if(tmp == "") {break;}
        istringstream ss(tmp);
        while(getline(ss,sbf,','))
        {
            sv.emplace_bacK();
        }
        sv.shrink_to_fit();
 
        xv.emplace_back(stod(sv[0]));
        yv.emplace_back(stod(sv[1]));
        sigv.emplace_back(stod(sv[2]));
 
        sv.clear();
        sv.shrink_to_fit();
    }
    xv.shrink_to_fit();
    yv.shrink_to_fit();
    sigv.shrink_to_fit();
 
    double x[xv.capacity()];
    double y[yv.capacity()];
    double sig[sig.capacity()];
 
    copy(xv.begin(), xv.end(), x);
    copy(yv.begin(), yv.end(), y);
    copy(sigv.begin(), sigv.end(), sig);
    */
 
    for(int i = 0; i < sizeof(sig)/sizeof(double); ++i)
    {
        sig[i] = 1/sig[i];
    }
 
    double a, b, cov00, cov01, cov11, chisq;
    gsl_fit_wlinear(x, 1, sig, 1, y, 1, sizeof(x)/sizeof(double), &b, &a, &cov00, &cov01, &cov11, &chisq);
    cout << "y = " << a << "x + "<< b << "\n";
}

각각의 결과는 다음과 같다 $$\text{Python:}\quad y = 2.09749(\pm 0.0771828)x + 0.838641(\pm 0.213449)$$ $$\text{C++:}\quad y = 2.09749(\pm 0.0771828)x + 0.838641(\pm 0.0.213449)$$ 유효숫자는 c++ cout의 기본 자릿수인 6자리에 맞췄다.

  • 전산물리학/선형_회귀_분석_linear_regression_analysis.txt
  • Last modified: 2022/01/17 16:12
  • by jonghoon