【Python】RANSACを用いてロバスト回帰モデルを実装｜scikit-learn・機械学習による線形回帰分析

# データ前処理
import numpy as np
import pandas as pd

# データ可視化
import matplotlib.pyplot as plt
import seaborn as sns
plt.style.use("ggplot")
# %matplotlib inline

# グラフの日本語表記対応
from matplotlib import rcParams
rcParams["font.family"]     = "sans-serif"
rcParams["font.sans-serif"] = "Hiragino Maru Gothic Pro"

# データセット読込
from sklearn.datasets import load_boston
boston = load_boston()

# DataFrame作成
df = pd.DataFrame(boston.data)
df.columns = boston.feature_names
df["MEDV"] = boston.target

from sklearn.linear_model import LinearRegression
from sklearn.linear_model import RANSACRegressor
from sklearn.model_selection import train_test_split

# 変数定義
X = df[['RM']].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)

# RANSACを用いたロバスト回帰モデル
model_ransac = RANSACRegressor(
                 base_estimator = LinearRegression(),  # ベースモデル
                 min_samples=50,                       # ランダムに選択されたデータサンプルの最小数
                 residual_threshold=8.0,               # 正常値として分類されるデータサンプルの最大残差（ここを大きくすると正常値として許容されるサンプル数が増える）
                 max_trials=200,                       # イテレーションの最大回数
                 stop_n_inliers=np.inf,                # 左記指定以上の正常値が見つかった場合は、イテレーション停止
                 stop_score=np.inf,                    # スコアが左記閾値よりも大きい場合の上位のイテレーション停止
                 stop_probability=0.99,                # トレーニングデータの少なくとも1つの外れ値のないセットがRANSACでサンプリングされると、RANSACのイテレーションは停止
                 loss='absolute_loss',                 # 絶対誤差(absolute_loss)または二乗誤差(squared_error)が選択可能
                 random_state=0                        # ランダムシード
                )

# モデル学習
model_ransac.fit(X_train, y_train)


inlier  = model_ransac.inlier_mask_    # 正常値の場合Trueを返す配列
outlier = np.logical_not(inlier)       # 外れ値の場合Trueを返す配列

x_lim       = np.arange(2, 11, 1)                         # 2~11までの配列を作成
pred_ransac = model_ransac.predict(x_lim[:, np.newaxis])  # 上記の配列を整数値に直し、モデルに適用

# 正常値プロット
plt.scatter(X_train[inlier],    # 横軸
            y_train[inlier],    # 縦軸
            c='blue',           # プロットカラー
            edgecolor="black",  # プロットカラー（枠線）
            s=20,               # プロットサイズ
            marker='o',         # プロット形状
            label='正常値')      # ラベル

# 外れ値プロット
plt.scatter(X_train[outlier],   # 横軸
            y_train[outlier],   # 縦軸
            c='green',          # プロットカラー
            edgecolor='black',  # プロットカラー（枠線）
            s=30,               # プロットサイズ
            marker='o',         # プロット形状
            label='外れ値'       # ラベル
           )

# 予測値プロット
plt.plot(x_lim, 
         pred_ransac, 
         color='red', 
         lw=3,
         label="回帰直線",
        )   

# グラフの書式設定
plt.xlim([3,10])
plt.ylim([y.min()-15,y.max()+10])
plt.xlabel('RM')
plt.ylabel('MEDV')
plt.legend(loc='upper left')
plt.show()

# 直線情報
print('傾き: %.2f' % model_ransac.estimator_.coef_[0])
# 出力結果（傾き）
# 11.42
print('切片: %.2f' % model_ransac.estimator_.intercept_)
# 出力結果（切片）
# -49.01

model_ransac = RANSACRegressor(
                    base_estimator=None,
                    min_samples=None,
                    residual_threshold=None,
                    max_trials=100, 
                    max_skips=inf, 
                    stop_n_inliers=inf, 
                    stop_score=inf, 
                    stop_probability=0.99,
                    loss='absolute_error',
                    random_state=None
                    )

from sklearn.metrics import r2_score            # 決定係数
from sklearn.metrics import mean_squared_error  # RMSE

# 予測値(Train）
y_train_pred = model_ransac.predict(X_train)
# 予測値（Test)
y_test_pred = model_ransac.predict(X_test)


# 平均平方二乗誤差(RMSE)
print('RMSE 学習: %.2f, テスト: %.2f' % (
        mean_squared_error(y_train, y_train_pred, squared=False), # 学習
        mean_squared_error(y_test, y_test_pred, squared=False)    # テスト
      ))

