【機械学習入門④】線形回帰・重回帰モデルをPythonの実装付きで理解！

前回の記事『【機械学習入門③】最小二乗法とは？手計算とPythonを使った実装を解説【初心者向け】』では、単回帰モデルの導出方法である最小二乗法を学びました。

今回は、もう一度モデルの学習に戻って、単回帰モデルよりも精度を上げやすい重回帰モデルについて触れていきたいと思います。

この記事で学べる内容は、以下のとおりです。

線形回帰・重回帰モデルについて
Pythonを使った重回帰モデルの実装方法

新しいモデルではありますが、単回帰モデルと最小二乗法を分かっていれば、ぶっちゃけ重回帰モデルはカンタンです。

逆に、単回帰モデルと最小二乗法が怪しい人は、以下の記事で復習しましょう！

この記事について

この記事は、機械学習を体系的に学べる「機械学習入門講座」の1つです。

機械学習を勉強しようと思ったとき、以下のような悩みを持ったことはないですか？

機械学習を学んでみたいけど、何から勉強したら良いのか分からない
機械学習の本を買ってみたけど、よく分からなくて挫折した
そもそも本を買うのではなく、まずは「無料」で学習してみたい

上記のような、「機械学習を勉強し始める前の悩み」を解決するために、この機械学習講座を作成しました。

「プログラミングスクールは高いし、参考書は途中で飽きちゃうし…」といった方は、ぜひ第1回から学習して頂けたらと思います(*’ω’*)

一覧ページへ≫

※ブログ記事なので、無料です。でも、作成に手は抜いていません。

Contents

1 線形回帰・重回帰モデルとは？単回帰モデルとの違い
2 Pythonを使った重回帰モデルの実装方法【Scikit-learn】
3 【機械学習入門④】重回帰モデルまとめ

線形回帰・重回帰モデルとは？単回帰モデルとの違い

重回帰モデルは、以下のような式で表されます。

$$y = w_0 + w_1 x_1 + w_2 x_2 + \cdots + w_N x_N$$

要するに、単回帰モデル$y = w_0 + w_1 x_1$の「特徴量が2つ以上ある場合」が、重回帰モデルです。

単回帰モデルでは、1つの特徴量を使ってターゲットの予測をしましたが、重回帰モデルでは複数の特徴量を使って予測をします。

具体例で考える重回帰モデル

以前に考えていたアイスクリームの例で、もう少し具体的に見ていきましょう。

単回帰モデルでは、気温からアイスクリームの売上個数を予測していました。

▼単回帰モデルで使っていた表

気温	アイスクリームの売上
28	10
26	9
20	4
24	7
22	5

ただ、気温だけでなく、1時間当たりの雨量によっても、アイスクリームの売上は左右するはずです。

雨量が0(=つまり雨が降っていない)なら、アイスクリームは売れるだろうし、雨量が多ければ外に出てアイスクリームを買いに行くのも渋ってしまいますよね。

これを変数にして表を作成すると、以下のようになります。

$x_1$ : 気温
$x_2$ : 雨量
$y$ : アイスクリームの売上個数

気温	雨量	アイスクリームの売上個数
28	5	10
26	5	9
20	6	4
24	10	7
22	7	5

数字は適当ですが、だいたいこんな感じの表が作れますね。

重回帰モデルでは、気温と雨量を使って、アイスクリームの売上を予測していくことになります。

そして、重回帰式$y = w_0 + w_1 x_1 + w_2 x_2$ の、$w_0$と$w_1$と$w_2$を求めれば、重回帰モデルを求められます。

単回帰モデルと同じですね。

予測データを使った重回帰モデル

具体的に、Pythonの機械学習ライブラリscikit-learnに入っている、『ボストン住宅価格データセット』を使って考えてみましょう。

ボストン住宅価格データセットの中身

