GitHubコードから読み解く！GradientBoostingClassifier@scikit-learnのfitアルゴリズム
はじめに
機械学習における基本タスクとして、2クラス分類があります。
データ属性を、上手く2分してくれる法則性を、機械に探索してもらうタスクです。
そして、探索してもらった法則性を用いることで、2クラス分類予測が、自動的に行えるようになります。
例えば、以下のイラストのようにです。
2クラス分類を実現してくれる機械学習手法は幾つか存在します。
中でも、人気のある手法が「Boosting」です。
特に、Boostingの改良版である「XGBoost」や「LightGBM」は、データ分析コンペサイト「Kaggle」でも、頻繁に用いられる手法であるようです。
理由は勿論、発揮精度が高い為です。
Boostingの考え方を、抽象的に把握
Boostingは複数の識別モデルを用いて、1つの判断を下すアルゴリズムです。
複数の識別モデルを用いて、1つの判断を下すアルゴリズムのことを、アンサンブルといいますが、Boostingもアンサンブル手法の1つです。
また、最も有名なアンサンブル手法は、Random Forestに代表されるような、Baggingになります。
BaggingとBoostingは似ていますが、個々の識別モデルの扱い方が、少し異なる形となっています。
Baggingは、独立した複数モデルの多数決で予測を行う方法です。
各々の識別モデルは、自分以外のモデルを配慮せず、別々に並列に生成されます。
すなわち、結果的に連携が生まれる形です。
一方、Boostingは、複数モデルの連携によって、1つの予測結果を出力する方法です。
各識別モデルは、自分以外のモデルに配慮しながら、相互関連しつつ直列的に生成されます。
連携することを前提に、各識別モデルが作られるのです。
具体的には、最初に作られた識別モデルの予測ミスを、2番目に作られた識別モデルが補填するような形で生成され、更に、そこで埋められなかった予測ミスを3番目の識別モデルが補填するような形で生成され、更に、そこで埋められなかった予測ミスを…という形で連携していきます。
繰り返していくと、いずれピッタリとFittingできる寸法です。
Boostingのアルゴリズムを、具体的に理解したい
さらに具体的なアルゴリズムの説明に、踏み込みたいと思います。
scikit-learnのプログラムをリバースして、具体的なアルゴリズムを追いかけてみます。
（プログラムのリバースは、より具体的な理解が捗るので、個人的に好きなエクササイズです。）
さて、上記の実装コードを読み進めていくと、コード内に「References」の記載があり、そこから以下論文に辿り着きました。
出典元を明らかにしてくれているのは有り難いですね。
また、この論文はかなり神がかったもので、これを読めばscikit-learnでのBoostingが概ね理解できる、というものでした。
ちなみに、一口にBoostingと言っても多種が存在する為、今回の記事で説明させていただく対象は、以下に絞りました。
① 導入理解向きな「Gaussian」解法によるGradient Boostingアルゴリズム
② sciki-learnのDefault実装である「Bernoulli」解法によるGradient Boostingアルゴリズム
という訳で、先ずは、①のGaussian解法を説明していきます。
GradientBoostingClassifier@scikit-learnアルゴリズムの「Gaussian」解法（①式の意味）
Gaussian解法によるGradient Boostingアルゴリズムは、以下式に基づいて、学習を行うものです。
Gaussian解法は、二乗誤差をコスト関数として、その最小化を目指すものです。
上記式でいうと、「Deviance」がコスト関数であり、目的変数「yᵢ」と、予測値「f(xᵢ)」との差の二乗の平均、つまり、分散を表しています。
この分散を最小化する「f(xᵢ)」探索がミッションとなります。
尚、「wᵢ」は一旦考えないで下さい。
これは、サンプルごとの重み付け値で、間違えて欲しくないサンプルには大きめの重みを、そうでないサンプルには小さめの重みを付与して、誤り傾向を制御するものです。
scikit-learnのオプションでいうと「sample_weight」です。
しかし、一般には「wᵢ」を割り当てないケースの方が多い為、全て1を代入することで、式から消してしまいましょう。
（※注意 : Σwᵢは、Nになります。）
「wᵢ」は、必要となったときに、再考すれば良いでしょう。
「yᵢ」は目的変数で、「xᵢ」は説明変数です。
「f(xᵢ)」は、「f(xᵢ) = f₀(xᵢ) + (η×f₁(xᵢ)) + (η×f₂(xᵢ))+…」という形で複数モデルの足し算で成り立っているものです。
この式での「η」は、モデル補填値の足し込み度合いで、scikit-learn上は「learning_rate」というパラメータになっています。
「Initial Value」の「f(x)」は、「f₀(xᵢ)」のことで、これは「Σ(yᵢ) / N」となります。
（※式に癖がありますが、scikit-learnをリバースすると、「yᵢ」の平均でした。2クラス分類の場合、ALL0.5決め打ちでも、そう問題は無いと思います。）
「zᵢ」は重要な項目で、「Gradient」になります。
「Gradient」式は、「Deviance」式を、「f(x)」について微分したものです。
「Deviance」が2乗の式なので、「Gradient」がZEROである時に、「Deviance」は最小値を取ります。
その為、「Deviance」の最小化は、「Gradient」をZEROに近付けることで実現できます。
GradientBoostingClassifier@scikit-learnアルゴリズムの「Gaussian」解法（②モデル生成Step）
式の意味が分かったところで、次にモデル生成Stepについて説明していきます。
結論から言うと、論文内の以下フローチャートになります。
これを日本語に置き換えると、以下となります。
（define）
i = 1〜N （サンプル数）（bias生成）
① 「Initial Value：f₀(xᵢ)」を、「f₀(xᵢ) = Σ(yᵢ)/N」という感じで求める。
（※「f₀(xᵢ)」はサンプル数次元のベクトルで、一律同じ値がセットされる）
② 「f(xᵢ)」に、「f₀(xᵢ)」を足しこむ。
（※ この時点では、f(xᵢ) = f₀(xᵢ) ）（1つ目の回帰木作成）
③ 「Gradient：zᵢ」を、「zᵢ = yᵢ - f(xᵢ)」という感じで求める。
（※「zᵢ」はサンプル数次元のベクトル）
④ 「zᵢ」を目的変数、「xᵢ」を説明変数として、回帰木「f₁(xᵢ)」を生成する。
⑤ 「f(xᵢ)」に、「η×f₁(xᵢ)」を足しこむ。
（※ f(xᵢ) = f₁(xᵢ) + (η × f₁(xᵢ)) ）（2つ目の回帰木作成）
⑥ 「Gradient：zᵢ」を、「zᵢ = yᵢ - f(xᵢ)」という感じで求める。
（※「zᵢ」はサンプル数次元のベクトル、③よりもΣzᵢが小さくなっているハズ）
⑦ 「zᵢ」を目的変数、「xᵢ」を説明変数として、回帰木「f₂(xᵢ)」を生成する。
⑧ 「f(xᵢ)」に、「η×f₁(xᵢ)」を足しこむ。
（※ f(xᵢ) = f₁(xᵢ) + (η × f₁(xᵢ)) + (η × f₂(xᵢ)) ）（3つ目の回帰木作成）
⑨ 「Gradient：zᵢ」を、「zᵢ = yᵢ - f(xᵢ)」という感じで求める。
（※「zᵢ」はサンプル数次元のベクトル、⑥よりもΣzᵢが小さくなっているハズ）
⑩ 「zᵢ」を目的変数、「xᵢ」を説明変数として、回帰木「f₃(xᵢ)」を生成する。
⑪ 「f(xᵢ)」に、「η×f₁(xᵢ)」を足しこむ。
（※ f(xᵢ) = f₁(xᵢ) + (η × f₁(xᵢ)) + (η × f₂(xᵢ)) + (η × f₃(xᵢ)) ）（※以降、指定数の回帰木を作り終えるまで、繰り返し…）
以上。
要するには、残差を予測し、補填する回帰木を、作っては足し込んでいく形です。
そうすることで、残差をどんどん減らしていく形となります。
GradientBoostingClassifier@scikit-learnアルゴリズムの「Gaussian」解法（③scikit-learnコードリバース）
先ほどの工程をscikit-learnがどう実装しているか、リバースさせていただこうと思います。
工程を、大きく以下3ブロックに分けて、説明していきます。
【Block ① ライブラリのインポート、データ準備】
【Block ② 学習コード】
【Block ③ テストコード】
Jupyter Notebookにて、1ブロック－1セルと対応させて実行すると、確認しやすいかもしれません。
【Block ① ライブラリのインポート、データ準備】
先ず、リバースコードを流すための準備です。
データは、scikit-learnに用意されている2クラス識別問題「breast cancer」を用います。
（※このブロックはリバースコードではありません）
# import lib
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import pydotplus
import sklearnfrom sklearn import datasets, metrics
from sklearn.model_selection import train_test_split
from sklearn.tree import DecisionTreeRegressor
from scipy.special import expit
from IPython.display import Image, display_png
from graphviz import Digraph# load sample data from sklearn
breast_cancer_data = datasets.load_breast_cancer()
columns = breast_cancer_data.feature_names# prepare X and y
X = breast_cancer_data.data
y = breast_cancer_data.target# print size of X and y
print('np.shape(X) = (%d, %d)' % np.shape(X))
print('np.shape(y) = (%d)' % np.shape(y))# separate train data and test data
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.5, random_state=0)print('np.shape(X_train) = (%d, %d)' % np.shape(X_train))
print('np.shape(y_train) = (%d)' % np.shape(y_train))
print('np.shape(X_test) = (%d, %d)' % np.shape(X_test))
print('np.shape(y_test) = (%d)' % np.shape(y_test))
以上で、説明変数「X」と目的変数「y」が準備できました。
また、必要なライブラリのimportも行いました。
ポイントとしては、回帰木クラスである「DecisionTreeRegressor」がimportされていることです。
これを用いて、scikit-learnのBoostingを再現します。
【Block ② 学習コード】
次に、scikit-learnを踏襲する形で、Gaussian解法によるGradient Boostingの、学習コードを記載します。
# ref : http://www.saedsayad.com/docs/gbm2.pdf# set param
gbm_n_estimators = 50
gbm_learning_rate = 0.8
tree_max_depth = 2
tree_criterion = 'mse'# option
tree_view = False # set true, if you want to watch tree structure# prepare work
TREE_LEAF = -1 # follow the specification of sklearn
estimators_ = np.empty((gbm_n_estimators, 1), dtype=np.object)
N = len(y_train)# ---------- [step 1] set initial value ----------
# fit initial model, init prediction (LogOddsEstimator)
pos = np.sum(y_train)
neg = y_train.shape[0] - pos# view number of pos and neg
print('pos = %d, neg = %d' % (pos, neg))# calc init prediction value
k = 1
prior = pos / (pos + neg)# prepare work variable for prediction value, and set
y_train_hat = np.empty((X_train.shape[0], k), dtype=np.float64)
y_train_hat.fill(prior)
# ---------------------------------------------# prepare work
residual_mean_stock = np.empty(gbm_n_estimators)
deviance_stock = np.empty(gbm_n_estimators)# set random seed
np.random.seed(0)# fit the boosting stages
for i in range(gbm_n_estimators): # ----- [step 2] calculation residual -----
# prediction by current predictor
residual = y_train - y_train_hat.ravel()
deviance = (1 / N) * np.sum(residual**2)
# -----------------------------------------
# stock
residual_mean_stock[i] = np.mean(residual)
deviance_stock[i] = deviance
print('[i = %3d] residual = %.3f, Deviance = %.3f' % (i, np.mean(np.abs(residual)), deviance)) # ----- [step 3] make tree for regression residual -----
# induce regression tree on residuals
tree = DecisionTreeRegressor(criterion=tree_criterion,
splitter='best',
max_depth=tree_max_depth,
random_state=i)
# fit next stage of trees
tree.fit(X_train, residual)
# ------------------------------------------------------
# if you want to visualize tree
if (tree_view):
dot_data = sklearn.tree.export_graphviz(tree,
out_file=None,
feature_names=columns,
label='all',
filled=True,
impurity=True,
node_ids=True,
proportion=True,
rotate=True,
rounded=True,
special_characters=True)
graph = pydotplus.graph_from_dot_data(dot_data)
display_png(Image(graph.create_png()))
plt.show() # compute leaf for each sample in X
terminal_regions = tree.apply(X_train) # add tree to ensemble
estimators_[i, 0] = tree
# update predictions (both in-bag and out-of-bag)
y_train_hat[:, 0] += (gbm_learning_rate * tree.tree_.value[:, 0, 0].take(terminal_regions, axis=0))# calc FPR, TPR and AUC
fpr_train, tpr_train, _ = metrics.roc_curve(y_train, y_train_hat)
auc_train = metrics.auc(fpr_train, tpr_train)# plot
fig = plt.figure(figsize=(5, 5), dpi=100)
plt.plot(fpr_train, tpr_train, alpha=0.5, linewidth=2, label=('ROC curve (area = %.2f)' % auc_train))
plt.scatter(fpr_train, tpr_train, alpha=0.5)
plt.legend(loc='best')
plt.title('ROC curve')
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.grid(True)
上記プログラムを実行すると、以下のROCカーブが出力されます。
上手く学習できてそうです。
【Block ③ テストコード】
最後に、作ったモデルのテストをするコードです。
学習に使用していないデータで評価します。
# ---------- [step 1] set initial value ----------
# prepare work variable for prediction value, and set
y_test_hat = np.empty((X_test.shape[0], k), dtype=np.float64)
y_test_hat.fill(prior)
# ---------------------------------------------# fit the boosting stages
for i in range(gbm_n_estimators): # ----- [step 2] calculation residual -----
# prediction by current predictor
residual = y_test - y_test_hat.ravel()
deviance = (1 / N) * np.sum(residual**2)
print('[i = %3d] residual = %.3f, Deviance = %.3f' % (i, np.mean(np.abs(residual)), deviance))
# -----------------------------------------
#
tree = estimators_[i, 0]
terminal_regions = tree.apply(X_test) # update predictions (both in-bag and out-of-bag)
y_test_hat[:, 0] += (gbm_learning_rate * tree.tree_.value[:, 0, 0].take(terminal_regions, axis=0))# calc FPR, TPR and AUC
fpr_test, tpr_test, _ = metrics.roc_curve(y_test, y_test_hat)
auc_test = metrics.auc(fpr_test, tpr_test)# plot
fig = plt.figure(figsize=(5, 5), dpi=100)
plt.plot(fpr_test, tpr_test, alpha=0.5, linewidth=2, label=('ROC curve (area = %.2f)' % auc_test))
plt.scatter(fpr_test, tpr_test, alpha=0.5)
plt.legend(loc='best')
plt.title('ROC curve')
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.grid(True)
上記コードを実施すると、以下のROCカーブが出力されます。
テストデータでも精度が出てそうです。
以上が、Gaussian解法によるGradient Boostingのリバースコードとなります。
最小二乗法が分かっている人であれば、割と理解しやすいかなと思います。
Gaussian解法よりも、少しトリッキーなBernoulli解法
次に、Bernoulli解法によるGradient Boostingの説明です。
Gaussian解法と比べて、若干理解が難しい手法となります。
ポイントは、2次微分値を用いる点です。
尚、この辺りの説明を書くに当たっては、以下の記事を参考にさせていただきました。
非常に神がかっている素晴らしい記事です。（有難や…😭）
GradientBoostingClassifier@scikit-learnアルゴリズムの「Bernoulli」解法（①式の意味）
「Bernoulli」解法は、以下式に基づいてBoostingを行います。
Gaussian解法に比べて、少し難しいですが、頑張って見ていきましょう。
先ずは、Bernoulli解法においても、「wᵢ」は一旦考えず、全て1に置き換えてしまいましょう。
それで式がスッキリします。
文字の意味はGaussian解法と同様で、「yᵢ」が目的変数、「xᵢ」が説明変数、「f(xᵢ)」は「f(xᵢ) = f₀(xᵢ) + (η×f₁(xᵢ)) + (η×f₂(xᵢ))+…」という形式の連結モデルです。
「Initial Value」の「f₀(xᵢ)」は、Gaussian解法同様に「yᵢ = {0, 1}」の比率によって求めます。
コスト関数の「Deviance」は、かなりトリッキーな形をしています。
どう最適化してもZEROにならなかったり、「yᵢ = {0, 1}」で加算されるコストが同一でなかったりします。
恐らく、微分式「Gradient」を導出したいがために、強引に「Deviance」式を形成したものと思われます。
つまり、「Deviance」自体に本質的な意味は無いと思われます。
大事なのは、微分式、及び、残差の式です。
（「そうじゃない！」という意見あれば、是非教えて下さい🙇♂️）
「pᵢ」は、所謂「sigmoid関数」です。
グラフにすると以下のような形。
目的変数「yᵢ」の値は「0 or 1」ですから、「pᵢ」の値がそれと合致する際には、「f(xᵢ)」の絶対値が10よりも大きな値を取っていそうです。
最適化ですが、Gaussian解法同様の「Gradient」式をZEROに近付けるアプローチを、Bernoulliでは行いません。
Bernoulli解法では、2次微分値による最適化「Newton-Raphson法」を用いて、直接「Deviance」を最小化させていきます。
それでは、Newton-Raphson法の適用方法を説明していきます。
先ず「Deviance」式を抽象的に、以下式表現とします。
この式を「f₀~𝑡-₁(xᵢ)」近辺にて、2次近似までのTaylor展開をします。
すると、以下式が導出されます。
ここで、「f𝑡(xᵢ)」は微妙なズレとして扱っています。
実際に、識別モデルを重ねるほどに「f𝑡(xᵢ)」の値は小さくなっていくので、識別モデル数が多くなることでTaylor展開の近似信頼度も上がっていく段取りです。
上記式ですが、「L'」を「Gᵢ」、「L''」を「Hᵢ」と置くと、以下のようにまとめられます。
この式に対して、以下変換を行います。
「const」は、 「f𝑡(xᵢ)」の値に左右されない項で、かつ、必ずZER0以上となる項をまとめたものです。
これらは、式全体を最小化する際、考慮から省いて問題ありません。
よって、この式から、「f𝑡(xᵢ)」は「Gᵢ / Hᵢ」の値を予測するような回帰木を生成すれば、「Deviance」の最小化を促せることになります。
尚、「Gᵢ / Hᵢ」とは「1次微分値 / 2次微分値」のことであり、具体的には、「Terminal node estimates」として書かれている式の導出値になります。
GradientBoostingClassifier@scikit-learnアルゴリズムの「Bernoulli」解法（②モデル生成Step）
次は、モデル生成Stepを説明していきます。
先ほど導出した式に基づいて、「Gᵢ / Hᵢ」を目的変数にした回帰木「f𝑡(xᵢ)」生成していく…と言いたいところですが、実はscikit-learnの実装は、そういうアルゴリズムになっていません。
なんと、「Gradient」式の「zᵢ」を目的変数に回帰木を作った後、葉の値を「Gᵢ / Hᵢ」で更新するということをしています。
論文中に記載されている識別モデル生成Step（フローチャート）にも、scikit-learnと同様の手順が書かれていました。
それが以下です。
尚、上記における、「subsampling rate」は、Gradient Boostingの基本アルゴリズムとは別の、過学習を防ぐための追加テクニックですので、今回の説明では取り上げないようにします。
上記のフローチャートを日本語で書き直すと、以下のようになります。
（define）
i = 1〜N （サンプル数）① 「Initial Value：f₀(xᵢ)」を、「f₀(xᵢ) = log(pos/neg)」という感じで求める。
（※「f₀(xᵢ)」はサンプル数次元のベクトルで、一律同じ値がセットされる）
② 「f(xᵢ)」に、「f₀(xᵢ)」を足しこむ。
（※ f(xᵢ) = f₀(xᵢ) ）③ 「Gradient：zᵢ」を、「zᵢ = yᵢ - sigmoid(f(xᵢ))」という感じで求める。
（※「zᵢ」はサンプル数次元のベクトル）
④ 「zᵢ」を目的変数、「xᵢ」を説明変数として、回帰木「f₁(xᵢ)」を生成する。
⑤ 回帰木の葉の値を、「Terminal node estimates」式の導出値で更新する。
（※葉に属するデータ群にて、Σを取る）
⑥ 「f(xᵢ)」に、「η × f₁(xᵢ)」を足しこむ。
（※ f(xᵢ) = f₀(xᵢ) + (η × f₁(xᵢ)) ）⑦ 「Gradient：zᵢ」を、「zᵢ = yᵢ - sigmoid(f(xᵢ))」という感じで求める。
（※「zᵢ」はサンプル数次元のベクトル）
⑧ 「zᵢ」を目的変数、「xᵢ」を説明変数として、回帰木「f₁(xᵢ)」を生成する。
⑨ 回帰木の葉の値を、「Terminal node estimates」式の導出値で更新する。
（※葉に属するデータ群にて、Σを取る）
⑩ 「f(xᵢ)」に、「η × f₂(xᵢ)」を足しこむ。
（※ f(xᵢ) = f₀(xᵢ) + (η × f₁(xᵢ)) + (η × f₂(xᵢ)) ）（※以降、指定数の回帰木を作り終えるまで、繰り返し…）
なぜ直接「Gᵢ / Hᵢ」を目的変数としないのかを、ソースコードをリバースしながら考察したところ、その場合にはZERO divide等の計算上不都合が挙げられそうでした。
サンプル別に「Gᵢ / Hᵢ」を計算すると、「Hᵢ」がZEROになりやすいのです。
試しに、その不都合に上手く対処しながら、「Gᵢ / Hᵢ」を目的変数にして回帰木でのBoostingを行ってみたのですが、収束が早いのは上記フローチャートの方式でした。
ここは非常に面白いポイントです。
GradientBoostingClassifier@scikit-learnアルゴリズムの「Bernoulli」解法（③scikit-learnコードリバース）
最後に、scikit-learnのリバースコードを記載します。
【Block ① ライブラリのインポート、データ準備】
先ず、リバースコードを流すための準備です。
データは、scikit-learnに用意されている2クラス識別問題「breast cancer」を用います。
（※このブロックはリバースコードではありません）
# import lib
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import pydotplus
import sklearnfrom sklearn import datasets, metrics
from sklearn.model_selection import train_test_split
from sklearn.tree import DecisionTreeRegressor
from scipy.special import expit
from IPython.display import Image, display_png
from graphviz import Digraph# load sample data from sklearn
breast_cancer_data = datasets.load_breast_cancer()
columns = breast_cancer_data.feature_names# prepare X and y
X = breast_cancer_data.data
y = breast_cancer_data.target# print size of X and y
print('np.shape(X) = (%d, %d)' % np.shape(X))
print('np.shape(y) = (%d)' % np.shape(y))# separate train data and test data
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.5, random_state=0)print('np.shape(X_train) = (%d, %d)' % np.shape(X_train))
print('np.shape(y_train) = (%d)' % np.shape(y_train))
print('np.shape(X_test) = (%d, %d)' % np.shape(X_test))
print('np.shape(y_test) = (%d)' % np.shape(y_test))
以上で、説明変数「X」と目的変数「y」が準備できました。
【Block ② 学習コード】
次に、scikit-learnからリバースしたBernoulli解法によるGradient Boosting学習コードを記載します。
# ref : http://www.saedsayad.com/docs/gbm2.pdf# set param
gbm_n_estimators = 50
gbm_learning_rate = 0.6
tree_max_depth = 2
tree_criterion = 'mse'# option
tree_view = False # set true, if you want to watch tree structure# prepare work
TREE_LEAF = -1 # follow the specification of sklearn
estimators_ = np.empty((gbm_n_estimators, 1), dtype=np.object)
N = len(y_train)# ---------- [step 1] set initial value ----------
# fit initial model, init prediction (LogOddsEstimator)
pos = np.sum(y_train)
neg = y_train.shape[0] - pos# view number of pos and neg
print('pos = %d, neg = %d' % (pos, neg))# calc init prediction value
k = 1
scale = 1.0
prior = scale * np.log(pos / neg) # if (pos == neg) then np.log(pos / neg) = 0# prepare work variable for prediction value, and set
y_train_hat = np.empty((X_train.shape[0], k), dtype=np.float64)
y_train_hat.fill(prior)
# ---------------------------------------------# prepare work
residual_mean_stock = np.empty(gbm_n_estimators)
deviance_stock = np.empty(gbm_n_estimators)# set random seed
np.random.seed(0)# fit the boosting stages
for i in range(gbm_n_estimators): # ----- [step 2] calculation residual -----
# prediction by current predictor
fx = expit(y_train_hat.ravel())
residual = y_train - fx
deviance = (-2 / N) * np.sum((y_train * fx) - np.log(1 + np.exp(fx)))
# -----------------------------------------
# stock
residual_mean_stock[i] = np.mean(residual)
deviance_stock[i] = deviance
print('[i = %3d] residual = %.3f, Deviance = %.3f' % (i, np.mean(np.abs(residual)), deviance)) # ----- [step 3] make tree for regression residual -----
# induce regression tree on residuals
tree = DecisionTreeRegressor(criterion=tree_criterion,
splitter='best',
max_depth=tree_max_depth,
random_state=i)
# fit next stage of trees
tree.fit(X_train, residual)
# ------------------------------------------------------
# if you want to visualize tree
if (tree_view):
dot_data = sklearn.tree.export_graphviz(tree,
out_file=None,
feature_names=columns,
label='all',
filled=True,
impurity=True,
node_ids=True,
proportion=True,
rotate=True,
rounded=True,
special_characters=True)
graph = pydotplus.graph_from_dot_data(dot_data)
display_png(Image(graph.create_png()))
plt.show() # ----- [step 4] update regression value -----
# update tree leaves (loss.update_terminal_regions)
# compute leaf for each sample in X
terminal_regions = tree.apply(X_train) # update each leaf (= perform line search)
for leaf in np.where(tree.tree_.children_left == TREE_LEAF)[0]: # Make a single Newton-Raphson step.
# our node estimate is given by:
# sum(w * (y - prob)) / sum(w * prob * (1 - prob))
# we take advantage that: y - prob = residual
terminal_region = np.where(terminal_regions == leaf)[0]
residual_ = residual.take(terminal_region, axis=0)
y_train_ = y_train.take(terminal_region, axis=0) numerator = np.sum(residual_)
denominator = np.sum((y_train_ - residual_) * (1 - y_train_ + residual_)) # prevents overflow and division by zero
if abs(denominator) < 1e-150:
tree.tree_.value[leaf, 0, 0] = 0.0
else:
# print('%.3f, %.3f, %.3f' % (tree.tree_.value[leaf, 0, 0], (numerator / denominator), (tree.tree_.value[leaf, 0, 0] - (numerator / denominator))))
tree.tree_.value[leaf, 0, 0] = numerator / denominator
# add tree to ensemble
estimators_[i, 0] = tree
# --------------------------------------------
# update predictions (both in-bag and out-of-bag)
y_train_hat[:, 0] += (gbm_learning_rate * tree.tree_.value[:, 0, 0].take(terminal_regions, axis=0))# calc FPR, TPR and AUC
fpr_train, tpr_train, _ = metrics.roc_curve(y_train, y_train_hat)
auc_train = metrics.auc(fpr_train, tpr_train)# plot
fig = plt.figure(figsize=(5, 5), dpi=100)
plt.plot(fpr_train, tpr_train, alpha=0.5, linewidth=2, label=('ROC curve (area = %.2f)' % auc_train))
plt.scatter(fpr_train, tpr_train, alpha=0.5)
plt.legend(loc='best')
plt.title('ROC curve')
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.grid(True)
上記プログラムを実行すると、以下のROCカーブが出力されます。
上手く学習が出来ていそうです。
【Block ③ テストコード】
最後に、作ったモデルのテストをするコードです。
学習に使用していないデータで評価します。
# prepare work variable for prediction value, and set
y_test_hat = np.empty((X_test.shape[0], k), dtype=np.float64)
y_test_hat.fill(prior)# fit the boosting stages
for i in range(gbm_n_estimators): #
fx = expit(y_test_hat.ravel())
residual = y_test - fx
deviance = (-2 / N) * np.sum((y_test * fx) - np.log(1 + np.exp(fx)))
print('[i = %3d] residual = %.3f, Deviance = %.3f' % (i, np.mean(np.abs(residual)), deviance))
#
tree = estimators_[i, 0]
terminal_regions = tree.apply(X_test)# update predictions (both in-bag and out-of-bag)
y_test_hat[:, 0] += (gbm_learning_rate * tree.tree_.value[:, 0, 0].take(terminal_regions, axis=0))# calc FPR, TPR and AUC
fpr_test, tpr_test, _ = metrics.roc_curve(y_test, y_test_hat)
auc_test = metrics.auc(fpr_test, tpr_test)# plot
fig = plt.figure(figsize=(5, 5), dpi=100)
plt.plot(fpr_test, tpr_test, alpha=0.5, linewidth=2, label=('ROC curve (area = %.2f)' % auc_test))
plt.scatter(fpr_test, tpr_test, alpha=0.5)
plt.legend(loc='best')
plt.title('ROC curve')
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.grid(True)
上記コードを実施すると、以下のROCカーブが出力されます。
テストデータでも精度が出ていそうです。
若干、Gaussian解法よりも、Bernoulli解法の方が高い精度が出ました。
さすがdefault採用されるだけある、と言えるかもしれません。
おわりに
という訳で、以上がGradientBoostingClassifier@scikit-learnの学習アルゴリズムの説明となります。
ここまで読んで下さった方、ありがとうございます。
やはり、リバースしてキッチリ確認するのは気持ちが良いですね。
次は、是非ともXGBoostのリバースをしてみたいと思います。
