Python: Boston データセットで線形回帰分析を学ぶ
今回は実践機械学習システムの第七章を参考にして、線形回帰分析について学んでみる。 使用する Boston データセットというのは、ボストンの物件の価格にその物件の人口統計に関する情報が付随したものだ。 つまり、ある物件の人口統計に関する情報を元に、その物件の価格を予測することになる。 線形回帰というのは説明変数を元に応答が一次関数で表されると仮定したモデルを回帰に使用するものだ。 ようするに、例えば部屋の数が多い住宅はその分価格も高いのではないか、という感じで説明変数が目的変数といい感じに比例する関係にあるんじゃねと考える。
今回は線形回帰分析に scipy と scikit-learn を、グラフ描画に matplotlib を使う。
$ pip install scipy scikit-learn matplotlib
まずは、一つの説明変数を元にした (一次元の) 回帰から取り掛かる。 先ほどの予測を元に、データセットの中の部屋の数から住宅価格を回帰分析してみよう。 以下のサンプルコードでは x 軸に部屋の数、そして y 軸に住宅の価格をとって、そこから最小二乗法で最も誤差の小さくなる直線を得ている。 また、得られた直線からこのアルゴリズムを評価するために平均二乗平方根誤差 (Root Mean Squared Error: RMSE) という値を計算している。 訓練誤差の RMSE は訓練に使った教師信号と、得られた直線との間にどれだけの違いがあるかを示している。
#!/usr/bin/env python # -*- coding: utf-8 -*- import numpy as np from sklearn import datasets from matplotlib import pyplot as plt def main(): # ボストンデータセットを読み込む boston = datasets.load_boston() # 部屋の数を取り出す rooms = boston.data[:, 5] # 家の値段を取り出す house_prices = boston.target # 部屋の数と家の値段の関係をプロットする plt.scatter(rooms, house_prices, color='r') # 最小二乗法で誤差が最も少なくなる直線を得る x = np.array([[v] for v in rooms]) y = house_prices slope, total_error, _, _ = np.linalg.lstsq(x, y) # 得られた直線をプロットする plt.plot(x, slope * x) # 訓練誤差の平均二乗平方根誤差 (Root Mean Squared Error: RMSE) rmse = np.sqrt(total_error[0] / len(x)) msg = 'RMSE (training): {0}'.format(rmse) print(msg) # グラフを表示する plt.xlabel('部屋の数') plt.ylabel('家の値段 (単位: 1000 ドル)') plt.autoscale() plt.grid() plt.show() if __name__ == '__main__': main()
上記を実行すると以下のようなグラフが得られる。 見た感じ直線がデータに上手くフィットしていないようだ。
RMSE の結果は以下の通り。
RMSE (training): 7.64268509308749
先ほどの回帰分析では単純に部屋の数と住宅価格が比例すると考えていた。 しかし、部屋の数が 0 の家があったら住宅価格も 0 になるのかというとそうではないだろう。 つまり、部屋の数に関わりなく最低限必要な価格もあるんじゃないかということになる。 次は、その最低限必要な価格をバイアス項として計算に付加してみることにする。
#!/usr/bin/env python # -*- coding: utf-8 -*- import numpy as np from sklearn import datasets from matplotlib import pyplot as plt def main(): # ボストンデータセットを読み込む boston = datasets.load_boston() # 部屋の数を取り出す rooms = boston.data[:, 5] # 家の値段を取り出す house_prices = boston.target # 部屋の数と家の値段の関係をプロットする plt.scatter(rooms, house_prices, color='r') # 最小二乗法で誤差が最も少なくなる直線を得る x = np.array([[v, 1] for v in rooms]) # バイアス項を追加する y = house_prices (slope, bias), total_error, _, _ = np.linalg.lstsq(x, y) # 得られた直線をプロットする plt.plot(x[:, 0], slope * x[:, 0] + bias) # 訓練誤差の RMSE rmse = np.sqrt(total_error[0] / len(x)) msg = 'RMSE (training): {0}'.format(rmse) print(msg) # グラフを表示する plt.xlabel('部屋の数') plt.ylabel('家の値段 (単位: 1000 ドル)') plt.grid() plt.show() if __name__ == '__main__': main()
上記を実行した際に得られるグラフは以下の通り。 先ほどに比べると直線がデータに上手くフィットしているように見える。
RMSE は以下の通り。 先ほどと比べると、訓練誤差の RMSE が小さくなっていることが分かる。
RMSE (training): 6.6030713892225625
ここまでは部屋の数という一次元の説明変数を元にした回帰を扱ってきた。 次は Boston データセットに含まれる全ての説明変数を一度に使って、多次元 (具体的には 14 次元) の回帰を試してみよう。
#!/usr/bin/env python # -*- coding: utf-8 -*- import numpy as np from sklearn import datasets def main(): # ボストンデータセットを読み込む boston = datasets.load_boston() # 家の値段を取り出す house_prices = boston.target # 最小二乗法で誤差が最も少なくなる直線を得る x = np.array([np.concatenate((v, [1])) for v in boston.data]) # バイアス項を追加する y = house_prices _slope, total_error, _, _ = np.linalg.lstsq(x, y) # 訓練誤差の RMSE rmse = np.sqrt(total_error[0] / len(x)) msg = 'RMSE (training): {0}'.format(rmse) print(msg) # 多次元 (14 次元) の回帰なのでグラフで可視化できない if __name__ == '__main__': main()
上記を実行して得られるグラフは…と言いたいところだけど 14 次元はグラフに可視化できない。 RMSE に関しては以下の通り。 先ほどに比べると訓練誤差の RMSE がまた一段と小さくなっている。
RMSE (training): 4.679506300635518
尚、ここまでは scipy を使って線形回帰しているが、scikit-learn を使っても同様のことができる。 上記を scikit-learn を使って書きなおしたものが以下のサンプルコードだ。 また、ここからはK分割交差検証を使うことでアルゴリズムの汎化性能 (未知のデータに対処する能力) も測ってみよう。 K分割交差検証では、教師信号の一部をテストデータとして学習に使用せずにおく。 そして、学習した結果をテストデータを使って検証することで、アルゴリズムの未知のデータに対処する能力を測ることができる。
#!/usr/bin/env python # -*- coding: utf-8 -*- import numpy as np from sklearn import datasets from sklearn import linear_model from sklearn import cross_validation def main(): # ボストンデータセットを読み込む boston = datasets.load_boston() # 家の値段を取り出す house_prices = boston.target # 説明変数 + バイアス項 x = np.array([np.concatenate((v, [1])) for v in boston.data]) # 家の値段 y = house_prices # 線形回帰する lr = linear_model.LinearRegression(fit_intercept=True) lr.fit(x, y) # 訓練誤差の RMSE p = np.array([lr.predict(xi) for xi in x]) error = p - y total_error = np.sum(error * error) rmse_train = np.sqrt(total_error / len(p)) # K-分割交差検証誤差の RMSE kf = cross_validation.KFold(len(x), n_folds=10) error = 0 for training, test in kf: lr.fit(x[training], y[training]) p = np.array([lr.predict(xi) for xi in x[test]]) e = p - y[test] error += np.sum(e * e) rmse_10cv = np.sqrt(error / len(x)) # 結果を表示する msg = 'RMSE (training): {0}'.format(rmse_train) print(msg) msg = 'RMSE (10-fold CV): {0}'.format(rmse_10cv) print(msg) if __name__ == '__main__': main()
実行結果は以下の通り。 訓練誤差の RMSE に比べるとK分割交差検証の RMSE が大きいようだ。
RMSE (training): 4.679506300635516 RMSE (10-fold CV): 5.881925072430106
線形回帰分析では、説明変数が多く教師信号が少ないほど過学習を起こしやすいらしい。 過学習が起こると訓練誤差は小さくなる反面汎化性能は落ちることになる。 そこで、罰則項 (または正則化項) という考えを導入することで過学習を防ぐ回帰アルゴリズムとして Ridge 回帰、Lasso 回帰、そして Elastic net という手法が存在する。 以下のサンプルコードでは各アルゴリズムを使うとどのような結果になるか試してみた。
#!/usr/bin/env python # -*- coding: utf-8 -*- import numpy as np from sklearn import datasets from sklearn import linear_model from sklearn import cross_validation def main(): # ボストンデータセットを読み込む boston = datasets.load_boston() # 説明変数 x = boston.data # 家の値段 y = boston.target # 線形回帰する models = [ linear_model.ElasticNetCV(fit_intercept=True), linear_model.RidgeCV(fit_intercept=True), linear_model.LassoCV(fit_intercept=True), ] for model in models: # 訓練データにそのままフィットさせる model.fit(x, y) # 訓練誤差の RMSE p = np.array([model.predict(xi) for xi in x]) error = p - y total_error = np.sum(error * error) rmse_train = np.sqrt(total_error / len(p)) # K-分割交差検証誤差の RMSE kf = cross_validation.KFold(len(x), n_folds=10) error = 0 for training, test in kf: model.fit(x[training], y[training]) p = np.array([model.predict(xi) for xi in x[test]]) e = p - y[test] error += np.sum(e * e) rmse_10cv = np.sqrt(error / len(x)) # 結果を表示する msg = 'Model: {0}'.format(model.__class__.__name__) print(msg) msg = 'RMSE (training): {0}'.format(rmse_train) print(msg) msg = 'RMSE (10-fold CV): {0}'.format(rmse_10cv) print(msg) print() if __name__ == '__main__': main()
実行結果は以下の通り。 LassoCV では単純な線形回帰よりも汎化性能が落ちてしまっているが、RidgeCV では同程度、ElasticNetCV では向上している。 注目すべきは、ElasticNetCV において、訓練誤差の値は先ほどより悪化している点だろう。 訓練誤差が小さいということは、より教師信号にフィットしているという意味になる。 つまり、ElasticNetCV では単純な線形回帰よりも教師信号にはフィットしない反面、未知のデータに対処する能力を上げることができている。 訓練誤差がいかに小さくとも、実際のデータにアルゴリズムを適用する際に重要なのはあくまでも汎化性能なので、Elastic netのような手法は貴重と言えるだろう。
Model: ElasticNetCV RMSE (training): 5.267971563621094 RMSE (10-fold CV): 5.697230789950297 Model: RidgeCV RMSE (training): 4.679889836546443 RMSE (10-fold CV): 5.870035920569123 Model: LassoCV RMSE (training): 5.012392434381019 RMSE (10-fold CV): 6.181841135046689
