Python入門トップページ

---

線形回帰分析をしてみよう-2 (StatsModels OLS)

問題

次の表はあるスーパーマーケットの売り場面積と売上高の実績値です．ここでは，線形回帰を StatsModels というパッケージを使って行い，新規開店予定の店舗 (P, Q) の売上高を予測しよう．

Python で線形回帰

まずは必要なモジュールをインポートします．

モジュールのインポートコピー
import numpy as np
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt

# 高解像度ディスプレイ用
from IPython.display import set_matplotlib_formats
# from matplotlib_inline.backend_inline import set_matplotlib_formats # バージョンによってはこちらを有効に
set_matplotlib_formats('retina')

データは GitHub のリポジトリで公開しているので，これを Pandas のデータフレームに直接読み込みます．

問題のデータを読み込むコピー
url = "https://github.com/rinsaka/sample-data-sets/blob/master/store-space-sales.csv?raw=true"
df = pd.read_csv(url)
df

Pandas のデータフレーム df から「space」列と「sales」列を取り出してそれぞれ NumPy 配列に変換します．

Pandas データフレームを NumPy 配列に変換コピー
x = df.loc[:, 'space'].values
y = df.loc[:, 'sales'].values
print(x)
print(y)

[26 28 25 26 28 33 32 29 30 29 38 40 32 26 34]
[65 65 62 59 69 73 79 71 73 71 86 88 71 63 80]

回帰分析の前に散布図を描いてみよう．

散布図コピー
fig, ax = plt.subplots(1, 1, figsize=(6, 4))
ax.scatter(x, y)
ax.set_xlabel('space')
ax.set_ylabel('sales')
ax.set_xlim(20,45)
ax.set_ylim(40,100)
plt.show()

y 切片も含めて回帰分析をしたいので，定数項をxに加えます．加えなければ当てはめた直線は原点 (0, 0) を通ることになります．

定数項を加えるコピー
x1 = sm.add_constant(x)

線形回帰モデルを定義します．

モデルの定義コピー
model = sm.OLS(y, x1)

当てはめ（パラメータ推定）を行う

当てはめコピー
results = model.fit()

推定したパラメータ（y切片，傾き）を表示します．

結果の表示コピー
print(results.params)

[15.92603266  1.83357349]

詳細な結果を表示します．

結果の詳細表示コピー
print(results.summary2())

                 Results: Ordinary least squares
=================================================================
Model:              OLS              Adj. R-squared:     0.903
Dependent Variable: y                AIC:                73.7815
Date:               2021-12-15 22:53 BIC:                75.1976
No. Observations:   15               Log-Likelihood:     -34.891
Df Model:           1                F-statistic:        131.8
Df Residuals:       13               Prob (F-statistic): 3.54e-08
R-squared:          0.910            Scale:              7.0803
-------------------------------------------------------------------
            Coef.    Std.Err.      t      P>|t|    [0.025    0.975]
-------------------------------------------------------------------
const      15.9260     4.9034    3.2480   0.0064   5.3329   26.5192
x1          1.8336     0.1597   11.4810   0.0000   1.4886    2.1786
-----------------------------------------------------------------
Omnibus:               0.733        Durbin-Watson:          2.200
Prob(Omnibus):         0.693        Jarque-Bera (JB):       0.721
Skew:                  -0.339       Prob(JB):               0.697
Kurtosis:              2.168        Condition No.:          219
=================================================================

グラフを表示しよう．

グラフを描画コピー
fig, ax = plt.subplots(1, 1, figsize=(6, 4))
ax.scatter(x, y)
ax.plot(x, results.predict(x1))
ax.set_xlabel('space')
ax.set_ylabel('sales')
ax.set_xlim(20,45)
ax.set_ylim(40,100)
ax.set_title('OSL results')
# plt.savefig('ols-2-2021.png', dpi=300, facecolor='white')
plt.show()

新規店舗の売上を予測します．

予測コピー
print(f'P (27) ==> {results.params[0]  + results.params[1] * 27:.2f}')
print(f'Q (38) ==> {results.params[0]  + results.params[1] * 38:.2f}')

P (27) ==> 65.43
Q (38) ==> 85.60

パラメータを取得するとき，このような書き方もできます．

予測コピー
b, a = results.params  # b:切片， a:傾き
# 予測する
print(f'P (27) ==> {b + a * 27:.2f}')
print(f'Q (38) ==> {b + a * 38:.2f}')

P (27) ==> 65.43
Q (38) ==> 85.60

---