機械学習（回帰問題）

#データの作成

from sklearn.model_selection import train_test_split

target_col = ''
exclude_cols = ['', '', '', '']
feature_cols = []
for col in dataset.columns:
    if col not in exclude_cols:
        feature_cols.append(col)
        
X = dataset[feature_cols]
y = dataset[target_col]           

X_train, X_test, y_train, y_test = train_test_split(X , y , test_size=0.3, random_state=1234)

print('X_test Features Shape: ', X_test.shape)
print('y_test Target Shape: ', y_test.shape)
print('X_train Features Shape: ', X_train.shape)
print('y_train Target Shape: ', y_train.shape)

回帰分析

線形回帰

誤差が最小になるように切片($\beta_0$)と傾き($\beta_1,\cdots,\beta_K$)を求める。

Y = \beta_0 + \beta_1 X_1 + \cdots + \beta_K X_K + \epsilon

　　　　　　　

from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error

lm = LinearRegression()
lm.fit(X_train, y_train)
y_pred = lm.predict(X_test)

lm_mse = mean_squared_error(y_test, y_pred)
print('LinerRegression RMSE: ', round(np.sqrt(lm_mse), 3))
print('決定係数:', lm.score(X_test, y_test))

#回帰係数
df_coef = pd.DataFrame({'feature':X.columns, 'coefficient':lm.coef_})
sns.barplot(x='coefficient', y='feature', 
            data=df_coef.sort_values('coefficient', ascending=False))
plt.show()

Lasso回帰

損失関数にL1正則化項を足し合わせたものを使用。
不要と判断される特徴量の係数$\beta_k$が0になるため、特徴量を削減することができる。

\begin{align}

　(損失関数) &= MSE + \alpha ||w||_1  \\
          &= \frac{1}{n} \displaystyle \sum_{ i = 1 }^{ n } (\hat{y_i} - (\beta_0+\cdots+\beta_i x_{i,k}+\cdots))^2 +\alpha \displaystyle \sum_{ k = 1 }^{ K } |\beta_k|\\

\end{align}

from sklearn.linear_model import LassoLars
from sklearn.metrics import mean_squared_error

lasso =  LassoLars(alpha=1.0)
lasso.fit(X_train, y_train)
y_pred = lasso.predict(X_test)

lasso_mse = mean_squared_error(y_test, y_pred)
print('LassoLars RMSE: ', round(np.sqrt(lasso_mse), 3))
print('決定係数:', lasso.score(X_test, y_test))

#回帰係数
df_coef = pd.DataFrame({'feature':X.columns, 'coefficient':lm.coef_})
sns.barplot(x='coefficient', y='feature', 
            data=df_coef.sort_values('coefficient', ascending=False))
plt.show()

#パラメーターの調整
from sklearn.model_selection import GridSearchCV

lasso = LassoLars()
params = {"alpha": [0.1*i for i in range(1,10,1)]}
gscv = GridSearchCV(lasso, 
                    param_grid=params,
                    verbose=1,  #詳細を表示
                    cv=3,  #交差検証の分割数
                    scoring='neg_mean_squared_error',
                    n_jobs=-1  #全てのプロセッサを使用する。
                    )
gscv.fit(X_train, y_train)
print('best params:', gscv.best_params_)
print('best MSE:', gscv.best_score_)

#best_pramsの可視化
gscv_df = pd.DataFrame(gscv.cv_results_['mean_test_score'],
                       index=np.round(params['alpha'],1))
plt.plot(gscv_df)
plt.xlabel('alpha')
plt.ylabel('mean_test_score')
plt.grid(True)
plt.show()

Ridge回帰

損失関数にL2正則化項を足し合わせたものを使用。
係数$\beta_k$は完全に0になることはない性質を持つ。

\begin{align}

　(損失関数) &= MSE + \alpha ||w||_2  \\
          &= \frac{1}{n} \displaystyle \sum_{ i = 1 }^{ n } (\hat{y_i} - (\beta_0+\cdots+\beta_i x_{i,k}+\cdots))^2 +\alpha \displaystyle \sum_{ k = 1 }^{ K } |\beta_k|^2\\

