linear regression_python

선형회귀는 회귀계수 w=(w1, w2, ..., wp)를 가진 선형모형을 생성하는 것입니다. 이는 반응변수의 관찰치와 선형 모형에 의한 추정치 사이의 잔차제곱합이 최소가 되도록 설정합니다.
$$ \text{min} \sum^n_{i=1} {(x_i w- y_i)^2}$$

회구모형을 구축하기 위해 sklearn 모듈을 사용하는데 제공되는 모든 함수는 list 또는 np.array 형을 대상으로 합니다. 그러므로 아래와 같이 pandas 모듈의 함수를 사용하여 엑셀로 부터 데이터를 호출하는 경우 호출되 데이터는 DataFrame 형태이므로 이들의 형식을 np.array 형으로 변경하여야 합니다.

>>> import numpy as np
>>> import pandas as pd
>>> from datetime import datetime, timedelta
>>> import matplotlib.pyplot as plt
>>> from sklearn import linear_model

>>> kl=pd.read_excel('C:\\~~.xlsx', sheet_name="kl")
>>> startd=pd.Timestamp(datetime(2017,12, 1))
>>> kl1=kl.ix[kl.index>=startd, :]

>>>  kl1
             Open   High    Low  Close    Volume
2017-12-01  17550  17580  17345  17425   9696255
2017-12-04  17550  17750  17400  17750   8287180
2017-12-05  17610  17955  17530  17875   9387273

....

위 데이터에서 Open, High, Low를 독립변수로  Close를 반응변수로 합니다. 이 경우 어제의 독립변수를 당일의 반응변수와 한행에 놓이도록 재배치 합니다. 

>>> ind=np.array(kl1.iloc[:, :3].shift(1))[1:]

>>> de=np.array(kl1.iloc[1:,3])

>>> new=np.array(kl1.iloc[(n-1), :3])

위의 new는 모형에서 최종적으로 독립변수로 사용될 데이터들로서 kl1의 마지막 행에서 추출한 변수들입니다. 

이들에 대한 선형회귀 모형을 구축하기 위해 다음의 절차를 갖습니다. 

sklearn에서 LinearRegression은 클래스 입니다. 그러므로 선형모형에서 적용할 수 있는 함수들은 메소드 형태로 제공되기 때문에 모형을 위한 클래스 객체를 먼저 생성하여야 합니다. 

>>> m=linear_model.LinearRegression()

>>> m

LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normalize=False)

위에서 생성된 클래스 객체의 메소드를 적용합니다. 먼저 독립변수와 반응변수를 전달하여 그에 적합한 즉 최소자승법(관찰치와 추정치의 차이의 제곱의 합이 최소가 되게 선형회귀선을 생성)에 의해 모형을 구축합니다. 

>>> m.fit(ind, de)

LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normalize=False)

이렇게 생성된 모형의 절편과 회귀계수는 다음과 같이 호출할 수 있습니다. 이 예의 경우 독립변수는 3개의 변수들로 구성됩니다. 그러므로 각각에 대응되는 회귀계수 역시 3개가 됩니다. 

>>> m.intercept_ #절편

11216.163169900159

>>> m.coef_ #회귀계수 

array([-0.10222574, -0.1565665 ,  0.62119229])

이 모형에 대한 추정값은 predict() 메소드를 사용합니다. 다음은 new에 대한 추정입니다.

>>> m.predict(new)

array([ 17264.47689106])

이 모형에 대한 tss, sse, mse, sigma, r^2를 계산합니다. 

tss: total sum error로 모든 독립변수에 대한 추정값이 관찰치라고 가정하는 경우의 각 관찰치와 추정치의 차이에 대한 제곱의 합

$$text{TSS} =\\sum^n_{i=1} (y_i - \bar y)$$

sse: sum square error 모형에 의한 추정치과 관찰치의 차이에 대한 제곱합 

$$text{TSS} =\\sum^n_{i=1} (y_i - \hat y)$$

MSE : $\frac{SSE}{n-df}$ 

  n은 자료의 수 df 자유도

sigma : 표준편차의 불편추정치 $\\sigma = sqrt{MSE}$

$R^2 = \frac{TSS-SSE}{TSS}$ 로서 메소드 score()를 사용하여 계산할 수 있습니다. 

>>> tss=np.sum((de-np.mean(de))**2) ; tss

1022423.3333333333

>>> sse=np.sum((de-m.predict(ind))**2); sse

904045.02190034126

>>> mse=sse/(n-1); mse

60269.668126689416

>>> sigma=mse**(0.5); sigma

245.49881491911404

>>> r2_1=(tss-sse)/tss; r2_1

0.11578209101219521

>>> r2=m.score(ind, de); r2

0.11578209101219517

위의 결과에서 결정계수는 매우 낮다. 즉, 독립변수와 반응변수의 관계가 낮다는 것을 의미한다.
위에서 사용한 각 변수에 대한 F 검정을 실시해보면
>>> from sklearn import feature_selection
>>> feature_selection.f_regression(ind, de)
(array([ 0.12716965,  0.36784891,  1.48167167]), ---> F value
 array([ 0.72710626,  0.55462331,  0.24514876])) ---> P value
위 모형은 적합하지 않음을 의미합니다.
그러나 이 분석은 모집단이 정규분포를 따른다는 가정하에 적합합니다.
다음 결과는 분석에 사용된 반응변수의 정규성을 나타내기 위한 것입니다.

이 그림의 왼쪽은 히스토 그램과 정규확률밀도함수 그래프와는 상이한 차이를 보입니다. 그러나 오른쪽 그래프의 경우 정규분포상의 이론 값들과 실제값들 사이에 밀접한 관계가 있습니다.
실제로 위의 그림은 다음 사용자 정의 함수를 적용한 것으로 오른쪽 그래프(QQ plot)의 상관계수와 p vlaue를 반환합니다. 이 결과는 두 변수는 즉, QQ plot의 x, y축의 변수는 밀접한 관계임을 나타냅니다.

위 결과를 보면 R2이 매우 낮습니다. 이는 독립변수와 반응변수 사이의 상관성이 작다는 것을 의미합니다. 위 데이터의 각 변수간의 상관성분석을 실시해 봅니다. 이 분석은 scipy모듈의 peasonr() 함수를 사용하는 것으로 상관계수와 p-value를 반환합니다.

>>> from scipy import stats
>>> [stats.pearsonr(ind[:,i], de) for i in range(3)]
[(0.098425172641939479, 0.72710626404008671),
 (0.16588379948908261, 0.55462330607877863),
 (0.3198649403704249, 0.24514876092784188)]

반응변수에 대응되는 각 독립변수의 상관계수는 매우 낮습니다. 그러므로 위 데이터로 부터의 선형회귀모형에 의해 유의한 결과를 도출하는 것은 어렵습니다.
이 데이터의 샘플수를 증가시켜봅니다. 다음 데이터는 총 1900여개로 구성됩니다. 이전의 데이터에 비해 매우 증가한 수로서 다음의 결과를 보입니다.
>>> len(klarr)
Out[25]: 1939

>>> m1=linear_model.LinearRegression()
>>> m1.fit(klarr[:,:3], klarr[:,3])
LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normalize=False)
>>> m1.intercept_
13.44051965610015

>>> m1.coef_
array([-0.57760303,  0.82535053,  0.7506489 ])

>>> m1.predict(new)
array([ 16932.76236813])

>>> m1.score(klarr[:,:3], klarr[:,3])
0.99873913782587753

>>> [stats.pearsonr(klarr[:,i], klarr[:,3]) for i in range(3)]
[(0.9965192878679775, 0.0),
 (0.99832737246267556, 0.0),
 (0.99852703764332618, 0.0)]

>>> feature_selection.f_regression(klarr[:,:3], klarr[:,3])
(array([ 276795.75255055,  577576.74481848,  656066.10736815]),
 array([ 0.,  0.,  0.]))

다음 분석 당일을 기준으로 75개의 행만으로 분석한 것입니다.
>>> startd=pd.Timestamp(datetime(2017,9, 1))
>>> klarr=np.array(kl.iloc[kl.index>=startd,:])
>>> len(klarr)
75
>>> m1=linear_model.LinearRegression()
>>> m1.fit(klarr[:,:3], klarr[:,3])
LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normalize=False)

>>> m1.intercept_
-27.417147115982516

>>> m1.coef_
array([-0.62246893,  0.62869076,  0.99812866])

>>> m1.predict(new)
array([ 16941.88535165])

>>> m1.score(klarr[:,:3], klarr[:,3])
0.99314889702931741

>>> [stats.pearsonr(klarr[:,i], klarr[:,3]) for i in range(3)]
[(0.97888753904626236, 4.3938524460377899e-52),
 (0.98952175356209637, 4.1505370086800028e-63),
 (0.99336371499953235, 2.5513971041206362e-70)]

>>> feature_selection.f_regression(klarr[:,:3], klarr[:,3])
(array([ 1674.28153332,  3428.75340488,  5445.37592932]),
 array([  4.39385245e-52,   4.15053701e-63,   2.55139710e-70])