変数名	データの説明
CRIM	犯罪発生率
ZN	住居区画の密集度
INDUS	非小売業の土地割合
CHAS	チャールズ川の周辺か否か
NOX	NOx濃度
RM	住居に含まれる平均部屋数
AGE	1940年より前に建てられた物件の割合
DIS	ボストン市の5つの雇用施設からの重み付き距離
RAD	大きな道路へのアクセスのしやすさ
TAX	1万ドルあたりの不動産税率の総計
PTRATIO	教師あたりの生徒数
B	対象地域に住む黒人の比率
LSTAT	低所得者の割合

上記です。

たとえば、『住居に含まれる平均部屋数=RM』と『1940年より前に建てられた物件の割合=AGE』を使って予測するなら、重回帰式は以下になります。

$${\rm MEDV} = w_0 + w_1 {\rm RM} + w_2 {\rm AGE}$$

※MEDV = 予測対象の住宅価格

この式を使って、単回帰モデルと同様に、各データに対して「なるべく近くなる位置を通る場所」を求めることになります。

重回帰モデルのイメージ

単回帰モデルのときは、以下の図のように、各点の中央を通るような直線を求めました。

単回帰式は$y = w_0 + w_1 x_1$のように書けたので、上の図のように二次元で書けるのは、直感的にも分かります。

では、$y = w_0 + w_1 x_1 + w_2 x_2$のような重回帰式だと、どのように表せば良いのでしょうか？

そもそも$x_1$と$x_2$の2変数、それに加えて予測対象$y$を同時に表すには、どうしたら良いでしょうか？

先に結論を言うと、以下の図のようになります。

単回帰モデルではX-Y軸上で「直線」を求めましたが、重回帰モデルの場合はX-Y-Zの3軸になるので「平面」を求めることになります。

ちょっと意味不明かもしれないですが、「重回帰モデルを考えるときは平面を求めるんだ〜」と捉えて頂ければと。

ちなみに、今回は変数が全部で3つ($x_1, x_2, y$)だから、表示するには3次元でした。

ということは、もっと変数が増えると、それだけ次元が増えるということです。

4次元以上の可視化は難しいので、重回帰モデルを使うときは「平面を求める」ことだけ、おさえておきましょう。

Pythonを使った重回帰モデルの実装方法【Scikit-learn】

重回帰モデルの概要が分かったところで、さっそく実装をしていきましょう。

実装は、単回帰モデルとほとんど同じで、以下の手順で進めていきます。

STEP① : ライブラリの読み込み
STEP② : データの読み込み
STEP③ : データの準備
STEP④ : 重回帰モデルの作成と評価

先にコードを見せてしまうと、以下のようになっています。

# ## 2. 線形回帰・重回帰モデルの実装
### 2.1 ライブラリの読み込み
import numpy as np
import pandas as pd

from sklearn.datasets import load_boston
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split

import matplotlib.pyplot as plt
import seaborn as sns
from mpl_toolkits.mplot3d import Axes3D

sns.set()
sns.set_style("ticks")

# ### 2.2 データの読み込み
boston = load_boston()
df = pd.DataFrame(boston.data, columns=boston.feature_names)
df['MEDV'] = boston.target

print(df.head())

# 散布図を使った可視化
X = df.RM
Y = df.LSTAT
Z = df.MEDV

fig = plt.figure()
ax = Axes3D(fig)

ax.set_xlabel('RM')
ax.set_ylabel('LSTAT')
ax.set_zlabel('MEDV')

ax.scatter(X,Y,Z, s=30, c='blue', alpha=0.2)

plt.show()

# ### 2.3 データの準備
X = df[['RM', 'LSTAT']]
y = df['MEDV']

print('-'*10 + '特徴量とターゲットに分割' + '-'*10)
print('X:', X.shape)
print('y:', y.shape)

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=0)

print('元データ', X.shape)
print('学習用データ', X_train.shape)
print('テスト用データ', X_test.shape)

# ### 2.4 重回帰モデルの作成と評価
lm = LinearRegression()
lm.fit(X_train, y_train)