\end{align}

from sklearn.linear_model import Ridge
from sklearn.metrics import mean_squared_error

ridge = Ridge(alpha=1.0)
ridge.fit(X_train, y_train)
y_pred = ridge.predict(X_test)

ridge_mse = mean_squared_error(y_test, y_pred)
print('Ridge RMSE: ', round(np.sqrt(ridge_mse), 3))
print('決定係数:', ridge.score(X_test, y_test))

#回帰係数
df_coef = pd.DataFrame({'feature':X.columns, 'coefficient':ridge.coef_})
sns.barplot(x='coefficient', y='feature',
            data=df_coef.sort_values('coefficient', ascending=False))
plt.show()

#パラメーターの調整
from sklearn.model_selection import GridSearchCV

ridge = Ridge()
params = {"alpha": [0.1*i for i in range(1,10,1)]}
gscv = GridSearchCV(ridge, 
                    param_grid=params,
                    verbose=1,  #詳細を表示
                    cv=3,  #交差検証の分割数
                    scoring='neg_mean_squared_error',
                    n_jobs=-1  #全てのプロセッサを使用する。
                    )
gscv.fit(X_train, y_train)
print('best params:', gscv.best_params_)
print('best MSE:', gscv.best_score_)

#best_pramsの可視化
gscv_df = pd.DataFrame(gscv.cv_results_['mean_test_score'],
                       index=np.round(params['alpha'],1))
plt.plot(gscv_df)
plt.xlabel('alpha')
plt.ylabel('mean_test_score')
plt.grid(True)
plt.show()

決定木

回帰問題の場合、MSEが最小となるように分類し、各クラスの平均値を出力する。

重要度の求め方
ある特徴量で分岐させることでどのくらい値を下げられるかを示す。
回帰問題の場合はMSE、分類問題の場合はジニ係数を使用している。
（例）回帰問題

\begin{align}
　(OverallQualの重要度) &= \frac{I(OverallQual)_1+I(OverallQual)_2+I(OverallQual)_3}{\sum{I(X)}}\\
\\
I(OverallQual)_1 &= (699*6055931186.045)-((559*2297939697.554)+(110*7185735311.169))  \\
I(OverallQual)_2 &= (559*2297939697.554)-((405*1262133046.748)+(154*1982912431.17))  \\
I(OverallQual)_3 &= (110*7185735311.169)-((78*2987618758.878)+(32*9164978854.358))  \\

\end{align}

from sklearn.tree import DecisionTreeRegressor, plot_tree
from sklearn.metrics import mean_squared_error

dt = DecisionTreeRegressor(random_state=1234, 
                           max_depth=None  #全てのリーフがPureになるまで分岐する
                           )
dt.fit(X_train, y_train)
y_pred = dt.predict(X_test)

dt_mse = mean_squared_error(y_test, y_pred)
print('Decision Tree RMSE: ', round(np.sqrt(dt_mse), 3))
print('決定係数:', dt.score(X_test, y_test))

#特徴量の重要度
df_feature_importance = pd.DataFrame({'feature':X.columns, 'importance':dt.feature_importances_})
sns.barplot(x='importance', y='feature',
            data=df_feature_importance.sort_values('importance', ascending=False))
plt.show()

#決定木の可視化
plot_tree(dt, feature_names=X.columns, filled=True, rounded=True)
plt.show()

#パラメータ調整
from sklearn.model_selection import GridSearchCV

dt = DecisionTreeRegressor(random_state=1234)
params = {"max_depth": [i for i in range(4,10,1)]}
gscv = GridSearchCV(dt, param_grid=params, verbose=1, cv=3,scoring='neg_mean_squared_error', n_jobs=-1)
gscv.fit(X_train, y_train)
print('best params:', gscv.best_params_)
print('best MSE:', gscv.best_score_)

#グリッドサーチの可視化
gscv_df = pd.DataFrame(gscv.cv_results_['mean_test_score'],
                       index=params['max_depth'])