# 決定係数(R^2)
print('R^2 学習: %.2f, テスト: %.2f' % (
        r2_score(y_train, y_train_pred), # 学習
        r2_score(y_test, y_test_pred)    # テスト
      ))

# 出力結果
# RMSE 学習: 6.66, テスト: 7.12
# R^2  学習: 0.48, テスト: 0.39

# 予測値と残差をプロット（学習データ）
plt.scatter(y_train_pred,             # グラフのx値(予測値)  
            y_train_pred - y_train,   # グラフのy値(予測値と学習値の差)
            c='blue',                 # プロットの色
            marker='o',               # マーカーの種類
            s=40,                     # マーカーサイズ
            alpha=0.7,                # 透明度
            label='学習データ')         # ラベルの文字


# 予測値と残差をプロット（テストデータ）
plt.scatter(y_test_pred,            
            y_test_pred - y_test, 
            c='red',
            marker='o', 
            s=40,
            alpha=0.7,
            label='テストデータ')

# グラフの書式設定
plt.xlabel('予測値')
plt.ylabel('残差')
plt.legend(loc='upper left')
plt.hlines(y=0, xmin=-20, xmax=60, lw=2, color='black')
plt.xlim([-20, 60])
plt.ylim([-50, 40])
plt.tight_layout()
plt.show()

from sklearn.linear_model import LinearRegression
from sklearn.linear_model import RANSACRegressor
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)

# ロバスト回帰モデルのインスタンス
model_ransac = RANSACRegressor(
                 base_estimator = LinearRegression(),  # ベースモデル
                 min_samples=50,                       # ランダムに選択されたデータサンプルの最小数
                 residual_threshold=8.0,               # 正常値として分類されるデータサンプルの最大残差（ここを大きくすると正常値として許容されるサンプル数が増える）
                 max_trials=200,                       # イテレーションの最大回数
                 stop_n_inliers=np.inf,                # 少なくともこの数の正常値が見つかった場合は、イテレーション停止
                 stop_score=np.inf,                    # スコアが左記閾値よりも大きい場合の上位のイテレーション停止
                 stop_probability=0.99,                # トレーニングデータの少なくとも1つの外れ値のないセットがRANSACでサンプリングされると、RANSACのイテレーションは停止
                 loss='absolute_loss',                 # 絶対誤差(absolute_loss)または二乗誤差(squared_error)が選択可能
                 random_state=0                        # ランダムシード
                )

# モデル学習
model_ransac.fit(X_train, y_train)

from sklearn.metrics import r2_score            # 決定係数
from sklearn.metrics import mean_squared_error  # RMSE


# 予測値(Train）
y_train_pred = model_ransac.predict(X_train)
# 予測値（Test)
y_test_pred = model_ransac.predict(X_test)

# 平均平方二乗誤差(RMSE)
print('RMSE 学習: %.2f, テスト: %.2f' % (
        mean_squared_error(y_train, y_train_pred, squared=False), # 学習
        mean_squared_error(y_test, y_test_pred, squared=False)    # テスト
      ))

# 決定係数(R^2)
print('R^2 学習: %.2f, テスト: %.2f' % (
        r2_score(y_train, y_train_pred), # 学習
        r2_score(y_test, y_test_pred)    # テスト
      ))


# 出力結果
# RMSE 学習: 4.90, テスト: 5.95
# R^2  学習: 0.72, テスト: 0.57

# 予測値と残差をプロット（学習データ）
plt.scatter(y_train_pred,             # グラフのx値(予測値)  
            y_train_pred - y_train,   # グラフのy値(予測値と学習値の差)
            c='blue',                 # プロットの色
            marker='o',               # マーカーの種類
            s=40,                     # マーカーサイズ
            alpha=0.7,                # 透明度
            label='学習データ')         # ラベルの文字


# 予測値と残差をプロット（テストデータ）
plt.scatter(y_test_pred,            
            y_test_pred - y_test, 
            c='red',
            marker='o', 
            s=40,
            alpha=0.7,
            label='テストデータ')

# グラフの書式設定
plt.xlabel('予測値')
plt.ylabel('残差')
plt.legend(loc='upper left')
plt.hlines(y=0, xmin=-20, xmax=60, lw=2, color='black')
plt.xlim([-20, 60])
plt.ylim([-50, 40])
plt.tight_layout()
plt.show()