解析環境セットアップ

ライブラリーの準備

このノート全体にわたって、共通して使用するパッケージ、関数やデータなどをここで呼び出して、すぐに使える状態にする。

In [1]:
%matplotlib inline

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

データセットの準備

このノートでは scikit-learn パッケージに保存されているボストン住宅価格データセットを使って、回帰分析を行う例を示す。このデータセットは、各区画における住宅の価格、犯罪率、非小売業の割合、全額固定資産税率などを集計したデータとなっています。このページではこのデータセット使って、犯罪率、非小売業の割合、全額固定資産税率などの特徴量（説明変数）を使って、住宅の価格（目的変数）を予測する回帰分析を行う例を示す。

In [2]:
import sklearn.datasets
boston = sklearn.datasets.load_boston()

X = boston.data
y = boston.target

print(X.shape)

(506, 13)
In [3]:
boston_df = pd.DataFrame(np.concatenate([y.reshape(y.shape[0], 1), X], axis=1),
                         columns=['y' if i == 0 else 'X' + str(i) for i in range(X.shape[1] + 1)])
sns.pairplot(boston_df)
Out[3]:
<seaborn.axisgrid.PairGrid at 0x91dff7d50>
In [4]:
sns.pairplot(boston_df, y_vars=['y'], x_vars=boston_df.columns[1:])
Out[4]:
<seaborn.axisgrid.PairGrid at 0x1a26188c50>

線形回帰

scikit-learn で回帰分析（単回帰・重回帰）を行うには、線形モデルモジュール sklearn.linear_model 中の LinearRegression 関数を使用する。データの分布を眺めて、6 番目の特徴量と 13 番目の特徴量が、目的変数との相関が高そうであるから、まず、これらの特徴量で単回帰モデルを作成してみる。

In [5]:
import sklearn.model_selection
import sklearn.metrics
import sklearn.linear_model

X_subset = X[:, [5]]
X_train, X_test, y_train, y_test =\
    sklearn.model_selection.train_test_split(X_subset, y, test_size=0.3, random_state=2020)

reg = sklearn.linear_model.LinearRegression()
reg.fit(X_train, y_train)

y_pred = reg.predict(X_test)

教師ラベルを縦軸に、予測値を横軸に散布図で可視化して、モデルのテスト結果を確認してみる。グラフをみると、だいたい予測できているが、一部のサンプルで、実際の値よりも高く予測されているのがわかる。

In [6]:
plt.scatter(y_pred, y_test)
plt.plot([0, 50], [0, 50])
plt.grid()
plt.ylabel('true')
plt.xlabel('predicted')
plt.show()

教師ラベル（真の値）と予測値の差を残差とよび、各サンプルの残差をプロットした図を残差プロットという。上のモデルの予測結果を利用して残差プロットを描いてみる。残差プロットをみると、全体的に、予測した値は、実際の値よりも大きくなっているのを確認できる。

In [7]:
plt.scatter(X_test, y_test - y_pred)
plt.ylabel('residuals')
plt.xlabel('X')
plt.show()

この残差から MSE を計算して、モデル良さを測る指標を計算する。

In [8]:
mse = sklearn.metrics.mean_squared_error(y_test, y_pred)
print(mse)

51.699563454809955

次に 13 番目の特徴量だけを使って、住宅価格を予測してみる。

In [9]:
import sklearn.model_selection
import sklearn.metrics
import sklearn.linear_model
import sklearn.pipeline

X_subset = X[:, [12]]
X_train, X_test, y_train, y_test =\
    sklearn.model_selection.train_test_split(X_subset, y, test_size=0.3, random_state=2020)

reg = sklearn.linear_model.LinearRegression()
reg.fit(X_train, y_train)

y_pred = reg.predict(X_test)

mse = sklearn.metrics.mean_squared_error(y_test, y_pred)
print(mse)

35.69436606917185

MSE をみると一つ前のモデルよりも良いモデルができるこを確認できる。上と同様に、残差プロットを描いてみる。この残差プロットを見ても、全体的に、予測値が実際の値よりも大きくなる傾向が見られる。

In [10]:
plt.scatter(X_test, y_test - y_pred)
plt.ylabel('residuals')
plt.xlabel('X')
plt.show()

次に、この 13 番目の特徴量を対数化してからもう一度モデルを構築してみる。MSE をみると、対数化することで、MSE が大きく改善しているのが確認できる。また、残差プロットをみると、X の値が小さ時に誤差が大きいことがわかる。一つの特徴量でばんばるにはこれで限界かもしれない。

In [11]:
import sklearn.model_selection
import sklearn.metrics
import sklearn.linear_model
import sklearn.pipeline

X_subset = np.log10(X[:, [12]])
X_train, X_test, y_train, y_test =\
    sklearn.model_selection.train_test_split(X_subset, y, test_size=0.3, random_state=2020)

reg = sklearn.linear_model.LinearRegression()
reg.fit(X_train, y_train)

y_pred = reg.predict(X_test)
mse = sklearn.metrics.mean_squared_error(y_test, y_pred)
print(mse)

