Ridge回帰の実装

概要

実装はsklearnを利用する。

環境

sklearn 0.21.2
matplotlib 3.1.1

実装

データ

ボストンデータを利用する。

from sklearn.datasets import load_boston
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler

data = load_boston()
x_train, x_test, y_train, y_test = train_test_split(data['data'], data['target'], test_size=0.2, random_state=0)
scaler = StandardScaler()
scaler.fit(data['data'])
x_train = scaler.transform(x_train)
x_test = scaler.transform(x_test)

モデル

from sklearn.linear_model import Ridge

model = Ridge()
model

デフォルトのパラメータ

Ridge(alpha=1.0, copy_X=True, fit_intercept=True, max_iter=None,
      normalize=False, random_state=None, solver='auto', tol=0.001)

ハイパーパラメータチューニング

今回のハイパーパラメータはalphaのみ。グリットサーチクロスバリデーションで最適化する。5-fold cross validationをする。

#grid_search cross validation
n_splits = 5
params = {'alpha' : [1.0]}

グリッドサーチクロスバリデーション

グリッドサーチクロスバリデーションを定義する。訓練データでも評価したいのでGridSearchCVを使わずに自分で定義する。

from sklearn.model_selection import train_test_split, KFold, ParameterGrid
from statistics import mean, stdev
import copy

def gscv(datas, grid_params, n_splits=5):
    x_train, y_train = datas
    kf = KFold(n_splits=n_splits, shuffle=False, random_state=0)
    results = {
        'params' : [],
        'validation_score_li' : [],
        'validation_score_mean' : [],
        'validation_score_std' : [],
        'train_score' : []
    }
    for params in ParameterGrid(grid_params):
        results['params'].append(params)
        for k, v in params.items():
            results.setdefault(f'param_{k}', []).append(v)
        score_li = []
        model.set_params(**params)
        for i_train, i_test in kf.split(x_train):
            x_cv_train = x_train[i_train]
            y_cv_train = y_train[i_train]
            x_cv_test = x_train[i_test]
            y_cv_test = y_train[i_test]
            model.fit(x_cv_train, y_cv_train)
            score = model.score(x_cv_test, y_cv_test)
            score_li.append(score)
        results['validation_score_li'].append(score_li)
        results['validation_score_mean'].append(mean(score_li))
        results['validation_score_std'].append(stdev(score_li))
        model.fit(x_train, y_train)
        score = model.score(x_train, y_train)
        results['train_score'].append(score)
    return results

次に、validationスコアの平均と訓練データでのスコアをプロットする関数を定義する。

import matplotlib.pyplot as plt

def ShowLiner(results, x_title):
    x_title = f'param_{x_title}'
    y_title = 'score'
    x_ticks = results[x_title]
    x = range(0, len(x_ticks))
    y_cv = results['validation_score_mean']
    y_train = results['train_score']
    fig = plt.figure(figsize=(10, 6))
    ax = fig.add_subplot(111)
    ax.plot(x, y_cv, label='cv_score')
    ax.plot(x, y_train, label='train_score')
    ax.set_xlabel(x_title)
    ax.set_ylabel(y_title)
    ax.set_xticks(x)
    ax.set_xticklabels(x_ticks)
    plt.grid()
    plt.legend()
    plt.show()

最後に、ハイパーパラメータをチューニングするパイプラインを関数で定義する。ハイパーパラメータをチューニングし、グラフを表示し、もっともよかったパラメータを自動的に更新する。

def search_param1(param_name, param_values, params):
    new_params = params.copy()
    new_params[param_name] = param_values
    print(f'set new_params: {new_params}')
    results = gscv([x_train, y_train], new_params)
    cv_score_li = results['validation_score_mean']
    best_index = cv_score_li.index(max(cv_score_li))
    print(f"cv score: {results['validation_score_mean'][best_index]:.2f}+-{results['validation_score_std'][best_index]:.3f}")
    print(f"best params: {results['params'][best_index]}")
    ShowLiner(results, param_name)
    params[param_name] = [results['params'][best_index][param_name]]

実行する

search_param1('alpha', [0, 1, 10, 100, 1000], params)

set new_params: {'alpha': [0, 1, 10, 100, 1000]}
cv score: 0.50+-0.211
best params: {'alpha': 10}

青が検証データ、黄色が訓練データである。alphaが10のときに一番スコアがよかった。

訓練と予測

#train
model.set_params(**params)
model.fit(x_train, y_train)

#predict
y_pred = model.predict(x_test)

結果と分析

決定係数は組み込み関数で計算できる。

r2 = model.score(x_test, y_test)
print(r2)

0.5796111714164923

予測した値をプロットしてみる。

import matplotlib.pyplot as plt

fig = plt.figure(figsize=(6,6))
ax = fig.add_subplot(111)
ax.scatter(y_test, y_pred)
ax.set_xlabel('y_test')
ax.set_ylabel('y_pred')
plt.show()

切片と係数を取得し、係数の大きさを確認してみる。

w_vec = model.coef_
b = model.intercept_

import matplotlib.pyplot as plt
import numpy as np

w_vec = model.coef_
names = data['feature_names']
mask = np.where(w_vec > 0, True, False)

w_vec = np.abs(w_vec)
order = np.argsort(w_vec)
mask = mask[order]

fig = plt.figure(figsize=(6, 10))
ax = fig.add_subplot(111)
line_n = ax.barh(names[order], w_vec[order], color='orange')
line_p = ax.barh(names[order][mask], w_vec[order][mask], color='blue')
line_p.set_label('Positive number')
line_n.set_label('Negative number')
ax.legend()
plt.show()

参考文献

https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.Ridge.html