Differences

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

--- 전산물리학:선형_회귀_분석_linear_regression_analysis [2022/01/16 17:16] – jonghoon
+++ 전산물리학:선형_회귀_분석_linear_regression_analysis [2022/01/17 16:12] – jonghoon
@@ Line 9: / Line 9: @@
 $$ \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를 사용한다. scipy를 쓸 때는 표준오차가 포함되어 있는 경우에는 absolute_sigma 옵션을 True로 줘야 한다.
+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)해보자.
@@ Line 35: / Line 39: @@
 | 4.75 | 10.519384785632700 | 0.47616567472190800 |
+===== Python =====
 먼저 Python 코드는 아래와 같다.
-<Code:Python>
+<code:Python | curve_fit.py>
 import numpy as np
 from scipy.optimize import curve_fit
@@ Line 60: / Line 65: @@
 print("y = {0}x + {1}".format(popt[0],popt[1])
-</Code>
+</code>
+===== C++ =====
 다음은 C++에서 gsl을 이용한 코드이다.
-<Code:C++>
+<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>
-// #include <deque>
-#include <gsl/gsl_fit.h>
 //using std::ifstream;
@@ Line 76: / Line 88: @@
 //using std::string;
 //using std::stod;
-//using std::deque;
+//using std::copy;
 int main()
@@ Line 90: / Line 102: @@
     fin.open("line.dat");
-    vector<double> x;
+    vector<double> xv;
-    vector<double> y;
+    vector<double> yv;
-    vector<double> sig;
+    vector<double> sigv;
     vector<string> sv;
+    string sbf;
     while(!fin.eof())
@@ Line 100: / Line 113: @@
         if(tmp == "") {break;}
         istringstream ss(tmp);
-        string sbf;
         while(getline(ss,sbf,','))
         {
             sv.emplace_bacK();
-            sv.
         }
+        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>
-</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자리에 맞췄다.