print('Train Score：{:.2f}' .format(lm.score(X_train, y_train)))
print('Test Score：{:.2f}' .format(lm.score(X_test, y_test)))

print('バイアス', lm.intercept_)
print('重み', lm.coef_)

# ### 2.5 学習結果の可視化

a1, a2 = lm.coef_ #係数
b = lm.intercept_ #切片

X = df.RM
Y = df.LSTAT
Z = df.MEDV

fig = plt.figure()
ax = Axes3D(fig)
ax.scatter(X,Y,Z, s=30, c='blue', alpha=0.2)

X_3d, Y_3d = np.meshgrid(np.arange(0, 10, 1), np.arange(0, 40, 1))
Z_3d = a1 * X_3d + a2 * Y_3d + b
ax.plot_surface(X_3d, Y_3d, Z_3d)
ax.set_xlabel('RM')
ax.set_ylabel('LSTAT')
ax.set_zlabel('MEDV')

ax.view_init(elev=60, azim=45)
plt.show()

だいたい30行くらいで書けますね。

STEP① : ライブラリの読み込み

import numpy as np
import pandas as pd

from sklearn.datasets import load_boston
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split

import matplotlib.pyplot as plt
import seaborn as sns
from mpl_toolkits.mplot3d import Axes3D

sns.set()
sns.set_style("ticks")

まずは、重回帰モデルを実装するためのライブラリを読み込みます。

今回必要になるライブラリをインポートします。

Numpy : 配列計算用ライブラリ
Pandas : データ操作用ライブラリ
Scikit-learn(sklearn) : データの準備・機械学習モデルの作成用ライブラリ
Matplotlib : 作図用ライブラリ①
Seaborn : 作図用ライブラリ②
mpl_toolkits : 3D作図の準備ライブラリ

STEP② : データの読み込み

boston = load_boston()
df = pd.DataFrame(boston.data, columns=boston.feature_names)
df['MEDV'] = boston.target

print(df.head())

# 散布図を使った可視化
X = df.RM
Y = df.LSTAT
Z = df.MEDV

fig = plt.figure()
ax = Axes3D(fig)

ax.set_xlabel('RM')
ax.set_ylabel('LSTAT')
ax.set_zlabel('MEDV')

ax.scatter(X,Y,Z, s=30, c='blue', alpha=0.2)

plt.show()

単回帰モデルのときと同じく、「ボストンの住宅価格」データを使って予測していきます。

データの中身は、以下の表のとおりです。

再掲 : ボストン住宅価格データセットの中身

変数名	データの説明
CRIM	犯罪発生率
ZN	住居区画の密集度
INDUS	非小売業の土地割合
CHAS	チャールズ川の周辺か否か
NOX	NOx濃度
RM	住居に含まれる平均部屋数
AGE	1940年より前に建てられた物件の割合
DIS	ボストン市の5つの雇用施設からの重み付き距離
RAD	大きな道路へのアクセスのしやすさ
TAX	1万ドルあたりの不動産税率の総計
PTRATIO	教師あたりの生徒数
B	対象地域に住む黒人の比率
LSTAT	低所得者の割合

上記の表に入っているデータを使って、住宅価格=MEDVを予測します。

今回は、RM(=住宅に含まれる平均部屋数)とLSTAT(=低所得者の割合)を使って、MEDV(=住宅価格)を予測しましょう。

散布図を使って可視化しておくと、以下のようになります。

X = df.RM
Y = df.LSTAT
Z = df.MEDV

fig = plt.figure()
ax = Axes3D(fig)

ax.set_xlabel('RM')
ax.set_ylabel('LSTAT')
ax.set_zlabel('MEDV')

ax.scatter(X,Y,Z, s=30, c='blue', alpha=0.2)

plt.show()

RM : 5
LSTAT : 20
MEDV : 10

上記の場所に集中しているように見えますね。

STEP③ : データの準備