plt.scatter(X_test, y_test - y_pred)
plt.ylabel('residuals')
plt.xlabel('X')
plt.show()

29.22244616392045

同じ関数で単回帰・重回帰モデルも構築することは可能である。重回帰分析の場合も、コードを変更せずに、特徴量を大きな行列として与えればよい。もちろん、特徴量が多いほど、バリアンスの大きいモデル（複雑なモデル）が構築される傾向が見えるので、予測誤差が小さくなっていきます。また、重回帰分析のとき、各変数間のスケールが異なるので、スケールを揃えるために標準化を行う必要がある。そのため、これ以降パイプラインを使ってモデルを構築していくことにする。

In [12]:
import sklearn.model_selection
import sklearn.metrics
import sklearn.linear_model
import sklearn.pipeline


X_train, X_test, y_train, y_test =\
    sklearn.model_selection.train_test_split(X, y, test_size=0.3, random_state=2020)


pipe = sklearn.pipeline.Pipeline([
     ('scaler', sklearn.preprocessing.StandardScaler()),
     ('reg', sklearn.linear_model.LinearRegression())
])


pipe.fit(X_train, y_train)
y_pred = pipe.predict(X_test)

mse = sklearn.metrics.mean_squared_error(y_test, y_pred)
print(mse)

23.44040354365431

13 番目の特徴量を使うだけで MSE が 29.222 であるのに対して、特徴量を 13 個全部使っての MSE が 23.44 であり、それほど改善しているように見えない。この回帰モデルの各係数は、次のようにして取得することができる。

回帰モデル中の各パラメーター（各係数）は、pipe.intercept_ および pipe.coef_ からアクセスすることができる。

In [13]:
print(pipe['reg'].intercept_)
print(pipe['reg'].coef_)

22.90706214689265
[-0.9116096   1.01328621  0.07548381  0.83939401 -1.72435422  3.09130242
 -0.28391525 -3.00957316  2.52848895 -1.84778323 -1.96638394  0.77217247
 -3.64295986]

残差プロットを描いて、残差の分布を確認してみる。特徴量全部使っても、予測値が、真の値よりもやや高くなっているのがわかる。おそらく、学習データの中に外れ値的なものが存在して、それに影響されていると思われる。

In [14]:
plt.scatter(np.arange(0, len(y_test)), y_test - y_pred)
plt.ylabel('residuals')
plt.xlabel('test sample id')
plt.show()

スパース回帰

LASSO を使用すると、変数選択を自動的に行ってくれる。例えば、次のように 13 特徴量全体を使ってモデルを構築することからスタートして、LASSO 実行して、その係数を見てみる。

In [15]:
import sklearn.model_selection
import sklearn.metrics
import sklearn.linear_model
import sklearn.pipeline


X_train, X_test, y_train, y_test =\
    sklearn.model_selection.train_test_split(X, y, test_size=0.3, random_state=2020)


pipe = sklearn.pipeline.Pipeline([
     ('scaler', sklearn.preprocessing.StandardScaler()),
     ('reg', sklearn.linear_model.Lasso(alpha=0.01))
])


pipe.fit(X_train, y_train)
y_pred = pipe.predict(X_test)

mse = sklearn.metrics.mean_squared_error(y_test, y_pred)
print(mse)

23.463664332788266
In [16]:
print(pipe['reg'].coef_)
print(pipe['reg'].intercept_)

[-0.88414383  0.97532977  0.          0.83950241 -1.66881949  3.09948388
 -0.26148148 -2.95331587  2.36582953 -1.69476384 -1.94743695  0.76402207
 -3.64182994]
22.90706214689265
In [17]:
import sklearn.model_selection
import sklearn.metrics
import sklearn.linear_model
import sklearn.pipeline


X_train, X_test, y_train, y_test =\
    sklearn.model_selection.train_test_split(X, y, test_size=0.3, random_state=2020)


pipe = sklearn.pipeline.Pipeline([
     ('scaler', sklearn.preprocessing.StandardScaler()),
     ('reg', sklearn.linear_model.Lasso(alpha=1))
])


pipe.fit(X_train, y_train)
y_pred = pipe.predict(X_test)

mse = sklearn.metrics.mean_squared_error(y_test, y_pred)
print(mse)

print(pipe['reg'].coef_)
print(pipe['reg'].intercept_)

29.488389719749218
[-0.          0.         -0.          0.16468261 -0.          3.14865468
 -0.         -0.         -0.         -0.         -1.3069594   0.1489821
 -3.30386633]
22.907062146892642

LASSO のハイパーパラメーターとして alpha がある。この alpha を動かして、最小 MSE を与える alpha を求めて、そのときのモデルの性能を評価してみよう。ハイパーパラメーターを探す時は、5-fold 交差検証を使用してみるとよい。

In [18]:
import sklearn.model_selection
import sklearn.metrics
import sklearn.linear_model
import sklearn.pipeline


X_train, X_test, y_train, y_test =\
  sklearn.model_selection.train_test_split(X, y, test_size=0.4, random_state=2020)


pipe = sklearn.pipeline.Pipeline([
         ('scaler', sklearn.preprocessing.StandardScaler()),
         ('reg', sklearn.linear_model.Lasso())
       ])