plt.plot(gscv_df)
plt.xlabel('max_depth')
plt.ylabel('mean_test_score')
plt.grid(True)
plt.show()

ランダムフォレスト

決定木をたくさん作り、多数決を行う。
回帰問題の場合は各決定木の平均値を出力する。
サンプリングはデータと特徴量両方で行われる。

重要度の求め方
OOBのデータを使用。
ある特徴量（列）をシャッフルし、ランダムフォレストの精度の差を求める。

(重要度)=(シャッフル前の精度)-(シャッフル後の精度)

精度が下がる → その特徴量が重要であった。
精度があまり変わらない → その特徴量があまり必要ではない。

from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import mean_squared_error

rf = RandomForestRegressor(random_state=1234,
                           n_estimators=100,
                           max_depth=None  #全てのリーフがPureになるまで分岐する
                           )
rf.fit(X_train, y_train)
y_pred = rf.predict(X_test)

rf_mse = mean_squared_error(y_test, y_pred)
print('Random Forest RMSE: ', round(np.sqrt(rf_mse), 3))
print('決定係数:', rf.score(X_test, y_test))

#特徴量の重要度
df_feature_importance = pd.DataFrame({'feature':X.columns, 'importance':rf.feature_importances_})
sns.barplot(x='importance', y='feature',
            data=df_feature_importance.sort_values('importance', ascending=False))
plt.show()

#ハイパーパラメータ調整
from sklearn.model_selection import GridSearchCV

rf = RandomForestRegressor(random_state=1234)
params = {"n_estimators":[10*i for i in range(8,12,1)], "max_depth":[5*i for i in range(1,5,1)]}
gscv = GridSearchCV(rf, 
                    param_grid=params,
                    verbose=1,  #詳細を表示
                    cv=3,  #交差検証の分割数
                    scoring='neg_mean_squared_error',
                    n_jobs=-1  #全てのプロセッサを使用する。
                    )
gscv.fit(X_train, y_train)
print('best params:', gscv.best_params_)
print('best MSE:', gscv.best_score_)

#グリッドサーチの可視化
gscv_df = pd.DataFrame(np.reshape(gscv.cv_results_['mean_test_score'],(len(params['n_estimators']),len(params['max_depth']))),
                       index=gscv.param_grid['n_estimators'],
                       columns=gscv.param_grid['max_depth'])
sns.heatmap(gscv_df, annot=True, fmt="1.2f")
plt.show()

XGBoost

勾配Boosting木の一種で精度が高く、使いやすい。

勾配Boosting木(Gradient Boosting Decision Tree)
〜特徴〜
・特徴量は数値　　　　　　　　　→　特徴量がある値より大きいか小さいかで分岐するため
・欠損値を扱うことができる　　　→　欠損値かどうかで分岐もすることができる
・変数間の相互作用が反映される　→　分岐の繰り返しによりできる
・スケーリングする必要がない　　→　値の大小関係のみが重要
・カテゴリー変数をone-hot-encodingする必要がない
　　　　　　　　　　　　　　　　→　決定木の分岐の際に、勝手に特徴量を抽出してくれる
・疎行列への対応がサポートされている

〜仕組み〜
目的変数$y_i$と予測値$\hat{y_i}$の差について学習し、逐次的に決定木を作成する。
ハイパーパラメータである決定木の本数分、上記を繰り返す。
最終的に各決定木のウェイト$W_m$の総和が予測値となる。

木の作成の仕方
Level-wise tree growth･･･深さ単位で木を成長させる。

重要度の求め方
gain → 特徴量によって得られた損失の減少を示し、基本的な手法。
他にも、"weight","cover","total_gain","total_cover"がある。

#sklearnAPI ver.
from xgboost import XGBRegressor
from sklearn.metrics import mean_squared_error

xgb =XGBRegressor(random_state=1234,
                  learning_rate=0.3,  # 学習率
                  max_depth=6,  #木の深さ（最も重要)
                  importance_type='gain'  #特徴量の算出方法
                  )
