10. 回帰分析---連続値を取る目的変数の予測¶
- 本章は「教師あり学習」の一つの分野「回帰分析(regression analysis)」を取り上げる
- 回帰モデル： 連続値をとる目的変数の予測
そのため、変数間の関係を理解、データの傾向を評価、予報を作成するという形での応用
例：企業の今後数カ月の月間売り上げの予測 - 内容
- データセットの探求と可視化
- 線形回帰モデルの実装のためのアプローチ
- 外れ値に対する頑健な回帰モデルの訓練
- 回帰モデルの評価とモデルに共通する問題の診断
- 回帰モデルを非線形データにあてはめる
In [1]:
from IPython.display import Image
%matplotlib inline
10.1 線形回帰¶
10.1.1 単線形回帰¶
- 単純な（単変量の）線形回帰の目的：単一の特徴量（説明変数x）と連続値の応答(目的変数y)との関係をモデルとして表現すること
- 説明変数が一つの線形モデルの方程式：$y = w_0+ w_1x$
重み $w_0$ はy軸の切片、$w_1$は説明変数$x$の係数
目標は、この1次式の重みを学習すること⇒未知データに対する応答を予測 - 1次式を前提とすると、線形回帰は「サンプル点を通過する直線のうち最も適合する直線を見つけ出すこと」
In [2]:
Image(filename='images/10_01.png', width=500)
Out[2]:
- 回帰直線(regression line): 最も適合する直線
オフセット、残渣(residual): 回帰直線からサンプル点への縦線―予測の誤差を表す
- 重回帰(multiple linear regression): 複数の説明変数による線形回帰($x_0 = 1$とする)
$\displaystyle y = w_0x_0 + w_1x_1 + \cdots + w_mx_m = \sum^m_{i=0} w_ix_i = \boldsymbol{w}^T\boldsymbol{x}$
In [3]:
Image(filename='images/10_15.png', width=500)
Out[3]:
10.2 Housingデータセットの探求¶
10.2.1 Housingデータセットをデータフレームに読み込む¶
- Harrison & Rubinfeld (1978): Housingデータセット
506個のボストン近郊の住宅情報、1個の目的変数(MEDV, 住宅価格の中央値)、13個の説明変数：
- CRIM per capita crime rate by town
- ZN proportion of residential land zoned for lots over
25,000 sq.ft. - INDUS proportion of non-retail business acres per town
- CHAS Charles River dummy variable (= 1 if tract bounds
river; 0 otherwise) - NOX nitric oxides concentration (parts per 10 million)
- RM average number of rooms per dwelling
- AGE proportion of owner-occupied units built prior to 1940
- DIS weighted distances to five Boston employment centres
- RAD index of accessibility to radial highways
- TAX full-value property-tax rate per $10,000
- PTRATIO pupil-teacher ratio by town
- B 1000(Bk - 0.63)^2 where Bk is the proportion of blacks
by town - LSTAT % lower status of the population
- MEDV Median value of owner-occupied homes in $1000s </pre>
- UCIリポジトリからpandasのDataFrameオブジェクトとして取り出す：
In [1]:
import pandas as pd
df = pd.read_csv('https://raw.githubusercontent.com/rasbt/'
'python-machine-learning-book-2nd-edition'
'/master/code/ch10/housing.data.txt',
header=None,
sep='\s+')
df.columns = ['CRIM', 'ZN', 'INDUS', 'CHAS',
'NOX', 'RM', 'AGE', 'DIS', 'RAD',
'TAX', 'PTRATIO', 'B', 'LSTAT', 'MEDV']
df.head()
Out[1]:
10.2.2 データセットの重要な特性を可視化する¶
- 探索的データ解析(EDA, Exploratory Data Analysis): 機械学習モデルの訓練を行う前の、最初の重要なステップ
- 散布図行列の作成：データセットの特徴量のペアに対する相関関係を一つの平面で可視化可能
seabornライブラリのpairplot関数を使用（installが必要、conda install seaborn）:
In [2]:
import matplotlib.pyplot as plt
import seaborn as sns
In [3]:
cols = ['LSTAT', 'INDUS', 'NOX', 'RM', 'MEDV']
sns.pairplot(df[cols], size=2.5)
plt.tight_layout()
# plt.savefig('images/10_03.png', dpi=300)
plt.show()
この散布図行列から変数間の関係が確認できる：RMとMEDVが正の相関、MEDVは正規分布＋外れ値
- 相関行列の作成：標準化されたデータから計算された共分散行列に等しい
10.2.3 相関行列を用いて関係を調べる¶
- 相関行列: ピアソンの積率相関係数--ピアソンの$r$(Pearson's r)--を成分とする正方行列
相関係数の値は-1から1の範囲、r=1なら完全に正の相関、r=-1なら負の相関:
$\displaystyle r = \frac{\sum^n_{i=1}[(x^{(i)}-\mu_x)(y^{(i)}-\mu_y)]}{\sqrt{\sum^n_{i=1}(x^{(i)}-\mu_x)^2}\sqrt{\sum^n_{i=1}(y^{(i)}-\mu_y)^2}} = \frac{\sigma_{xy}}{\sigma_x\sigma_y}$
ここで$\mu$は対応する特徴量の標本平均、$\sigma_{xy}$は特徴量$x$と$y$の共分散、$\sigma_x$は特徴量$x$の標準偏差
- Numpyの
corrcoef関数によりピアソンの積率相関係数を計算
seabornのheatmap関数を用いて、相関行列をヒートマップとしてプロット
特徴量の線形相関に基づき特徴量を選択するのに役立つ
In [9]:
import numpy as np
cm = np.corrcoef(df[cols].values.T)
#sns.set(font_scale=1.5)
hm = sns.heatmap(cm,
cbar=True,
annot=True,
square=True,
fmt='.2f',
annot_kws={'size': 15},
yticklabels=cols,
xticklabels=cols)
plt.tight_layout()
# plt.savefig('images/10_04.png', dpi=300)
plt.show()
- 線形回帰モデルの適合には、目的変数(MEDV)との相関が高い特徴量に注目
LSTATが-0.74 (散布図行列から非線形の関係)、RMが0.70
10.3 最小二乗線形回帰モデルの実装¶
- 最小二乗(OLS, Ordinary Least Squares)法: サンプル点に対する縦の距離（残差）の二乗の和を最小化するパラメータを推定する
10.3.1 勾配降下法を用いて回帰パラメータの回帰を求める¶
- ADALINE(2章)のコスト関数は、誤差平方和(SSE) $\displaystyle J(\boldsymbol{w}) = \frac{1}{2} \sum^n_{i=1}(y^{(i)}-\hat{y}^{(i)})^2$
最小二乗法(OLS)のコスト関数に等しい
最小二乗線形回帰は、基本的に単位ステップ関数のないADALINEとみなせる - これを示すために、ADALINEの勾配降下法の実装から単位ステップ関数を削除することで、線形回帰モデルを実装してみる
In [6]:
# 基本的な線形回帰モデル
class LinearRegressionGD(object):
def __init__(self, eta=0.001, n_iter=20):
self.eta = eta
self.n_iter = n_iter
def fit(self, X, y):
self.w_ = np.zeros(1 + X.shape[1])
self.cost_ = []
for i in range(self.n_iter):
output = self.net_input(X)
errors = (y - output)
self.w_[1:] += self.eta * X.T.dot(errors)
self.w_[0] += self.eta * errors.sum()
cost = (errors**2).sum() / 2.0
self.cost_.append(cost)
return self
def net_input(self, X):
return np.dot(X, self.w_[1:]) + self.w_[0]
def predict(self, X):
return self.net_input(X)
- LinearRegressionGD回帰器の動作確認のためHousingデータセットのMEDVを予測するモデルを学習してみる: 説明変数としてRMを使用
In [7]:
X = df[['RM']].values
y = df['MEDV'].values
In [8]:
y
Out[8]:
In [10]:
from sklearn.preprocessing import StandardScaler
sc_x = StandardScaler()
sc_y = StandardScaler()
X_std = sc_x.fit_transform(X)
y_std = sc_y.fit_transform(y[:, np.newaxis]).flatten()
In [11]:
lr = LinearRegressionGD()
lr.fit(X_std, y_std)
Out[11]:
- 勾配降下法などの最適化アルゴリズムを使用する場合は、コストの最小値に収束することを確認するため、コストをエポック数の関数としてプロットする
In [11]:
plt.plot(range(1, lr.n_iter+1), lr.cost_)
plt.ylabel('SSE')
plt.xlabel('Epoch')
#plt.tight_layout()
#plt.savefig('images/10_05.png', dpi=300)
plt.show()
勾配降下法のアルゴリズムは5回目のエポックの後に収束している
- 線形回帰の直線がどの程度訓練データに適合しているかを可視化してみる
In [12]:
def lin_regplot(X, y, model):
plt.scatter(X, y, c='steelblue', edgecolor='white', s=70)
plt.plot(X, model.predict(X), color='black', lw=2)
return
In [15]:
lin_regplot(X_std, y_std, lr)
plt.xlabel('Average number of rooms [RM] (standardized)')
plt.ylabel('Price in $1000s [MEDV] (standardized)')
#plt.savefig('images/10_06.png', dpi=300)
plt.show()
グラフからわかるように、部屋の数と住宅価格の相関が高いことを示している
ただし、説明しきれない場合が多いことも示している($y=3$でデータ列が一列に並んでいる---住宅価格がそこで打ち切られている?)
In [16]:
print('Slope: %.3f' % lr.w_[1])
print('Intercept: %.3f' % lr.w_[0])
上は、標準化されたデータを扱っているため、y切片が0になる例を示している
結果変数の尺度を元の尺度に戻す必要がある場合の操作（この例では価格に戻す）
In [17]:
num_rooms_std = sc_x.transform(np.array([[5.0]]))
price_std = lr.predict(num_rooms_std)
print("Price in $1000s: %.3f" % sc_y.inverse_transform(price_std))
部屋数が5の住宅の価格は10,840ドルとなる
10.3.2 scikit-learnを用いて回帰モデルの係数を推定¶
- scikit-learnのLinearRegressionオブジェクト：標準化されていない変数に適応する
LIBLINEARライブラリや最適化されたアルゴリズムを採用
In [13]:
from sklearn.linear_model import LinearRegression
slr = LinearRegression()
slr.fit(X, y)
y_pred = slr.predict(X)
print('Slope: %.3f' % slr.coef_[0])
print('Intercept: %.3f' % slr.intercept_)
- 10.3.1では標準化したため傾きが0.695, 切片が0.0であった
それに対して異なるモデル係数を出力
RMに対するMEDVをプロットすることで、前の結果と比較
In [19]:
lin_regplot(X, y, slr)
plt.xlabel('Average number of rooms [RM]')
plt.ylabel('Price in $1000s [MEDV]')
#plt.savefig('images/10_07.png', dpi=300)
plt.show()
参考: 最小二乗法には、連立1次方程式を利用する方法もある: 正規方程式$\boldsymbol{w}=(\boldsymbol{X}^T\boldsymbol{X})^{-1}\boldsymbol{X}^T\boldsymbol{y}$
- Pythonによる実装：
この方法は最適解を解析的に求めることが保証されているが、逆行列の計算コストが高いXb = np.hstack((np.ones((X.shape[0],1)),X) w = np.zeros(X.shape[1]) z = np.linalg.inv(np.dot(Xb.T, Xb)) # (X^T X)^{-1} w = np.dot(z, np.dot(Xb.T, y)) # Xb^T y print(‘Slope: %.3f’ % w[0]) # 9.102 print(‘Intercept: %.3f’ % w[0]) # -34.671
10.4 RANSACを用いたロバスト回帰モデルの学習¶
- 線形回帰モデルは外れ値によって影響を受ける
外れ値の検出や除去は専門知識が必要 - 外れ値の除去に代わる方法： RANSAC (RANdom Sample Consensus)アルゴリズム
回帰モデルにデータのサブセット（正常値）を学習させる- 正常値としてランダムな数のサンプルを選択し、モデルを学習させる
- 学習済みのモデルに対し、そのほかのすべてのデータ点を評価し、ユーザー指定の許容範囲となるデータ点を正常値に追加
- すべての正常値を用いてモデルを再学習
- モデルの性能を評価、ユーザー指定の閾値の条件を満たしているか、繰り返し回数が規定回数に達した場合は終了。さもなくばステップ１に戻る
- scikit-learnのRANSACRegressorオブジェクトを使用
In [14]:
from sklearn.linear_model import RANSACRegressor
ransac = RANSACRegressor(LinearRegression(),
max_trials=100,
min_samples=50,
loss='absolute_loss',
residual_threshold=5.0,
random_state=0)
ransac.fit(X, y)
inlier_mask = ransac.inlier_mask_
outlier_mask = np.logical_not(inlier_mask)
line_X = np.arange(3, 10, 1)
line_y_ransac = ransac.predict(line_X[:, np.newaxis])
plt.scatter(X[inlier_mask], y[inlier_mask],
c='steelblue', edgecolor='white',
marker='o', label='Inliers')
plt.scatter(X[outlier_mask], y[outlier_mask],
c='limegreen', edgecolor='white',
marker='s', label='Outliers')
plt.plot(line_X, line_y_ransac, color='black', lw=2)
plt.xlabel('Average number of rooms [RM]')
plt.ylabel('Price in $1000s [MEDV]')
plt.legend(loc='upper left')
#plt.savefig('images/10_08.png', dpi=300)
plt.show()
In [15]:
print('Slope: %.3f' % ransac.estimator_.coef_[0])
print('Intercept: %.3f' % ransac.estimator_.intercept_)
- RANSACを使えば、データセットの外れ値の影響は抑えられる
- 回帰モデルをいろいろな方法で評価することは重要 --- 次節
In [15]:
# 参考 residual_threshold=10.0に変更してみる
from sklearn.linear_model import RANSACRegressor
ransac = RANSACRegressor(LinearRegression(),
max_trials=100,
min_samples=50,
loss='absolute_loss',
residual_threshold=10.0,
random_state=0)
ransac.fit(X, y)
inlier_mask = ransac.inlier_mask_
outlier_mask = np.logical_not(inlier_mask)
line_X = np.arange(3, 10, 1)
line_y_ransac = ransac.predict(line_X[:, np.newaxis])
plt.scatter(X[inlier_mask], y[inlier_mask],
c='steelblue', edgecolor='white',
marker='o', label='Inliers')
plt.scatter(X[outlier_mask], y[outlier_mask],
c='limegreen', edgecolor='white',
marker='s', label='Outliers')
plt.plot(line_X, line_y_ransac, color='black', lw=2)
plt.xlabel('Average number of rooms [RM]')
plt.ylabel('Price in $1000s [MEDV]')
plt.legend(loc='upper left')
plt.show()
In [16]:
print('Slope: %.3f' % ransac.estimator_.coef_[0])
print('Intercept: %.3f' % ransac.estimator_.intercept_)
10.5 線形回帰モデルの性能評価¶
- モデルの性能を偏りなく推定するには、学習に使用されなかったデータでモデルを検証する必要がある
- 6章で述べたように、データセットを分割し、訓練データセットをモデルの学習に、テストデータセットを未知データに対する汎化性能評価に私用する
- ここではデータセットの変数すべてを使用して重回帰モデルを学習させてみる
In [16]:
from sklearn.model_selection import train_test_split
X = df.iloc[:, :-1].values
y = df['MEDV'].values
X_train, X_test, y_train, y_test = train_test_split(
X, y, test_size=0.3, random_state=0)
slr = LinearRegression()
slr.fit(X_train, y_train)
y_train_pred = slr.predict(X_train)
y_test_pred = slr.predict(X_test)
- 残差プロット_(residual plot): 線形回帰モデルを診断するため、予測された値に対する残差、つまり実際の値と予測値との差のプロット
回帰モデルを診断し、非線形性や外れ値を検出し、誤差がランダムに分布しているかどうかをチェックするためによく使用される
註: 線形回帰モデルでは以下の仮定が置かれており、データにモデルを適合させた後で残差がこの条件を満たしているかどうかの確認が必要:
1. 誤差の期待値が0
2. 誤差は互いに無相関
3. 誤差の分散が等しい(等分散、等質性)In [17]:
plt.scatter(y_train_pred, y_train_pred - y_train,
c='steelblue', marker='o', edgecolor='white',
label='Training data')
plt.scatter(y_test_pred, y_test_pred - y_test,
c='limegreen', marker='s', edgecolor='white',
label='Test data')
plt.xlabel('Predicted values')
plt.ylabel('Residuals')
plt.legend(loc='upper left')
plt.hlines(y=0, xmin=-10, xmax=50, color='black', lw=2)
plt.xlim([-10, 50])
plt.tight_layout()
# plt.savefig('images/10_09.png', dpi=300)
plt.show()
- 予測が完璧ならば、残差はちょうど0となる---現実のデータではありえない
- よい回帰モデルでは、誤差がランダムに分布し、残差が中央の直線のまわりに散らばる
何らかのパタンが見られる場合は、モデルが何らかの情報を補足できていないことを意味
この例でもそのようなパタンがある - 外れ値の検出も可能---中央の直線から大きく離れた点が外れ値
- 平均二乗誤差(Mean Squared Error, MSE): モデルの性能の数値化のもうひとつの効果的な手法
最小化した誤差平方和コスト関数の平均値、線形回帰モデルの学習に使用
回帰モデルの比較や、グリッドサーチと考査検証によるパラメタのチューニングに役立つ$\displaystyle MSE = \frac{1}{n}\sum^n_{i=1}(y^{(i)}-\hat{y}^{(i)})^2$
In [26]:
from sklearn.metrics import r2_score
from sklearn.metrics import mean_squared_error
print('MSE train: %.3f, test: %.3f' % (
mean_squared_error(y_train, y_train_pred),
mean_squared_error(y_test, y_test_pred)))
print('R^2 train: %.3f, test: %.3f' % (
r2_score(y_train, y_train_pred),
r2_score(y_test, y_test_pred)))
- 訓練データセットに対してMSEは19.96, テストデータセットに対しては27.20 --- かなり大きいため、過学習の発生を示唆
- 場合によっては決定係数を求めるほうが効果的な場合あり
- 決定係数$R^2$: 標準化された平均二乗誤差とみなせる。モデルによって補足された応答の分散の割合を表す。
$\displaystyle R^2 = 1-\frac{SSE}{SST}$ ここでSSEは誤差平方和、SST(Sum of Squared Total)は: $\displaystyle SST = \sum^n_{i=1} (y^{(i)} - \mu_y)^2$
つまりSSTは応答分散
- $R^2$は平均二乗誤差(MSE)の尺度を変換した版にすぎない:
$\displaystyle R^2 = 1 - \frac{SSE}{SST} = 1 - \frac{\frac{1}{n}\sum^n_{i=1 }(y^{(i)}-\hat{y}^{(i)})^2}{\frac{1}{n}\sum^n_{i=1}(y^{(i)}-\mu_y)^2} = 1 - \frac{MSE}{Var(y)}$
訓練データセットでは$0 \leq R^2 \leq 1$、テストデータセットでは負の値になる可能性がある
$R^2=1$の場合、モデルは$MSE=0$となり、完全に適合
In [27]:
print('MSE train: %.3f, test: %.3f' % (
mean_squared_error(y_train, y_train_pred),
mean_squared_error(y_test, y_test_pred)))
print('R^2 train: %.3f, test: %.3f' % (
r2_score(y_train, y_train_pred),
r2_score(y_test, y_test_pred)))
- 訓練データを用いた評価ではモデルの$R^2=0.765$で悪くはないが、テストデータに対しては$0.673$と低い
10.6 回帰に正則化手法を使用¶
- 正則化(3章): 過学習に対処する手法の一つ、モデルの極端なパラメタの重みにペナルティを課す
- 線形回帰への正則化手法: リッジ回帰(Ridge regression), LASSO(Least Absolute Shrinkage and Selection Operator), Elastic Net
- リッジ回帰: L2ペナルティつきのモデル。$\displaystyle J(\boldsymbol{w})_{Ridge}=\sum^n_{i=1}(y^{(i)} - \hat{y}^{(i)})^2 - \lambda \| \boldsymbol{w}\|^2_2$
ここでL2 : $\displaystyle \lambda \|\boldsymbol{w}\|^2_2 = \lambda \sum^n_{j=1} w_j^2$
- LASSO: 教師ありの特徴選択手法としても役立つ
$\displaystyle J(\boldsymbol{w})_{LASSO} = \sum^n_{i=1} (y^{(i)} - \hat{y}^{(i)})^2 + \lambda \|\boldsymbol{w}\|_1$
ここで L1: $\lambda \|\boldsymbol{w}\|_1 = \lambda \sum^m_{j=1}|w_j|$
- LASSOには制約があり$m>n$の場合(サンプルよりも特徴量が多い場合)、最大値の$n$を選択する必要がある
- リッジ回帰とLASSOの折衷案が Elastic Net --- Elastic NetではL1とL2の両方を使用
$\displaystyle J(\boldsymbol{w})_{ElasticNet} = \sum^n_{i=1}(y^{(i)}-\hat{y}^{(i)})^2 + \lambda_1 \sum^m_{j=1}w_j^2 + \lambda_2\sum^m_{j=1}|w_j|$
- scikit-learnでの正則化回帰モデルの使用法
In [28]:
# リッジ回帰モデル---正則化の強さをalphaパラメタで指定
from sklearn.linear_model import Ridge
ridge = Ridge(alpha=1.0)
In [29]:
# LASSO回帰モデル---正則化の強さをalphaパラメタで指定
from sklearn.linear_model import Lasso
lasso = Lasso(alpha=1.0)
In [30]:
# Elastic Netモデル---l1_ratioによりL1ペナルティとL2ペナルティの比率を指定
from sklearn.linear_model import ElasticNet
elanet = ElasticNet(alpha=1.0, l1_ratio=0.5)
10.7 多項式回帰：線形回帰モデルから曲線を見出す¶
- ここまでは説明変数と目的変数の関係が線形であることを前提
- この前提が外れる場合、多項式の項が追加された多項式モデルを使用($d$は多項式の次数):
$ y = w_0 + w_1x + w_2 x^2 + \cdots + w_dx^d $ 線形回帰の係数$w$については線形なので、重回帰モデルとみなされる
10.7.1 scikit-learnを用いて多項式の項を追加¶
- PolynomialFeatures変換器クラスを用いて、説明変数が1つだけの単回帰問題に対し、2次の項($d-2$)を追加し、線形回帰と多項式回帰とを比較
In [18]:
# 1. 多項式の2次の項を追加
from sklearn.preprocessing import PolynomialFeatures
X = np.array([258.0, 270.0, 294.0,
320.0, 342.0, 368.0,
396.0, 446.0, 480.0, 586.0])\
[:, np.newaxis]
y = np.array([236.4, 234.4, 252.8,
298.6, 314.2, 342.2,
360.8, 368.0, 391.2,
390.8])
# 線形回帰モデルクラスのインスタンス
lr = LinearRegression()
pr = LinearRegression()
# 2次の多項式特徴量クラスのインスタンス
quadratic = PolynomialFeatures(degree=2)
# データに適合させデータを変換
X_quad = quadratic.fit_transform(X)
In [19]:
# 2. 単回帰モデルの学習
lr.fit(X, y)
X_fit = np.arange(250, 600, 10)[:, np.newaxis]
y_lin_fit = lr.predict(X_fit)
In [20]:
# 3. 多項式回帰のため、変換された特徴量で重回帰モデルを学習
pr.fit(X_quad, y)
y_quad_fit = pr.predict(quadratic.fit_transform(X_fit))
In [35]:
# 4. 結果の表示
plt.figure(figsize=(12,6))
# plot results
plt.scatter(X, y, label='training points')
plt.plot(X_fit, y_lin_fit, label='linear fit', linestyle='--')
plt.plot(X_fit, y_quad_fit, label='quadratic fit')
plt.legend(loc='upper left')
plt.tight_layout()
#plt.savefig('images/10_10.png', dpi=300)
plt.show()
- 目的変数と説明変数の関係は、多項式回帰のほうがよりよく捉えているように見える
In [37]:
y_lin_pred = lr.predict(X)
y_quad_pred = pr.predict(X_quad)
print('Training MSE linear: %.3f, quadratic: %.3f' % (
mean_squared_error(y, y_lin_pred),
mean_squared_error(y, y_quad_pred)))
print('Training R^2 linear: %.3f, quadratic: %.3f' % (
r2_score(y, y_lin_pred),
r2_score(y, y_quad_pred)))
- MSEは線形回帰が570に対し2次の多項式回帰は61、$R^2$は0.832に対して0.982なので、2次の多項式モデルの方がよりよく適合している
10.7.2 Housingデータセットで非線形関係をモデル化¶
- より複雑な問題に多項式の概念を適用---2次と3次の多項式を用いt、MEDV(住宅価格の中央値）とLSTAT(低所得者の割合)の関係をモデル化して、線形回帰モデルと多項式モデルを比較
In [41]:
X = df[['LSTAT']].values
y = df['MEDV'].values
regr = LinearRegression()
# create quadratic features
quadratic = PolynomialFeatures(degree=2)
cubic = PolynomialFeatures(degree=3)
X_quad = quadratic.fit_transform(X)
X_cubic = cubic.fit_transform(X)
# fit features
X_fit = np.arange(X.min(), X.max(), 1)[:, np.newaxis]
regr = regr.fit(X, y)
y_lin_fit = regr.predict(X_fit)
linear_r2 = r2_score(y, regr.predict(X))
regr = regr.fit(X_quad, y)
y_quad_fit = regr.predict(quadratic.fit_transform(X_fit))
quadratic_r2 = r2_score(y, regr.predict(X_quad))
regr = regr.fit(X_cubic, y)
y_cubic_fit = regr.predict(cubic.fit_transform(X_fit))
cubic_r2 = r2_score(y, regr.predict(X_cubic))
# plot results
plt.scatter(X, y, label='training points', color='lightgray')
plt.plot(X_fit, y_lin_fit,
label='linear (d=1), $R^2=%.2f$' % linear_r2,
color='blue',
lw=2,
linestyle=':')
plt.plot(X_fit, y_quad_fit,
label='quadratic (d=2), $R^2=%.2f$' % quadratic_r2,
color='red',
lw=2,
linestyle='-')
plt.plot(X_fit, y_cubic_fit,
label='cubic (d=3), $R^2=%.2f$' % cubic_r2,
color='green',
lw=2,
linestyle='--')
plt.xlabel('% lower status of the population [LSTAT]')
plt.ylabel('Price in $1000s [MEDV]')
plt.legend(loc='upper right')
#plt.savefig('images/10_11.png', dpi=300)
plt.show()
In [42]:
print('Training R^2 linear:%.3f quadratic: %.3f cubic: %.3f' % (linear_r2, quadratic_r2, cubic_r2))
- このグラフからわかるように、3次の多項式回帰モデルは他のモデルよりもよりよくMEDVとLSTATの関係を捉えている
- ただ、多項式の特徴量を追加すればするほどモデルの複雑さが増す$\rightarrow$過学習の可能性が高くなる
- モデルの汎化性能を推定できるよう、別のテストデータセットを用いてモデルの性能を評価することが勧められる
- 非線形関係のモデル化において、多項式の特徴量は必ずしも最良の選択とは限らないことに注意
今の例では、LSTATの対数とMEDVの平方根との線形回帰モデルのほうが適している
In [44]:
X = df[['LSTAT']].values
y = df['MEDV'].values
# transform features
X_log = np.log(X)
y_sqrt = np.sqrt(y)
# fit features
X_fit = np.arange(X_log.min()-1, X_log.max()+1, 1)[:, np.newaxis]
regr = regr.fit(X_log, y_sqrt)
y_lin_fit = regr.predict(X_fit)
linear_r2 = r2_score(y_sqrt, regr.predict(X_log))
# plot results
plt.scatter(X_log, y_sqrt, label='training points', color='lightgray')
plt.plot(X_fit, y_lin_fit,
label='linear (d=1), $R^2=%.2f$' % linear_r2,
color='blue',
lw=2)
plt.xlabel('log(% lower status of the population [LSTAT])')
plt.ylabel('$\sqrt{Price \; in \; \$1000s \; [MEDV]}$')
plt.legend(loc='lower left')
plt.tight_layout()
#plt.savefig('images/10_12.png', dpi=300)
plt.show()
In [45]:
print('Training R^2 linear:%.3f quadratic: %.3f cubic: %.3f' % (linear_r2, quadratic_r2, cubic_r2))
この結果から、この変換により他のモデルよりもよりよく捉えていることがわかる($R^2= 0.69$)
10.7.3 ランダムフォレストを用いて非線形関係を扱う¶
- ランダムフォレスト回帰: 複数の決定木からなるアンサンブル---区分線形関数の和とみなせる、つまり決定木アルゴリズムにより、入力空間をより『管理しやすい」小さな領域に分割
決定木回帰¶
- 決定木アルゴリズムの利点：非線形データを扱うのに、特徴量を変換する必要がない
- 二分決定木の情報利得(IG): $\displaystyle IG(D_p,x_i)=I(D_p) - \frac{N_{left}}{N_p} I(D_{left}) - \frac{N_{right}}{N_p}I(D_{right})$
ここで$x$は分割対象の特徴量、$N_p$は親ノードのサンプル個数、$I$は不純度関数、$D_p$は親ノードの学習サンプルのサブセット
- 回帰に決定木を使う---連続値の変数に適した不純度指標が必要: ノード$t$の不純度指標として平均二乗誤差(MSE)を定義
$\displaystyle I(t)=MSE(t) = \frac{1}{N_t}\sum_{i \in D_t} (y^{(i)}-\hat{y}_t)^2$
ただし$N$は$t$の学習サンブルの個数、$D_t$はノード$t$の学習サブセット（のインデックス集合）、$y^{(i)}$は真の目的値、$\hat{y}_t$は目的値の予測値 $\hat{y}_t = \frac{1}{N_t}\sum_{i \in D_t} y^{(i)}$
註: $I(t)$はしたがって目的変数の予測値と真の値の分散
- 決定木の学習直線の確認
In [46]:
from sklearn.tree import DecisionTreeRegressor
X = df[['LSTAT']].values
y = df['MEDV'].values
tree = DecisionTreeRegressor(max_depth=3)
tree.fit(X, y)
sort_idx = X.flatten().argsort()
lin_regplot(X[sort_idx], y[sort_idx], tree)
plt.xlabel('% lower status of the population [LSTAT]')
plt.ylabel('Price in $1000s [MEDV]')
#plt.savefig('images/10_13.png', dpi=300)
plt.show()
ランダムフォレスト回帰¶
- ランダムフォレストは外れ値に影響を受けにくく、パラメータのチューニングをあまり必要としない
- 一般に実験が必要なのは、アンサンブルの決定木の個数のみ
- 3章の分類用のアルゴリズムと回帰用のアルゴリズムの違いは、個々の決定木の成長に平均二乗誤差の条件を私用することだけ。予測される目的変数はすべての決定木の予測の平均
- Housingデータセットのすべての特徴量を使用し、サンブルの60％で学習、40％で性能評価
In [21]:
X = df.iloc[:, :-1].values
y = df['MEDV'].values
X_train, X_test, y_train, y_test = train_test_split(
X, y, test_size=0.4, random_state=1)
In [48]:
from sklearn.ensemble import RandomForestRegressor
forest = RandomForestRegressor(n_estimators=1000,
criterion='mse',
random_state=1,
n_jobs=-1)
forest.fit(X_train, y_train)
y_train_pred = forest.predict(X_train)
y_test_pred = forest.predict(X_test)
print('MSE train: %.3f, test: %.3f' % (
mean_squared_error(y_train, y_train_pred),
mean_squared_error(y_test, y_test_pred)))
print('R^2 train: %.3f, test: %.3f' % (
r2_score(y_train, y_train_pred),
r2_score(y_test, y_test_pred)))
- テストは訓練データに比べ、精度がかなり低い---過学習の傾向にある
- 目的変数と説明変数の関係は、うまく説明できている
- 予測値の残差
In [49]:
plt.scatter(y_train_pred,
y_train_pred - y_train,
c='steelblue',
edgecolor='white',
marker='o',
s=35,
alpha=0.9,
label='training data')
plt.scatter(y_test_pred,
y_test_pred - y_test,
c='limegreen',
edgecolor='white',
marker='s',
s=35,
alpha=0.9,
label='test data')
plt.xlabel('Predicted values')
plt.ylabel('Residuals')
plt.legend(loc='upper left')
plt.hlines(y=0, xmin=-10, xmax=50, lw=2, color='black')
plt.xlim([-10, 50])
plt.tight_layout()
# plt.savefig('images/10_14.png', dpi=300)
plt.show()
- 上図から明らかなようにテストデータよりも訓練データのほうに適合---外れ値の存在によって示される
- 残差は原点の周囲に完全にランダムに分布しているようにはみえない---モデルが情報を完全には把握できていないことを示唆
- ただし、10.5節の線形モデルよりも改善されている
- 残差プロットの非ランダム性への対処の普遍的な手法は存在しない:
変数変換、ハイパーパラメータのチューニング、単純・複雑なモデルの採用、外れ値の除去、変数の追加によるモデルの改良、などの手法
註: SVMは非線形回帰にも適用可能。参考文献: Support Vector Machines for Classification and Regression, S. R. Gunn and others, ISIS technical report, 14, 1998.
scikit-learnにSVM 回帰器がある: http://scikit-learn.org/stable/ modules/generated/sklearn.svm.SVR.html#sklearn.svm.SVR.