Python入門トップページ

---

線形回帰分析をしてみよう-3 (sklearn)

問題

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

Python で線形回帰

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

モジュールのインポートコピー
import numpy as np
import pandas as pd
from sklearn import linear_model
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 配列に変換します．このとき，x に該当する「space」は1次元配列ではなく，2次元配列（つまり行列形式）になるように変換する必要があることに注意します．

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]

NumPy 配列の形式を確認します．

x の形式を確認するコピー
x.shape

(15, 1)

y の形式を確認するコピー
y.shape

(15,)

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

散布図コピー
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()

データの準備ができたので，モデルを準備します．

モデルの準備コピー
model = linear_model.LinearRegression()

回帰分析は2次元のNumPy配列に変換した x, y を引数に与えると良いでしょう．

回帰分析の実行コピー
results = model.fit(x, y)

結果を表示してみよう．

結果の表示コピー
results

LinearRegression()

具体的に results にはどのような結果が格納されているか，すべての中身を表示します．

結果の表示コピー
vars(results)

{'fit_intercept': True,
 'copy_X': True,
 'n_jobs': None,
 'positive': False,
 'n_features_in_': 1,
 'coef_': array([1.83357349]),
 'rank_': 1,
 'singular_': array([16.66133248]),
 'intercept_': 15.926032660902976}

係数を表示してみよう．

係数の表示コピー
model.coef_

array([1.83357349])

上の係数は array 形式になっているので，次のような方法で値を取得して表示してみよう．

係数の取得と表示コピー
a = model.coef_[0]
print(a)

1.8335734870317013

次は切片を表示してみよう．

切片の表示コピー
model.intercept_

15.926032660902955

切片も値を取得して表示してみよう．

切片の取得と表示コピー
b = model.intercept_
print(b)

15.926032660902955

決定係数 (R-squared) は score 関数で取得できます（もちろん，ここでの R-squared（係数，傾きなど）の結果と一致するはずです）．

決定係数の表示コピー
model.score(x, y)

0.9102297512020177

グラフを表示しよう．

グラフを描画コピー
fig, ax = plt.subplots(1, 1, figsize=(6, 4))
ax.scatter(x, y)
ax.plot(x, [a * t + b for t in x])
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('P (27) ==> ', a * 27 + b)
print('Q (38) ==> ', a * 38 + b)

P (27) ==>  65.4325168107589
Q (38) ==>  85.60182516810758

なお，sklearn では predict を使って予測ができます．

予測(1)コピー
model.predict([[27]])[0]

65.4325168107589

予測(2)コピー
model.predict([[38]])[0]

85.60182516810761

さらに，x が 27 と 38 という2つの新規店舗の売上高がそれぞれ 61 と 79 であることがわかっている状況で，これらのデータをテストデータとして，テストデータに対する当てはまりの良さを検証します．まず，NumPy の行列形式でテストデータを準備します．

テストデータの準備コピー
x_test = np.array([[27], [38]])
y_test = np.array([61, 79])
print(x_test)
print(y_test)

[[27]
 [38]]
[61 79]

テストデータを NumPy 配列で準備することで，予測もまとめて行うことができます．

コピー
y_pred = model.predict(x_test)
print(y_pred)

[65.43251681 85.60182517]

当てはまりの良さは決定係数で検証できます．これは model.score() 関数にテストデータを与えるだけで良いでしょう．

スコア：決定係数コピー
model.score(x_test, y_test)

0.6096833282227203

なお，この決定係数の計算方法はここでは次のように記載されています． \begin{eqnarray*} R^2 = 1 - \frac{u}{v} \end{eqnarray*} ここで \(u\) は残差の二乗和，\(v\) は標本値と標本平均値の偏差の二乗和です．これらを順に計算して決定係数が正しいことも確認しておこう．

まず，NumPy 配列形式のデータでは要素ごとの演算を簡単に行えるので，残差は次のように求めることができます．

残差コピー
y_test - y_pred

array([-4.43251681, -6.60182517])

残差を二乗します．

残差の二乗コピー
(y_test - y_pred)**2

array([19.64720528, 43.58409555])

これらの総和を求めるためには NumPy の sum() 関数が利用できます．これで，残差の二乗和 \(u\) を計算できます．

残差二乗和コピー
u = ((y_test - y_pred)**2).sum()
print(u)

63.2313008279193

次に，標本平均値は次のように y_test.mean() で求めることができます．

コピー
print(y_test)
print(y_test.mean())

[61 79]
70.0

平均値との偏差を求める．

コピー
print(y_test)
print(y_test.mean())
print(y_test - y_test.mean())

[61 79]
70.0
[-9.  9.]

偏差の二乗を計算します．

コピー
print((y_test - y_test.mean())**2)

[81. 81.]

偏差の二乗和 \(u\) は次のように計算できます．

偏差の二乗和コピー
v = ((y_test - y_test.mean())**2).sum()
print(v)

162.0

上の手順をまとめると，決定係数は次のように計算できます．

決定係数コピー
u = ((y_test - y_pred)**2).sum()
v = ((y_test - y_test.mean())**2).sum()
r2 = 1 - u/v
print(f'{r2:.4f}')

0.6097

したがって，score()で計算した値と一致することが確認できました．

スコア（決定係数）の計算：再掲コピー
print(f'{model.score(x_test, y_test):.4f}')

0.6097

なお，回帰直線は x と y それぞれの平均値の交点を通ることに注意してください．それぞれの平均値を表示したあと，回帰式に x の平均値を代入するとその結果は y の平均値と等しくなります．グラフでもそのことを確認しましょう．

コピー
print(f'面積の平均 : {np.mean(x):.2f}')
print(f'売上の平均 : {np.mean(y):.2f}')
print(f'平均面積の店舗の売上 ({np.mean(x)}) ==> {a * np.mean(x) + b:.2f}')

面積の平均 : 30.40
売上の平均 : 71.67
平均面積の店舗の売上 (30.4) ==> 71.67

---