xgb.fit(X_train, y_train)
y_pred = xgb.predict(X_test)

xgb_mse = mean_squared_error(y_test, y_pred)
print('XGBoost RMSE: ', round(np.sqrt(xgb_mse), 3))
print('決定係数:', xgb.score(X_test, y_test))

#特徴量の重要度
df_feature_importance = pd.DataFrame({'feature':X.columns, 'importance':xgb.feature_importances_})
sns.barplot(x='importance', y='feature',
            data=df_feature_importance.sort_values('importance', ascending=False))
plt.show()

#ハイパーパラメータ調整
from sklearn.model_selection import GridSearchCV

xgb =XGBRegressor(random_state=1234)
params = {"learning_rate":[0.05*i for i in range(0,5)],
          "max_depth":[i for i in range(6,10,1)]}
gscv = GridSearchCV(xgb, 
                    param_grid=params,
                    verbose=1,  #詳細を表示
                    cv=3,  #交差検証の分割数
                    scoring='neg_mean_squared_error',
                    n_jobs=-1  #全てのプロセッサを使用する。
                    )
gscv.fit(X_train, y_train)
print('best params:', gscv.best_params_)
print('best MSE:', gscv.best_score_)

#グリッドサーチの可視化
gscv_df = pd.DataFrame(np.reshape(gscv.cv_results_['mean_test_score'],(len(params['learning_rate']),len(params['max_depth']))),
                       index=np.round(gscv.param_grid['learning_rate'], 2),
                       columns=gscv.param_grid['max_depth'])
sns.heatmap(gscv_df, annot=True, fmt="1.2f")
plt.show()

LGBM

XGBoostと比較して、高速かつ精度が同等程度と考えられている。

木の作成の仕方
Leaf-wise tree growth･･･葉単位で木を成長させる

重要度の求め方
split → 分岐に使われた回数
他にも、"gain"を使うことができる。

#sklearnAPI ver.
from lightgbm import LGBMRegressor
from sklearn.metrics import mean_squared_error

lgbm = LGBMRegressor(random_state=1234,
                     num_leaves=31,  #葉の数
                     min_child_samples =20,  #各ノードに入る最小データ数
                     importance_type='split'  #特徴量の算出方法
                     )
lgbm.fit(X_train, y_train)
y_pred = lgbm.predict(X_test)

lgbm_mse = mean_squared_error(y_test, y_pred)
print('LGBM RMSE: ', np.sqrt(lgbm_mse), 3)
print('決定係数:', lgbm.score(X_test, y_test))

#特徴量の重要度
df_feature_importance = pd.DataFrame({'feature':X.columns, 'importance':lgbm.feature_importances_})
sns.barplot(x='importance', y='feature',
            data=df_feature_importance.sort_values('importance', ascending=False))
plt.show()

#ハイパーパラメータ調整
from sklearn.model_selection import GridSearchCV

lgbm = LGBMRegressor(random_state=1234)
params = {"num_leaves":[i for i in range(28,34,1)],
          "min_child_samples":[i for i in range(17,23,1)]}
gscv = GridSearchCV(lgbm, 
                    param_grid=params,
                    verbose=1,  #詳細を表示
                    cv=3,  #交差検証の分割数
                    scoring='neg_mean_squared_error',
                    n_jobs=-1  #全てのプロセッサを使用する。
                    )
gscv.fit(X_train, y_train)
print('best params:', gscv.best_params_)
print('best MSE:', gscv.best_score_)

#グリッドサーチの可視化
gscv_df = pd.DataFrame(np.reshape(gscv.cv_results_['mean_test_score'],(len(params['num_leaves']),len(params['min_child_samples']))),
                       index=gscv.param_grid['num_leaves'],
                       columns=gscv.param_grid['min_child_samples'])
sns.heatmap(gscv_df, annot=True, fmt="1.2f")
plt.show()

###SVM

from sklearn.svm import SVR
from sklearn.metrics import mean_squared_error