X = df[['RM', 'LSTAT']]
y = df['MEDV']

print('-'*10 + '特徴量とターゲットに分割' + '-'*10)
print('X:', X.shape)
print('y:', y.shape)

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=0)

print('元データ', X.shape)
print('学習用データ', X_train.shape)
print('テスト用データ', X_test.shape)

それでは、重回帰モデルの実装で必要なデータの準備をしていきます。

特徴量(=機械学習で使う変数)はRMとLSTAT、予測対象はMEDVですね。

重回帰式で表すと、以下のようになりますね。

$${\rm MEDV} = w_0 + w_1 {\rm CRIM} + w_2 {\rm ZN} + \cdots + w_{13} {\rm LSTAT}$$

特徴量とターゲットの組み合わせを、Scikit-learnに含まれているtrain_test_split()を使って学習用データとテスト用データに分割します。

元々506行あったデータを、学習:テスト = 8:2で分割したので、以下のデータ数になります。

学習用データ : 404
テスト用データ : 102

STEP④ : 重回帰モデルの作成と評価

lm = LinearRegression()
lm.fit(X_train, y_train)

print('-'*10 + '学習結果' + '-'*10)
print('バイアス:', lm.intercept_)
print('重み:', lm.coef_)
print('-'*10 + 'モデルの評価' + '-'*10)
print('Train Score:', lm.score(X_train, y_train))
print('Test Score:', lm.score(X_test, y_test))

重回帰モデルの作成には、単回帰同様にLinearRegression()を使います。

単回帰は、重回帰モデルの$N=1$の場合なので、当たり前といえば当たり前ですね。

モデルの学習は、lm.fit()を使うだけで実行できます。

lm = LinearRegression()
lm.fit(X_train, y_train)

これで重回帰モデルの作成が完了しました。

あとは、学習結果を見てみましょう。lm.score()を使えばモデルの決定係数を確認できます。

print('Train Score：{:.2f}' .format(lm.score(X_train, y_train)))
print('Test Score：{:.2f}' .format(lm.score(X_test, y_test)))

Train Score：0.66
Test Score：0.54

決定係数は、0〜1の値を取り、1に近い方が精度が良いという解釈になります。

単回帰モデルのときは、

Train Score：0.50
Test Score：0.42

だったので、重回帰モデルを使うことで、精度はかなりよくなりましたね。

ちなみに、「良い重回帰モデルを作成する」とは、単回帰のときと同じで、重回帰式$y = w_0 + w_1 x_1+ w_2 x_2$に含まれる、$w_0$と$w_1$と$w_2$で最適な値を求めることになります。

実際に中身を確認してみましょう。

print('バイアス', lm.intercept_)
print('重み', lm.coef_)

バイアス -1.2625286844374664
重み [ 5.10906846 -0.65494879]

バイアスと呼んでいるのが、$w_0$、重みと呼んでいるのが$w_1$と$w_2$です。

数式で表すと、以下のようなモデルが完成したことになります。

$$y = -1.2625 + 5.1091 x_1 – 0.6550 x_2$$

STEP⑤ : 学習結果の可視化

a1, a2 = lm.coef_
b = lm.intercept_

#　実測値(散布図)の描画
X = df.RM
Y = df.LSTAT
Z = df.MEDV

fig = plt.figure()
ax = Axes3D(fig)
ax.scatter(X,Y,Z, s=30, c='blue', alpha=0.2)

# 回帰平面を描画
X_3d, Y_3d = np.meshgrid(np.arange(0, 10, 1), np.arange(0, 40, 1))
Z_3d = a1 * X_3d + a2 * Y_3d + b
ax.plot_surface(X_3d, Y_3d, Z_3d)
ax.set_xlabel('RM')
ax.set_ylabel('LSTAT')
ax.set_zlabel('MEDV')

ax.view_init(elev=60, azim=45)

plt.show()

元データと、今回作成した重回帰モデルの予測結果を確認してみます。