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

Differences

This shows you the differences between two versions of the page.

Link to this comparison view

Next revision		 Previous revision
전산물리학:선형_회귀_분석_linear_regression_analysis [2022/01/15 00:09] – created jonghoon전산물리학:선형_회귀_분석_linear_regression_analysis [2023/09/05 15:46] (current) – external edit 127.0.0.1
Line 1: Line 1:
 선형 회귀 분석(Linear regression analysis)은 $x$와 $y$의 값을 알고 있을 때, $y = ax +b$의 꼴에서 $a$와 $b$의 값을 추측하는데 사용하는 방법이다. 선형 회귀 분석(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 =====
 +먼저 Python 코드는 아래와 같다.
 +
 +<code: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])
 +</code>
 +
 +===== C++ =====
 +다음은 C++에서 gsl을 이용한 코드이다.
 +<code:C++ | 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";
 +}
 +</code>
 +
 +===== 결과 =====
 +각각의 결과는 다음과 같다
 +$$\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.1642172945.txt.gz
  • Last modified: 2023/09/05 15:46
  • (external edit)