svm = SVR(kernel='rbf', #linear:直線的に分割
          C=1,  #正規化パラメータ
          gamma='scale')  #決定境界の複雑さ scale:1/(n_features*X.var()), auto:1/n_features
svm.fit(X_train, y_train)
y_pred = svm.predict(X_test)
svm_mse = mean_squared_error(y_pred, y_test)

print('SVM RMSE: ', round(np.sqrt(svm_mse), 3))
print('決定係数:', svm.score(X_test, y_test))

#ハイパーパラメータ調整
from sklearn.model_selection import GridSearchCV

svm = SVR(kernel='rbf')
params = {'C':[0.1*i for i in range(0,10,1)], "gamma":[0.01*i for i in range(0,10,1)]}
gscv = GridSearchCV(svm, 
                    param_grid=params,
                    verbose=1,  #詳細を表示
                    cv=3,  #交差検証の分割数
                    scoring='neg_mean_squared_error',
                    n_jobs=-1  #全てのプロセッサを使用する。
                    )
gscv.fit(X_train, y_train)
print('best params:', gscv.best_params_)
print('best logloss:', gscv.best_score_)

#グリッドサーチの可視化
gscv_df = pd.DataFrame(np.reshape(gscv.cv_results_['mean_test_score'],(len(params['C']),len(params['gamma']))),
                       index=np.round(gscv.param_grid['C'],1),
                       columns=np.round(gscv.param_grid['gamma'],2))
sns.heatmap(gscv_df, annot=True, fmt="1.2f")
plt.show()

損失関数（最小二乗法）

以下の尤度関数が最小となるように$,f(x_i)$を求める。

\begin{align}
\displaystyle \sum_{ i = 1 }^{ n } (y_i - f(x_i))^2
\end{align}

精度評価指標

\begin{align}

MSE(MeanSquaredError) &= \frac{1}{n} \displaystyle \sum_{ i = 1 }^{ n } (\hat{y_i} - y_i)^2\\

RMSE(RootMeanSquaredError) &= \sqrt{\frac{1}{n} \displaystyle \sum_{ i = 1 }^{ n } (\hat{y_i} - y_i)^2}  \quad  \cdots  \quad 一般的 \\
MAE(Mean Absolute Error) &= \frac{1}{n} \displaystyle \sum_{ i = 1 }^{ n } |\hat{y_i} - y_i|\\

\end{align}

	メリット	デメリット
MSE	誤差に対して罰則が強い。 通常の線形回帰モデルを使う場合は、目的関数そのままなので直接最適化可能	元の単位の2乗なので解釈しにくい。
RMSE	元の単位に合っているので解釈しやすい。 通常の線形回帰モデルを使う場合は、目的関数そのままなので直接最適化可能	外れ値に対する罰則は強い。
MAE	外れ値の罰則を緩めることができる。	通常の線形回帰モデルでは直接最適化しているわけではない。

決定係数

モデルが予測できている割合。
分類問題の場合は、Accuracyと等価。

R^2 = 1 - \frac{\displaystyle \sum_{ i = 1 }^{ n } (\hat{y_i} - y_i)^2}
{\displaystyle \sum_{ i = 1 }^{ n } (\bar{y} - y_i)^2}\\

y_i \cdots 実測値 \quad \hat{y_i} \cdots 予測値 \quad \bar{y} \cdots 実測値の平均値

#あとがき
KaggleやSignateのコンペティションに熱中している間に、モデルの内容や指標の解釈を忘れてしまうので、すぐ概要を確認できるようにまとめました。
詳細を突き詰めると、間違えている箇所はあるかもしれませんがご容赦ください。

＜参考文献＞
・scikit-learnホームページ(https://scikit-learn.org/stable/)
・XGBoostホームページ(https://xgboost.readthedocs.io/en/latest/)
・LightGBMホームページ(https://lightgbm.readthedocs.io/en/latest/)
・Kaggleで勝つデータ分析の技術（技術評論社/著：門脇大輔，阪田隆司，保坂桂佑，平松雄司）