parameters = {
    'reg__alpha': 10**np.linspace(-5, 2, 10)
}


# GrideSearchCV で scoring='mse' を指定できないので、
# (R)MSE を計算する関数を自分で定義して GrideSearchCV に代入することにする
def calc_rmse(predict, actual):
    return np.sqrt(sklearn.metrics.mean_squared_error(predict, actual))
rmse_score = sklearn.metrics.make_scorer(calc_rmse, greater_is_better = False)

gs = sklearn.model_selection.GridSearchCV(pipe, parameters, scoring=rmse_score, cv=10)
gs.fit(X_train, y_train)


best_pipe = gs.best_estimator_
print(best_pipe['reg'].coef_)

y_pred = best_pipe.predict(X_test)

mse = sklearn.metrics.mean_squared_error(y_test, y_pred)
print(mse)

[-0.70897493  0.78613098 -0.0181187   0.75523844 -1.23984228  3.33432731
 -0.         -2.45040548  1.44180597 -0.95908091 -1.7952367   0.8391347
 -3.84376342]
21.176559196279555

/opt/anaconda3/lib/python3.7/site-packages/sklearn/model_selection/_search.py:814: DeprecationWarning: The default of the `iid` parameter will change from True to False in version 0.22 and will be removed in 0.24. This will change numeric results when test-set sizes are unequal.
  DeprecationWarning)

次に、5-fold 交差検証を用いて、最適なハイパーパラメーターで elastic net 回帰モデルを構築し、その性能を評価してみよう。

In [19]:
import sklearn.model_selection
import sklearn.metrics
import sklearn.linear_model
import sklearn.pipeline


X_train, X_test, y_train, y_test =\
  sklearn.model_selection.train_test_split(X, y, test_size=0.4, random_state=2020)


pipe = sklearn.pipeline.Pipeline([
         ('scaler', sklearn.preprocessing.StandardScaler()),
         ('reg', sklearn.linear_model.ElasticNet())
       ])

parameters = {
    'reg__alpha': 10**np.linspace(-5, 2, 10),
    'reg__l1_ratio': [0.1, 0.2, 0.3, 0.4, 0.5],
}


def calc_rmse(predict, actual):
    return np.sqrt(sklearn.metrics.mean_squared_error(predict, actual))


rmse_score = sklearn.metrics.make_scorer(calc_rmse, greater_is_better = False)


gs = sklearn.model_selection.GridSearchCV(pipe, parameters, scoring=rmse_score, cv=10)
gs.fit(X_train, y_train)


best_pipe = gs.best_estimator_
print(best_pipe['reg'].coef_)

y_pred = best_pipe.predict(X_test)

mse = sklearn.metrics.mean_squared_error(y_test, y_pred)
print(mse)

[-0.7500703   0.76209903 -0.20621511  0.79441487 -1.10702076  3.30047096
 -0.06756445 -2.32221787  1.34371409 -0.88813681 -1.75604804  0.83846338
 -3.58392839]
21.164491603474023

/opt/anaconda3/lib/python3.7/site-packages/sklearn/model_selection/_search.py:814: DeprecationWarning: The default of the `iid` parameter will change from True to False in version 0.22 and will be removed in 0.24. This will change numeric results when test-set sizes are unequal.
  DeprecationWarning)

ロバスト推定

ロバスト推定アルゴリズムの一つ RANSAC は RANSACRegressor で実行することができる。この推定では、様々なハイパーパラメーターを設定していく必要がある。以下の例は、すべてのハイパーパラメーターをデフォルトのままでモデルを構築している。

In [20]:
import sklearn.model_selection
import sklearn.metrics
import sklearn.linear_model
import sklearn.pipeline


X_train, X_test, y_train, y_test =\
    sklearn.model_selection.train_test_split(X, y, test_size=0.3, random_state=2020)


pipe = sklearn.pipeline.Pipeline([
     ('scaler', sklearn.preprocessing.StandardScaler()),
     ('reg', sklearn.linear_model.RANSACRegressor(random_state=2020))
])


pipe.fit(X_train, y_train)
y_pred = pipe.predict(X_test)

mse = sklearn.metrics.mean_squared_error(y_test, y_pred)
print(mse)

32.33166450489739

scikit-learn に実装されているロバスト推定アルゴリズムを調べてみると、スライドで紹介したRANSAC 以外のも、HuberRegressor および HuberRegressor がある。これらのアルゴリズムの詳細は Google 検索で。

In [21]:
import sklearn.model_selection
import sklearn.metrics
import sklearn.linear_model
import sklearn.pipeline


X_train, X_test, y_train, y_test =\
    sklearn.model_selection.train_test_split(X, y, test_size=0.3, random_state=2020)


pipe = sklearn.pipeline.Pipeline([
     ('scaler', sklearn.preprocessing.StandardScaler()),
     ('reg', sklearn.linear_model.HuberRegressor())
])


pipe.fit(X_train, y_train)
y_pred = pipe.predict(X_test)

mse = sklearn.metrics.mean_squared_error(y_test, y_pred)
print(mse)

24.724846128957655
In [ ]: