2018/01/13
前回の続きです。今回は、Pythonで回帰分析を行います。
３．Python Pythonにおける「回帰分析」では、「scikit-learnを用いる方法」と「StatsModelsを用いる方法」の2パターンを試してみました。
まず、scikit-learnを用いる方法です。プログラムは以下の通りです。
import numpy as np
import matplotlib.pyplot as plt
from sklearn import linear_model
x = np.array([107,336,233,82,61,378,129,313,142,428,100,350,230,70,50,380,120,300,135,430]).reshape(-1, 1)
y = np.array([286,851,589,389,158,1037,463,563,372,1020,170,845,570,300,160,1050,450,550,380,1030]).reshape(-1, 1)
# Linear Regression (y=ax+b)
lrm = linear_model.LinearRegression()
lrm.fit_intercept = True
lrm.fit(x, y)
print('< y = ax + b >')
print('Slope = ', (lrm.coef_[0])[0] )
print('Intercept = ', lrm.intercept_[0])
print('R-squared = ', lrm.score(x, y))
plt.scatter(x, y)
plt.xlabel('x')
plt.ylabel('y')
plt.xlim(0, 500)
plt.ylim(0, 1200)
plt.plot(x, lrm.predict(x)) # Regression Line
plt.show()
# Linear Regression (y=ax)
lrm.fit_intercept = False
lrm.fit(x, y)
print('< y = ax >')
print('Slope = ', (lrm.coef_[0])[0] )
print('R-squared = ', lrm.score(x, y))
plt.scatter(x, y)
plt.xlabel('x')
plt.ylabel('y')
plt.xlim(0, 500)
plt.ylim(0, 1200)
plt.plot(x, lrm.predict(x)) # Regression Line
plt.show()
5、6行目でデータを入力しますが、reshape(-1, 1)で一次元配列に変換します。8~22行目は「切片≠0」の場合で、24~36行目は「切片=0」の場合のコードです。違いはfit_interceptパラメータをTrue(切片≠0)にするか、False(切片=0）にするかです。11、26行目でそれぞれ直線回帰を行い、グラフを描画しました。
「切片≠0」の結果は以下の通りです。
< y = ax + b >
Slope = 2.19202803591
Intercept = 82.2534685472
R-squared = 0.901455294047
結果は、
前々回のExcel、前回のRと同じ結果でした。
「切片=0」の結果は以下の通りです。
< y = ax >
Slope = 2.46992382474
R-squared = 0.881847019096
この決定係数は、Rの結果と異なり、怪しい値と言っていたExcelグラフの近似曲線から求めたものと同じ値になりました。いったい何が正しいのでしょうか？謎解きの前に、StatsModelsを用いる方法でも回帰分析を行いました。プログラムは以下の通りです。
import numpy as np
import matplotlib.pyplot as plt
from statsmodels import api as sm
import statsmodels.formula.api as smf
x = np.array([107,336,233,82,61,378,129,313,142,428,100,350,230,70,50,380,120,300,135,430])
y = np.array([286,851,589,389,158,1037,463,563,372,1020,170,845,570,300,160,1050,450,550,380,1030])
model = smf.OLS(y, x) # Intercept = 0
result = model.fit()
print(result.summary())
a = result.params[0]
print('Slope = ', a)
plt.scatter(x, y)
plt.xlabel('x')
plt.ylabel('y')
plt.xlim(0, 500)
plt.ylim(0, 1200)
plt.plot(x, a*x) # Regression Line
plt.show()
XX = sm.add_constant(x)
YY = y
model = smf.OLS(YY, XX)
result2 = model.fit()
print(result2.summary())
a = result2.params[0]
b = result2.params[1]
print('Slope = ', b)
print('Intercept = ', a)
plt.scatter(x, y)
plt.xlabel('x')
plt.ylabel('y')
plt.xlim(0, 500)
plt.ylim(0, 1200)
plt.plot(x, b*x+a) # Regression Line
plt.show()
StatsModelsでは、データをreshape(-1, 1)で一次元配列にする必要はなく処理できました。10行目で直線回帰をすると、結果は「切片=0」のものでした。切片ありにするには22行目のような「定数付き」のおまじないが必要でした。
先に「切片=0」の結果は以下の通りです。
OLS Regression Results
==============================================================================
Dep. Variable: y R-squared: 0.974
Model: OLS Adj. R-squared: 0.972
Method: Least Squares F-statistic: 704.8
Date: Mon, 08 Jan 2018 Prob (F-statistic): 1.75e-16
Time: 17:15:13 Log-Likelihood: -121.11
No. Observations: 20 AIC: 244.2
Df Residuals: 19 BIC: 245.2
Df Model: 1
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
x1 2.4699 0.093 26.549 0.000 2.275 2.665
==============================================================================
Omnibus: 2.296 Durbin-Watson: 2.224
Prob(Omnibus): 0.317 Jarque-Bera (JB): 1.207
Skew: -0.597 Prob(JB): 0.547
Kurtosis: 3.153 Cond. No. 1.00
==============================================================================
Rの結果のように詳しい結果が表示できました。この決定係数(0.974)はまたまた先ほどと(0.8818)異なるものでした。ますます混乱してきました・・。
「切片≠0」の結果は以下の通りです。
OLS Regression Results
==============================================================================
Dep. Variable: y R-squared: 0.901
Model: OLS Adj. R-squared: 0.896
Method: Least Squares F-statistic: 164.7
Date: Mon, 08 Jan 2018 Prob (F-statistic): 1.70e-10
Time: 17:15:13 Log-Likelihood: -119.29
No. Observations: 20 AIC: 242.6
Df Residuals: 18 BIC: 244.6
Df Model: 1
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const 82.2535 43.463 1.893 0.075 -9.058 173.565
x1 2.1920 0.171 12.832 0.000 1.833 2.551
==============================================================================
Omnibus: 1.585 Durbin-Watson: 2.360
Prob(Omnibus): 0.453 Jarque-Bera (JB): 1.022
Skew: -0.547 Prob(JB): 0.600
Kurtosis: 2.832 Cond. No. 498.
==============================================================================
こちらは問題なしです。
Excel、R、Pythonの順に結果を見てきましたが、まとめると以下の通りになります。
y=ax+bの形での回帰分析は三者とも同じ結果で問題なかったのですが、y=axの形の時の決定係数が2パターンに分かれてしまいました。気持ち悪さが残っていますので、時間をかけて調べてみようと思います。
