【機器學習筆記】——Bagging、Boosting、Stacking（RF / Adaboost / Boosting Tree / GBM / GBDT / XGBoost / LightGBM）

時間 2021-01-02

標籤集成學習 Bagging Boosting GBDT XGBoost 欄目 C&C++ 简体版

原文原文鏈接

集成學習

概念

首先，讓我們先來了解一下，什麼是集成學習（Ensemble Learning，寫作銀桑，讀作昂桑）。集成學習由訓練數據構建一組基分類器，然後通過對每個基分類器的預測進行投票來進行分類，對每個個體學習器的預測進行（加權）平均來進行迴歸，有時也被稱爲多分類器系統、基於委員會的學習等。如果把單個分類器比作一個決策者的話，集成學習的方法就相當於多個決策者共同進行一項決策。嚴格來說，集成學習並不算是一種分類器，而是一種分類器結合的方法。

集成學習將多個分類方法聚集在一起，以提高分類的準確率。通常一個集成分類器的分類性能會好於單個分類器

對於二分類問題，假設基分類器的錯誤率爲 $\epsilon$ ：
$P(h_i(\mathcal{x} \ne \mathcal{y})) = \epsilon$

通過投票我們得到集成學習的分類預測：
$H(\mathcal{x}) = sign\left(\sum_{i = 1}^{T} h_i(\mathcal{x})\right)$

假設分類器的錯誤率相互獨立，則由Hoeffding不等式有：
$P(H(\mathcal{x} \ne \mathcal{y})) = \sum_{k = 0}^{\lfloor T/2 \rfloor} \dbinom{T}{k}(1 - \epsilon)^k \epsilon^{T - k} \le exp\left(-\frac{1}{2}T(1-2\epsilon)^2 \right)$

即隨着集成中基分類器數目 $T$ 增大，集成的錯誤率會指數級下降，最終趨於0。當然實際中個體學習器錯誤率獨立的假設是不成立的

集成學習的很多理論研究都是針對弱分類器（精度略高於50%）進行的，但在實踐中往往會使用較強的分類器
個體學習器要有一定的「準確性」，並且要有「多樣性」。一般的，準確性很高之後，要增加多樣性就需犧牲準確性。如何產生並結合「好而不同」的個體學習器，恰是集成學習研究的核心

下面這個例子，圖(a)中每個學習器只有67%的精度，但集成學習卻達到了100%；圖(b)中三個學習器完全一樣，集成之後性能沒有提高；圖©中每個學習器只有33%的精度集成學習的效果更差。

思維導圖

Bagging算法

概念

Bagging算法（Bootstrap AGGregatING）是一種並行式集成算法。它基於自助採樣法（Bootstrap sampling）採樣出T個含有m個訓練樣本的採樣集，然後基於每個採樣集訓練出一個基學習器，再將這些基學習器進行結合。

算法的僞代碼如下：

從偏差-方差分解的角度看，降低一個估計的方差的方式是把多個估計平均起來，所以Bagging主要關注降低方差，因此它在不剪枝決策樹、神經網絡等易受樣本擾動的學習器上效用更爲明顯。同時，Bagging並不能減小模型的偏差，所以應儘量選擇偏差較小的基分類器，如未剪枝的決策樹。

編程（分類）

參考：https://github.com/vsmolyakov/experiments_with_python/blob/master/chp01/ensemble_methods.ipynb

我們用iris數據集來看一下不同基分類器Bagging的效果

import itertools
import numpy as np

import seaborn as sns
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec

from sklearn import datasets

from sklearn.tree import DecisionTreeClassifier
from sklearn.neighbors import KNeighborsClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import RandomForestClassifier

from sklearn.ensemble import BaggingClassifier
from sklearn.model_selection import cross_val_score, train_test_split

from mlxtend.plotting import plot_learning_curves
from mlxtend.plotting import plot_decision_regions

np.random.seed(0)

導入數據

iris = datasets.load_iris()
X, y = iris.data[:, 0:2], iris.target

X.shape, y.shape

((150, 2), (150,))

設置基分類器和集成方法

# 配置基分類器的參數
clf1 = DecisionTreeClassifier(criterion='entropy', max_depth=1)    # 深度爲1的決策樹
clf2 = KNeighborsClassifier(n_neighbors=1)    # k=1的k近鄰

# 每個Bagging有10個基分類器
bagging1 = BaggingClassifier(base_estimator=clf1, n_estimators=10, max_samples=0.8, max_features=0.8)
bagging2 = BaggingClassifier(base_estimator=clf2, n_estimators=10, max_samples=0.8, max_features=0.8)

查看單個學習器和Bagging後的效果

label = ['Decision Tree(depth = 1)', 'K-NN(k = 1)', 'Bagging Tree', 'Bagging K-NN']
clf_list = [clf1, clf2, bagging1, bagging2]

fig = plt.figure(figsize=(10, 8))
gs = gridspec.GridSpec(2, 2)    # 創建 2x2 的畫布
grid = itertools.product([0,1],repeat=2)    # 返回參數的笛卡爾積的元組，每個元素可以用 2 次

for clf, label, grd in zip(clf_list, label, grid):        
    scores = cross_val_score(clf, X, y, cv=3, scoring='accuracy')    # Array of scores of the estimator for each run of the cross validation.
    print("Accuracy: %.2f (+/- %.2f) [%s]" %(scores.mean(), scores.std(), label))
        
    clf.fit(X, y)
    ax = plt.subplot(gs[grd[0], grd[1]])
    fig = plot_decision_regions(X=X, y=y, clf=clf, legend=2)    # 繪製決策邊界
    plt.title(label)

plt.show()

Accuracy: 0.63 (+/- 0.02) [Decision Tree(depth = 1)]
Accuracy: 0.70 (+/- 0.02) [K-NN(k = 1)]
Accuracy: 0.66 (+/- 0.02) [Bagging Tree]
Accuracy: 0.61 (+/- 0.02) [Bagging K-NN]

可以發現決策樹的集成有一定效果，而k近鄰的效果反而變差了。這是因爲k近鄰是穩定學習器，不易受樣本擾動影響。

再來看看數據集的劃分對訓練結果的影響

#plot learning curves
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)
    
plt.figure()
plot_learning_curves(X_train, y_train, X_test, y_test, bagging1, print_model=False, style='ggplot')
plt.show()

雖然每次結果有較大的差異，但是我們還是能得到一些結論：對於給定的劃分比例，模型在訓練集上的錯分率基本在30%上下，而在測試集上，當訓練集較小時，模型在測試集上表現很差，當訓練集比例逐漸增加時，模型在測試集上的表現慢慢變好，一般在訓練集佔到80%左右時，模型在訓練集和測試集上都能取得較好的表現

接下來我們看一下學習器的個數對bagging效果的影響

#Ensemble Size
num_est = np.concatenate(([1],np.arange(5,105,5)))
bg_clf_cv_mean = []
bg_clf_cv_std = []

for n_est in num_est:    
    bg_clf = BaggingClassifier(base_estimator=clf1, n_estimators=n_est, max_samples=0.8, max_features=0.8)
    scores = cross_val_score(bg_clf, X, y, cv=3, scoring='accuracy')
    bg_clf_cv_mean.append(scores.mean())
    bg_clf_cv_std.append(scores.std())
    

plt.figure()
# 繪製誤差棒圖，橫軸爲num_est，縱軸爲scores.mean()，誤差爲scores.std()
(_, caps, _) = plt.errorbar(num_est, bg_clf_cv_mean, yerr=bg_clf_cv_std, c='blue', fmt='-o', capsize=5)
for cap in caps:
    cap.set_markeredgewidth(1)                                                           
plt.ylabel('Accuracy'); plt.xlabel('Ensemble Size'); plt.title('Bagging Tree Ensemble');
plt.show()

基於交叉驗證的結果，我們可以看到在基分類器數量達到20後，集成的效果不再有明顯的提升

隨機森林

隨機森林（Random Forest, RF）可以看成是改進的Bagging算法，它可以更好地解決決策樹的過擬合問題。隨機森林以CART樹作爲基分類器，對比決策樹的Bagging，它有一個明顯的特點——隨機選擇特徵，對於普通的決策樹，我們會在每個節點在當前的所有特徵（假設有n個）中尋找最優特徵，而隨機森林的決策樹會在每個節點在隨機選取的特徵集合（推薦 $n_{rf} = \log_2 n$ ）中尋找最優特徵。這樣做的好處是如果訓練集中的幾個特徵對輸出的結果有很強的預測性，那麼這些特徵會被每個決策樹所應用，這樣會導致樹之間具有相關性，這樣並不會減小模型的方差，但同時隨機森林會不可避免的增加偏差。

RF的主要優點有：

訓練可以高度並行化，對於大數據時代的大樣本訓練速度有優勢。個人覺得這是的最主要的優點。
由於可以隨機選擇決策樹節點劃分特徵，這樣在樣本特徵維度很高的時候，仍然能高效的訓練模型。
在訓練後，可以給出各個特徵對於輸出的重要性
由於採用了隨機採樣，訓練出的模型的方差小，泛化能力強。
相對於Boosting系列的Adaboost和GBDT， RF實現比較簡單。
對部分特徵缺失不敏感。

RF的主要缺點有：

在某些噪音比較大的樣本集上，RF模型容易陷入過擬合。
取值劃分比較多的特徵容易對RF的決策產生更大的影響，從而影響擬合的模型的效果。

擴展

RF還有很多變種，其目的都是爲了進一步降低方差。

Extremely randomized Trees

Extremely randomized Trees（以下簡稱Extra Trees）放棄了Boostrap sampling的方式，而選擇使用原始數據集，而且在特徵選擇時更加暴力的隨機選擇一個特徵直接作爲最優特徵

*Totally Random Trees Embedding

參考：http://www.cnblogs.com/pinard/p/6156009.html

Totally Random Trees Embedding(以下簡稱 TRTE)是一種非監督學習的數據轉化方法。它將低維的數據集映射到高維，從而讓映射到高維的數據更好的運用於分類迴歸模型。我們知道，在支持向量機中運用了核方法來將低維的數據集映射到高維，此處TRTE提供了另外一種方法。

TRTE在數據轉化的過程也使用了類似於RF的方法，建立T個決策樹來擬合數據。當決策樹建立完畢以後，數據集裏的每個數據在T個決策樹中葉子節點的位置也定下來了。比如我們有3顆決策樹，每個決策樹有5個葉子節點，某個數據特徵xx劃分到第一個決策樹的第2個葉子節點，第二個決策樹的第3個葉子節點，第三個決策樹的第5個葉子節點。則x映射後的特徵編碼爲(0,1,0,0,0, 0,0,1,0,0, 0,0,0,0,1), 有15維的高維特徵。這裏特徵維度之間加上空格是爲了強調三顆決策樹各自的子編碼。

映射到高維特徵後，可以繼續使用監督學習的各種分類迴歸算法了。

*Isolation Forest

Isolation Forest（以下簡稱IForest）是一種異常點檢測的方法。它也使用了類似於RF的方法來檢測異常點。

對於在T個決策樹的樣本集，IForest也會對訓練集進行隨機採樣,但是採樣個數不需要和RF一樣，對於RF，需要採樣到採樣集樣本個數等於訓練集個數。但是IForest不需要採樣這麼多，一般來說，採樣個數要遠遠小於訓練集個數？爲什麼呢？因爲我們的目的是異常點檢測，只需要部分的樣本我們一般就可以將異常點區別出來了。

對於每一個決策樹的建立， IForest採用隨機選擇一個劃分特徵，對劃分特徵隨機選擇一個劃分閾值。這點也和RF不同。

另外，IForest一般會選擇一個比較小的最大決策樹深度max_depth,原因同樣本採集，用少量的異常點檢測一般不需要這麼大規模的決策樹。

對於異常點的判斷，則是將測試樣本點 $x$ 擬合到 $T$ 顆決策樹。計算在每顆決策樹上該樣本的葉子節點的深度 $h_t(x)$ ，從而可以計算出平均高度 $h(x)$ 。此時我們用下面的公式計算樣本點xx的異常概率:

$s(x,m)=2^{−\frac{h(x)}{c(m)}}$

其中，m爲樣本個數。 $c(m)$ 的表達式爲：

$c(m)=2\ln (m−1)+\xi−2\frac{m−1}{m}$

$\xi$ 爲歐拉常數， $s(x,m)$ 的取值範圍是 $[0,1]$ ，取值越接近於1，則是異常點的概率也越大。

編程（分類）

import time

# DecisionTree
clf3 = DecisionTreeClassifier(criterion='entropy', max_depth=2)

# Tree Bagging
tree_bagging = BaggingClassifier(base_estimator=clf3, n_estimators=100)

# Random Forest
random_forest = RandomForestClassifier(n_estimators=100, criterion='entropy', max_depth = 2)

label = ['Decision Tree', 'Tree Bagging', 'Random Forest']
clf_list = [clf3, tree_bagging, random_forest]

fig = plt.figure(figsize=(15, 4))
gs = gridspec.GridSpec(1, 3)
grid = itertools.product([0], [0, 1, 2])

for clf, label, grd in zip(clf_list, label, grid):
    b = time.time()
    scores = cross_val_score(clf, X, y, cv=3, scoring='accuracy')
    runtime = time.time() - b
    print("Accuracy: %.2f (+/- %.2f) [%s] Runtime: %.2f s" %(scores.mean(), scores.std(), label, runtime))
        
    clf.fit(X, y)
    ax = plt.subplot(gs[grd[0], grd[1]])
    fig = plot_decision_regions(X=X, y=y, clf=clf, legend=2)
    plt.title(label)

plt.show()

Accuracy: 0.63 (+/- 0.04) [Decision Tree] Runtime: 0.01 s
Accuracy: 0.77 (+/- 0.03) [Tree Bagging] Runtime: 0.40 s
Accuracy: 0.76 (+/- 0.04) [Random Forest] Runtime: 0.36 s

從結果可以看到隨機森林的運行速度比Tree Bagging要快一點，這是因爲在選擇最優特徵時，隨機森林並沒有從全部特徵中進行選擇

rf_label = ['DecisionTree','1 Trees','5 Trees','10 Trees','20 Trees','100 Trees']

random_forest1 = RandomForestClassifier(n_estimators=1, criterion='entropy', max_depth = 2)
random_forest5 = RandomForestClassifier(n_estimators=5, criterion='entropy', max_depth = 2)
random_forest10 = RandomForestClassifier(n_estimators=10, criterion='entropy', max_depth = 2)
random_forest20 = RandomForestClassifier(n_estimators=20, criterion='entropy', max_depth = 2)
random_forest100 = RandomForestClassifier(n_estimators=100, criterion='entropy', max_depth = 2)

clf_list = [clf3, random_forest1, random_forest5, random_forest10, random_forest20, random_forest100]

fig = plt.figure(figsize=(15, 8))
gs = gridspec.GridSpec(2, 3)
grid = itertools.product([0, 1], [0, 1, 2])

for clf, label, grd in zip(clf_list, rf_label, grid):
    scores = cross_val_score(clf, X, y, cv=3, scoring='accuracy')
    print("Accuracy: %.2f (+/- %.2f) [%s]" %(scores.mean(), scores.std(), label))
        
    clf.fit(X, y)
    ax = plt.subplot(gs[grd[0], grd[1]])
    fig = plot_decision_regions(X=X, y=y, clf=clf, legend=2)
    plt.title(label)

plt.show()

Accuracy: 0.63 (+/- 0.04) [DecisionTree]
Accuracy: 0.59 (+/- 0.08) [1 Trees]
Accuracy: 0.69 (+/- 0.02) [5 Trees]
Accuracy: 0.71 (+/- 0.06) [10 Trees]
Accuracy: 0.74 (+/- 0.04) [20 Trees]
Accuracy: 0.79 (+/- 0.02) [100 Trees]

可以看到隨着樹的增加，隨機森林的決策邊界越來越光滑，精度越來越高。而且當樹爲1個時隨機森林色精度是比決策樹要差的，這很容易理解，有可能重要的特徵被排除在外了，所以同樣是一棵樹，隨機森林分類效果更差。我們再通過誤差棒圖看一下決策樹數量對精度的影響

num_est = np.concatenate(([1],np.arange(5,105,5)))
rf_clf_cv_mean = []
rf_clf_cv_std = []

for n_est in num_est:    
    rf_clf = RandomForestClassifier(n_estimators=n_est, criterion='entropy', max_depth = 2)
    scores = cross_val_score(rf_clf, X, y, cv=3, scoring='accuracy')
    rf_clf_cv_mean.append(scores.mean())
    rf_clf_cv_std.append(scores.std())
    

plt.figure()
(_, caps, _) = plt.errorbar(num_est, rf_clf_cv_mean, yerr=rf_clf_cv_std, c='blue', fmt='-o', capsize=5)
for cap in caps:
    cap.set_markeredgewidth(1)                                                           
plt.ylabel('Accuracy')
plt.xlabel('Trees Number')
plt.show()

可以看到隨着樹的數量增加，隨機森林一開始有明顯的提升，但是當樹到達10個後就基本趨於穩定了

下面我們比較一下隨着學習器的增加，Tree Bagging和Random Forest的差異

num_est = np.concatenate(([1],np.arange(10,1005,10)))
bg_arr = []
rf_arr = []

for n_est in num_est:    
    rf_clf = RandomForestClassifier(n_estimators=n_est, criterion='entropy', max_depth = 2)
    rf_scores = cross_val_score(rf_clf, X, y, cv=3, scoring='accuracy')
    rf_arr.append(rf_scores.mean())
    
    bg_clf = BaggingClassifier(base_estimator=clf3, n_estimators=n_est, max_samples=0.8, max_features=0.8)
    bg_scores = cross_val_score(bg_clf, X, y, cv=3, scoring='accuracy')
    bg_arr.append(bg_scores.mean())
    

plt.figure()
plt.plot(num_est, bg_arr,  color='blue', label='Bagging')
plt.plot(num_est, rf_arr, color='red', label='Random Forest')
plt.legend()
plt.ylabel('Accuracy')
plt.xlabel('Learner Number')
plt.show()

可以看到Bagging和隨機森林有相似的收斂性，但是隨着個體學習器數量的增加，隨機森林的泛化性能要更好。

爲什麼說Bagging通過減小方差來提升精度

因爲Bagging是同質集成並且數據集也相同，因此各個學習器的偏差 $bias(h_t (x))$ 和方差 $Var(h_t (x))$ 近似相等，但是各個學習器並不是獨立的，於是
$bias(H(x)) = bias \left(\frac{1}{T}\sum_{t = 1}^{T} h_t (x)\right) = bias(h_t (x))$

$Var(H(x)) = Var \left(\frac{1}{T}\sum_{t = 1}^{T} h_t (x)\right) = \frac{1}{T}Var(h_t (x))$

所以Bagging不能降低偏差，而是會降低方差（因爲學習器之間的相關性，方差不會下降到 $\frac{1}{T}Var(h_t (x))$ ），如果各個學習器完全相同，即有 $\frac{1}{T}\sum_{t = 1}^{T} h_t (x) = h_t (x)$ ，那麼方差也不會降低，這也就是爲什麼要求學習器之間要有「多樣性」。隨機森林通過隨機選取特徵的方式降低了各決策樹的相關性，使得方差進一步降低。

Boosting

不同於bagging方法，boosting方法通過分步迭代（stage-wise）的方式來構建模型，在迭代的每一步構建的弱學習器都是爲了彌補已有模型的不足。Boosting算法先從初始訓練集訓練出一個基學習器，再根據基學習器的表現對訓練樣本進行調整，使得先前基學習器做錯的訓練樣本在後續受到更多關注，然後基於調整後的樣本分佈來訓練下一個基學習器；如此重複進行，直至基學習器數目達到事先指定的值T，最終將這T個基學習器進行加權。從偏差—方差分解的角度看，Boosting主要關注降低偏差，因此Boosting能基於泛化性能相當弱的學習器構建出很強的集成。下圖以Adaboost分類爲例簡要說明了Boosting的機制

Adaboost

Boosting算法最著名的代表是Adaboost（Adaptive Boosting）算法。Adaboost不改變所給的訓練數據，而不斷改變訓練數據的權值分佈，使得訓練數據在基本分類器的學習中起不同的作用，然後利用基本分類器的線性組合構建最終的分類器。Adaboost的二類分類算法的描述如下圖所示，其中 $y_i\in\{-1,+1\}$ ， $i = 1,2,...,m$ ， $f$ 是真實函數， $\mathcal{D}_t$ 是調整後用於進行第 $t$ 次訓練的樣本分佈， $h_t$ 是基於分佈 $\mathcal{D}_t$ 從數據集 $D$ 中訓練出的分類器， $\epsilon_t$ 是 $h_t$ 誤差的估計， $\alpha_t$ 是分類器 $h_t$ 的權重， $Z_t$ 是規範化因子，以確保 $\mathcal{D}_{t+1}$ 是一個分佈， $t=1,2,...,T$ 。

Adaboost是一種比較有特點的算法，可以總結如下：

每次迭代改變的是樣本的分佈，而不是重複採樣（re weight)。
樣本分佈的改變取決於樣本是否被正確分類：總是分類正確的樣本權值低，總是分類錯誤的樣本權值高（通常是邊界附近的樣本）。
最終的結果是弱分類器的加權組合，權值表示該弱分類器的性能。

我們通過一個例子來感受一下Adaboost的機制

實例

我們有如下數據集，其中有紅色圓形和綠色三角形共10個樣本，希望建立一個分類器，能判斷出任意位置樣本屬於哪一類，我們指定基學習器個數 $T$ 爲3

我們初始化樣本權值分佈 $\mathcal{D}_1 = {0.1, 0.1, \cdots, 0.1}$ ，並基於該分佈從訓練集中訓練出分類器 $h_1$ （理論上 $h_1$ 應使分類錯誤率達到最小，這裏僅用一條線代表分類結果，並未給出這條線是怎麼得到的，但是不代表這條線是隨意做出的），之後估計 $h_1$ 的誤差 $\epsilon_1$ ，確定 $h_1$ 的權重 $\alpha_1$ ，並更新樣本的權值分佈

重複上面操作依次訓練出基學習器 $h_2$ 、 $h_3$ ，此時基學習器個數達到我們預先指定的值。可以按照基學習器的權重給出輸出公式。並據此預測其他樣本的分類

每個區域是屬於哪個類別，由這個區域所在分類器的權值綜合決定。比如左下角的區域的樣本，屬於綠色分類區的權重爲 $h_1$ 中的0.42和 $h_2$ 中的0.65，其和爲1.07；屬於紅色分類區域的權重爲 $h_3$ 中的0.92；屬於紅色分類區的權重小於屬於綠色分類區的權值，因此左下角屬於綠

從結果中看，即使簡單的分類器，組合起來也能獲得很好的分類效果。

模型推導

Adaboost算法有多種推導方式，一種比較容易理解的是基於「加性模型」，即基學習器的線性組合

$H(x) = \sum_{t = 1}^{T} \alpha_t h_t(x)$

來最小化指數損失函數

$Loss(H | \mathcal{D}) = \mathbb{E}_{\mathcal{D}}[e^{-yH(x)}] = e^{-H(x)}P(y = 1 | x) + e^{H(x)}P(y = -1 | x)$

●什麼是前向分步算法（forward stagewise algorithm）？

前向分步算法是求解加法模型損失函數（ $\sum_{i = 1}^{N}Loss(y_i, \alpha h(x_i))$ ）最小化的一種算法，簡單來說就是在每一次學習一個基學習器，然後逐步逼近目標函數：

$\begin{array}{lcl} f_0(x) & = & 0 \\ f_1(x) & = & f_0(x) + \alpha_1 h_1(x) & = & \alpha_1 h_1(x) \\ f_2(x) & = & f_1(x) + \alpha_2 h_2(x) & = & \alpha_1 h_1(x) + \alpha_2 h_2(x) \\ \cdots \\ f_M(x) & = & f_{M-1}(x) + \alpha_M h_M(x) & = & \sum_{m = 1}^{M}\alpha_m h_m(x) \end{array}$

●爲什麼是指數損失函數？

指數損失函數對 $H(x)$ 求偏導得

$\frac{\partial{Loss(H | \mathcal{D})}}{\partial{H(x)}} = -e^{-H(x)}P(y = 1 | x) + e^{H(x)}P(y = -1 | x)$

令偏導數爲0有

$H(x) = \frac{1}{2} \ln \frac{P(y = 1 |x)}{P(y = -1 |x)}$

於是有

$\begin{aligned} sign\left(H(x)\right) & = sign\left(\frac{1}{2} \ln \frac{P(y = 1 |x)}{P(y = -1 |x)}\right) \\ & = \begin{cases} 0& \text{if } P(y = 1 |x) \gt P(y = -1 |x)\\ 1& \text{if } P(y = 1 |x) \lt P(y = -1 |x) \end{cases} \\ & = \mathop{\arg\max}_{y \in \{+1,-1\}} P(y | x) \end{aligned}$

過程中忽略了 $P(y = 1 |x) = P(y = -1 |x)$ 的情況，另外西瓜書中對於概率的寫法是 $P(f(x) = y|x)$ ，其中 $f$ 是真實函數。從結果可知， $sign\left(H(x)\right)$ 達到了貝葉斯最優錯誤率（ $if P(\omega_1|x) \gt P(\omega_1|x),\; then\ x\in \omega_1,\;else\ x\in \omega_2$ ），即最小化指數損失函數等價於分類錯誤率最小化。因此我們可以用它替代0/1損失函數作爲最優目標。有了以上結論，我們便可以進行算法的推導。

●推導過程

在算法的第 $t$ 次循環中，我們基於分佈 $\mathcal{D}_t$ 產生了基分類器 $h_t(x)$ 和分類器權重 $\alpha_t$ ，那麼它們應該使得 $H_t(x) = \sum_{k = 1}^{t}\alpha_k h_k(x)$ 最小化損失函數

$\begin{aligned} Loss(H_t(x) | \mathcal{D}_t) & = Loss(\sum_{k = 1}^{t}\alpha_k h_k(x) \ |\ \mathcal{D}_k) \\ & = Loss(H_{t-1}(x) \ +\ \alpha_t h_t(x)\ |\ \mathcal{D}_k) \\ & = \mathbb{E}_{\mathcal{D}_t} e^{-y (H_{t-1}(x) \ +\ \alpha_t h_t(x))} \\ & = \sum_{i = 1}^{N} e^{-y_i H_{t-1}(x_i)} e^{-y_i\alpha_t h_t(x_i)} \end{aligned}$

我們令 $w_{t,i} = e^{-y_i H_{t-1}(x_i)}$ ，顯然， $w_{t,i}$ 與 $\alpha_t$ 和 $h_t$ 無關，所以與最小化損失也無關，於是有

$Loss(H_t(x) | \mathcal{D}_t) = \sum_{i = 1}^{N} w_{t,i} e^{-y_i\alpha_t h_t(x_i)}$

設 $h_t^*$ 和 $\alpha_t^*$ 是滿足要求的分類器和分類器權重，即

$(h_t^*,\alpha_t^*) = \mathop{\arg \min}_{\alpha_t, h_t}Loss(H_t(x) | \mathcal{D}_t) = \mathop{\arg \min}_{\alpha_t, h_t}\sum_{i = 1}^{N} w_{t,i} e^{-y_i\alpha_t h_t(x_i)}$

對任意 $\alpha_t \gt 0$ ，考慮

$h_t^*(x) = \mathop{\arg \min}_{h_t} \sum_{i = 1}^{N}w_{t,i}I(h_t(x_i) \ne y_i)$

這意味着 $h_t^*(x)$ 是使第t輪加權訓練數據分類誤差最小的基本分類器，即爲所求。再求 $\alpha_t^*$

$Loss(H_t(x) | \mathcal{D}_t) = \sum_{i = 1}^{N} w_{t,i} e^{-y_i\alpha_t h_t^*(x_i)} = \sum_{h_t^*(x_i) \ne y_i}w_{t,i}e^{\alpha_t}+\sum_{h_t^*(x_i) = y_i}w_{t,i}e^{-\alpha_t}$

損失函數對 $\alpha_t$ 求偏導得

$\frac{\partial{Loss(H_t(x) | \mathcal{D}_t)}}{\partial{\alpha_t}} = -\sum_{h_t^*(x_i) \ne y_i}w_{t,i}e^{\alpha_t}+\sum_{h_t^*(x_i) = y_i}w_{t,i}e^{-\alpha_t}$

令偏導數爲0有

$\alpha_t^* = \frac{1}{2}\ln \frac{1-\sum_{h_t^*(x_i) \ne y_i}w_{t,i}}{\sum_{h_t^*(x_i) \ne y_i}w_{t,i}} = \frac{1}{2}\ln \frac{1-\epsilon_t}{\epsilon_t}$

再來看樣本權值的更新， $H_t(x_i) = H_{t-1}(x_i) + \alpha_t h_t(x_i)$ 兩邊同時乘 $-y_i$ 並進行指數運算得

$w_{t+1,i} = e^{-y_iH_t(x_i)} = e^{-y_iH_{t-1}(x_i)}e^{-y_i\alpha_t h_t(x_i)} = w_{t,i}e^{-y_i\alpha_t h_t(x_i)}$

之後爲確保 $\mathcal{D}_{t+1}$ 是一個分佈，再對權值進行規範化即可

爲什麼說Boosting通過減小偏差來提升精度

Boosting算法等價於用前向分步算法來最小化損失函數，以Adaboost爲例，算法用這樣的方式順序地最小化指數損失函數 $Loss(y,H_t(x))$ ，偏差自然降低。但也正是因爲這樣造成了子模型之間的強相關性，因此不能顯著降低方差。所以所Boosting主要通過降低偏差來提升精度

*多分類任務

參考：http://www.javashuo.com/article/p-hwlwdvps-nz.html

目前比較常用的多分類方法有三種

adaboost M1方法

該方法的主要思路是利用多分類的基學習器，最終選擇最有可能的分類結果

$H(x) = \mathop{\arg \max}_{y \in \mathcal{Y}} \sum_{t:h_t(x) = y}\alpha_t$

adaboost MH方法

應用「one-vs-all」的原理重構樣本空間，最終選擇最有可能的結果

對多分類輸出進行二進制編碼

類似adaboost MH方法，不過不需要重構樣本空間，而是把Label映射到 $K$ 位的二進制碼上（ $K$ 位分類數），然後訓練 $K$ 個二分類分類器，在解碼時生成 $K$ 位的二進制數。從而對應到一個label上。

編程（分類）

仍然使用Iris數據集，決策樹算法作爲個體學習器算法

import itertools
import numpy as np

import seaborn as sns
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec

from sklearn import datasets

from sklearn.tree import DecisionTreeClassifier

from sklearn.ensemble import BaggingClassifier
from sklearn.ensemble import RandomForestClassifier
from sklearn.ensemble import AdaBoostClassifier
from sklearn.model_selection import cross_val_score, train_test_split

from mlxtend.plotting import plot_learning_curves
from mlxtend.plotting import plot_decision_regions

iris = datasets.load_iris()
X, y = iris.data[:, 0:2], iris.target
    
clf = DecisionTreeClassifier(criterion='entropy', max_depth=2)

num_est = [1, 2, 3, 10]
label = ['AdaBoost (n_est=1)', 'AdaBoost (n_est=2)', 'AdaBoost (n_est=3)', 'AdaBoost (n_est=10)']

fig = plt.figure(figsize=(10, 8))
gs = gridspec.GridSpec(2, 2)
grid = itertools.product([0,1],repeat=2)

for n_est, label, grd in zip(num_est, label, grid):     
    boosting = AdaBoostClassifier(base_estimator=clf, n_estimators=n_est)  
    
    scores = cross_val_score(boosting, X, y, cv=3, scoring='accuracy')
    print("Accuracy: %.2f (+/- %.2f) [%s]" %(scores.mean(), scores.std(), label))
    
    boosting.fit(X, y)
    ax = plt.subplot(gs[grd[0], grd[1]])
    fig = plot_decision_regions(X=X, y=y, clf=boosting, legend=2)
    plt.title(label)

plt.show()

Accuracy: 0.63 (+/- 0.04) [AdaBoost (n_est=1)]
Accuracy: 0.56 (+/- 0.11) [AdaBoost (n_est=2)]
Accuracy: 0.73 (+/- 0.04) [AdaBoost (n_est=3)]
Accuracy: 0.71 (+/- 0.08) [AdaBoost (n_est=10)]

再來看看不同訓練計劃分下Adaboost的表現

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

boosting = AdaBoostClassifier(base_estimator=clf, n_estimators=10)
        
plt.figure()
plot_learning_curves(X_train, y_train, X_test, y_test, boosting, print_model=False, style='ggplot')
plt.show()

下面是不同基學習器個數下的Adaboost的精度

#Ensemble Size
num_est = np.concatenate(([1],np.arange(5,105,5)))
ada_clf_cv_mean = []
ada_clf_cv_std = []

for n_est in num_est:
    ada_clf = AdaBoostClassifier(base_estimator=clf, n_estimators=n_est)
    scores = cross_val_score(ada_clf, X, y, cv=3, scoring='accuracy')
    ada_clf_cv_mean.append(scores.mean())
    ada_clf_cv_std.append(scores.std())
    


plt.figure()
(_, caps, _) = plt.errorbar(num_est, ada_clf_cv_mean, yerr=ada_clf_cv_std, c='blue', fmt='-o', capsize=5)
for cap in caps:
    cap.set_markeredgewidth(1)                                                                        
plt.ylabel('Accuracy'); plt.xlabel('Ensemble Size'); plt.title('AdaBoost Ensemble');
plt.show()

運行過程中我們可以明顯的發現其和Bagging與RF的不同，就是一旦參數確定了結果就是一定的，沒有隨機性。這是因爲Bagging採樣的隨機性和RF特徵選擇的隨機性。我們再來比較一下這三個模型的精度隨着基學習數目增加精度的變化。

num_est = np.concatenate(([1],np.arange(10,100,10),np.arange(100,1001,50)))
bg_arr = []
rf_arr = []
ada_arr = []

for n_est in num_est:    
    rf_clf = RandomForestClassifier(n_estimators=n_est, criterion='entropy', max_depth = 2)
    rf_scores = cross_val_score(rf_clf, X, y, cv=3, scoring='accuracy')
    rf_arr.append(rf_scores.mean())
    
    bg_clf = BaggingClassifier(base_estimator=clf, n_estimators=n_est, max_samples=0.8, max_features=0.8)
    bg_scores = cross_val_score(bg_clf, X, y, cv=3, scoring='accuracy')
    bg_arr.append(bg_scores.mean())
    
    ada_clf = AdaBoostClassifier(base_estimator=clf, n_estimators=n_est)
    scores = cross_val_score(ada_clf, X, y, cv=3, scoring='accuracy')
    ada_arr.append(scores.mean())
    

plt.figure()
plt.plot(num_est, bg_arr,  color='blue', label='Bagging')
plt.plot(num_est, rf_arr, color='red', label='Random Forest')
plt.plot(num_est, ada_arr, color='green', label='Adaboost')
plt.legend()
plt.ylabel('Accuracy')
plt.xlabel('Learner Number')
plt.show()

可見以決策樹爲基學習器的Adaboost模型表現並不是很好

*訓練誤差分析

Adaboost最基本的性質是它能在學習過程中不斷減少訓練誤差，下面我們來證明這一性質。我們先令 $g(x) = \sum_{t = 1}{T} \alpha_t h_t(x)$ ，令 $w_{t,i}$ 表示第 $t$ 輪中樣本 $i$ 的權值，顯然有 $w_{1,i} = \frac{1}{N}$ ，則Adaboost最終分類器 $H(x)$ 的訓練誤差界爲：

$bias(H(x)) \ =\ \frac{1}{N} \sum_{i = 1}^{N} I(H(x_i) \ne y_i) \ \le\ \frac{1}{N} \sum_{i = 1}^{N} e^{-y_i g(x_i)} \ =\ \prod_{t = 1}^{T} Z_t$

證明：

當 $H(x_i) \ne y_i$ 時， $y_i g(x_i) \lt 0$ ，所以 $e^{-y_i g(x_i)} \gt 1 = I(H(x_i) \ne y_i)$ ，於是不等式的前半部分得證，下面只需證明 $\frac{1}{N} \sum_{i = 1}^{N} e^{-y_i g(x_i)} \ =\ \prod_{t = 1}^{T} Z_t$ ：

$\begin{aligned} \frac{1}{N} \sum_{i = 1}^{N} e^{-y_i g(x_i)} & = \frac{1}{N} \sum_{i = 1}^{N} e^{-y_i \sum_{t = 1}^{T} \alpha_t h_t(x_i)} \\ & = w_{1,i} \sum_{i = 1}^{N} e^{-\sum_{t = 1}^{T} \alpha_t h_t(x_i)y_i } \\ & = \sum_{i = 1}^{N} w_{1,i} \prod_{t = 1}^{T} e^{- \alpha_t h_t(x_i)y_i } \\ & = Z_1 \sum_{i = 1}^{N} \frac{w_{1,i} \prod_{t = 1}^{T} e^{- \alpha_t h_t(x_i)y_i }}{Z_1} \\ & = Z_1 \sum_{i = 1}^{N} w_{2,i} \prod_{t = 2}^{T} e^{- \alpha_t h_t(x_i)y_i } \\ & = Z_1 Z_2 \sum_{i = 1}^{N} w_{3,i} \prod_{t = 3}^{T} e^{- \alpha_t h_t(x_i)y_i } \\ & = \cdots \\ & = Z_1 Z_2 \cdots Z_{t-1} \sum_{i = 1}^{N} w_{T,i} e^{- \alpha_T h_T(x_i)y_i } \\ & = \prod_{t = 1}^{T} Z_t \end{aligned}$ N1i=1∑Ne−yig(xi)=N1i=1∑Ne−yi∑t=1Tαtht(xi)=w1,ii=1∑Ne−∑t=1Tαtht(xi)yi=i=1∑Nw1,it=1∏Te−αtht(xi)yi=Z1i=1∑NZ1w1,i∏t=1Te−αtht(xi)yi=Z1i=1∑Nw2,it=2∏Te−αtht(xi)yi=Z1Z2i=1∑Nw3,it=3∏Te−αtht(xi)yi=⋯=Z1Z2⋯Zt−1i=1∑NwT,ie−αThT(xi)yi=t=1∏TZt

■

通過這個定理，我們可以在每一輪都適當修改 $H_t$ 使得 $Z_t$ 最小，從而使訓練誤差下降最快。特別的，對於二分類問題，我們有以下結論：

$\prod_{t =1}^{T} Z_t \ =\ \prod_{t = 1}^{T}[2\sqrt{\epsilon_t(1-\epsilon_t)} \ = \ \prod_{t = 1}^{T} \sqrt{(1 - 4\gamma_t^2)} \ \le \ exp\left(-2 \sum_{t = 1}^{T} \gamma_t^2 \right)$

其中， $\gamma_t = \frac{1}{2} - \epsilon_t \gt 0$

證明：

先證明等式的部分：

$\begin{aligned} Z_t & = \sum_{i=1}^{N} w_{t,i} exp(-\alpha_t y_i h_t(x_i)) \\ & = \sum_{y_i = h_t(x_i)}w_{t, i} e^{-\alpha_t} \ +\ \sum_{y_i \ne h_t(x_i)} w_{t, i} e^{\alpha_t} \\ & = e^{-\alpha_t}\sum_{y_i = h_t(x_i)}w_{t, i} \ +\ e^{\alpha_t}\sum_{y_i \ne h_t(x_i)} w_{t, i} \\ & = e^{-\alpha_t}(1 - \epsilon_t) \ +\ e^{\alpha_t}\epsilon_t \\ & = {\left(\frac{1-\epsilon_t}{\epsilon_t}\right)}^{-\frac{1}{2}}(1 - \epsilon_t) \ +\ {\left(\frac{1-\epsilon_t}{\epsilon_t}\right)}^{\frac{1}{2}}\epsilon_t \\ & = 2\sqrt{\epsilon_t(1-\epsilon_t)} \\ & = \prod_{t = 1}^{T} \sqrt{(1 - 4\gamma_t^2)} \end{aligned}$ Zt=i=1∑Nwt,iexp(−αtyiht(xi))=yi=ht(xi)∑wt,ie−αt + yi̸=ht(xi)∑wt,ieαt=e−αtyi=ht(xi)∑wt,i + eαtyi̸=ht(xi)∑wt,i=e−αt(1−ϵt) + eαtϵt=(ϵt1−ϵt)−21(1−ϵt) + (ϵt1−ϵt)21ϵt=2ϵt(1−ϵt) =t=1∏T(1−4γt2)

不等式部分可由 $\sqrt{1-x}$ 和 $e^x$ 在 $x = 0$ 處的泰勒展開式推導得出

■

進一步，若存在 $\gamma \gt 0$ 使得對所有 $t$ 有 $\gamma_t \ge \gamma$ ，那麼有：

$bias(H(x)) \le exp\left( -2T\gamma^2 \right)$

這表明Adaboost在這種條件下的訓練誤差是以指數速率下降的

*過擬合分析

這裏僅給出結論——Adaboost不會發生過擬合，這是Adaboost的另一個性質。「margin theory」可以比較直觀地解釋這一性質。

總結

●優點

adaboost是一種有很高精度的分類器。
可以使用各種方法構建子分類器，adaboost算法提供的是框架。
當使用簡單分類器時，計算出的結果是可以理解的。而且弱分類器構造極其簡單。
簡單，不用做特徵篩選。
不用擔心overfitting！
隨着迭代次數的增加，實際上錯誤率上界在下降。

●缺點

對異常值非常敏感
對損失函數形式有要求

●應用

用於二分類或多分類的應用場景
用於做分類任務的baseline
用於特徵選擇（feature selection)
用於修正badcase，

提升樹

提升樹（Boosting Tree）是以二叉分類樹或二叉迴歸樹爲基本分類器的提升方法，採用加法模型（即基函數的線性組合）與前向分步算法。提升樹有很多不同的算法，其主要區別在於使用的損失函數不同，包括用平方誤差損失函數的迴歸問題，用指數函數損失的分類問題，以及用一般損失函數的一般決策問題。

提升樹模型爲：

$f_M(x) = \sum_{m = 1}^{M}T(x; \Theta_m)$

其中 $T(x; \Theta_m)$ 爲決策樹； $\Theta_m$ 爲決策樹的參數； $M$ 爲樹的個數。首先確定初始模型：

$f_0(x) = 0$

第 $m$ 步的模型是：

$f_m(x) = f_{m-1}(x) + T(x; \Theta_m),\quad m = 1, 2, \cdots, M$

通過最小化損失函數確定決策樹 $T(x; \Theta_m)$ 的參數：

$\hat{\Theta}_m = \mathop{\arg \min}_{\Theta_m}Loss(y, f_m(x)) = \mathop{\arg \min}_{\Theta_m}Loss(y, f_{m-1}(x) + T(x; \Theta_m))$

以迴歸問題的提升樹爲例，其損失函數是平方誤差損失

$\begin{aligned} Loss(y, f_m(x)) & = Loss(y, f_{m-1}(x) + T(x; \Theta_m)) \\ & = {(y - f_{m-1}(x) - T(x; \Theta_m))}^2 \\ & = {(r_m - T(x; \Theta_m))}^2 \end{aligned}$

其中， $r_m$ 是當前模型擬合數據的殘差，因此我們只需根據 $(x, r_m)$ 來學習一個迴歸樹 $T(x; \Theta_m)$ ，進而確定 $f_m(x) = f_{m-1}(x) + T(x; \Theta_m)$ ，最終學得 $f_M(x) = \sum_{m = 1}^{M}T(x; \Theta_m)$

實例

考慮用如下數據集學習一個提升樹模型（深度爲1）

$x_i$	1	2	3	4	5	6	7
$y_i$	4.5	4.7	4.9	5.3	5.8	7.0	7.9

解第1步求 $f_1(x)$ 即迴歸樹 $T(x;\Theta_1)$

根據已知的數據，我們可以考慮如下切分點：

$1.5, \;2.5, \;3.5, \;4.5, \;5.5, \;6.5$

回憶二叉迴歸樹的解法，容易求出各切分點 $s$ 下

$g(s) = \mathop{\min}_{c_1}\sum_{x\in R_1}{(y_i - c_1)}^2 + \mathop{\min}_{c_2}\sum_{x\in R_2}{(y_i - c_2)}^2$

其中 $R_1 = \{x \le s\}$ ， $R_2 = \{x \gt s\}$ ， $c_1$ 、 $c_2$ 易知爲 $c_1 = \bar{y}I(x \in R_1)$ ， $c_2 = \bar{y}I(x \in R_2)$

當 $s = 1.5$ 時， $c_1 = 4.5$ ， $c_2 = (4.7+4.9+5.3+5.8+7+7.9)/6 = 5.93$ ， $g(s) = 0 + (4.7-5.93)^2+ (4.9-5.93)^2+ (5.3-5.93)^2+ (5.8-5.93)^2+ (7-5.93)^2+ (7.9-5.93)^2 = 8.0134$

當 $s = 2.5$ 時， $c_1 = (4.5+4.7)/2 = 4.6$ ， $c_2 = (4.9+5.3+5.8+7+7.9)/5 = 6.18$ ， $g(s) = (4.5-4.6)^2 + (4.7-4.6)^2+ (4.9-6.18)^2+ (5.3-6.18)^2+ (5.8-6.18)^2+ (7-6.18)^2+ (7.9-6.18)^2 = 6.208$

當 $s = 3.5$ 時， $c_1 = (4.5+4.7+4.9)/3 = 4.7$ ， $c_2 = (5.3+5.8+7+7.9)/4 = 6.5$ ， $g(s) = (4.5-4.7)^2 + (4.7-4.7)^2+ (4.9-4.7)^2+ (5.3-6.5)^2+ (5.8-6.5)^2+ (7-6.5)^2+ (7.9-6.5)^2 = 4.22$

當 $s = 4.5$ 時， $c_1 = (4.5+4.7+4.9+5.3)/4 = 4.85$ ， $c_2 = (5.8+7+7.9)/3 = 6.90$ ， $g(s) = (4.5-4.85)^2 + (4.7-4.85)^2+ (4.9-4.85)^2+ (5.3-4.85)^2+ (5.8-6.9)^2+ (7-6.9)^2+ (7.9-6.9)^2 = 2.57$

當 $s = 5.5$ 時， $c_1 = (4.5+4.7+4.9+5.3+5.8)/5 = 5.04$ ， $c_2 = (7+7.9)/2 = 7.45$ ， $g(s) = (4.5-5.04)^2 + (4.7-5.04)^2+ (4.9-5.04)^2+ (5.3-5.04)^2+ (5.8-5.04)^2+ (7-7.45)^2+ (7.9-7.45)^2 = 1.477$

當 $s = 6.5$ 時， $c_1 = (4.5+4.7+4.9+5.3+5.8+7)/6 = 5.37$ ， $c_2 = 7.9$ ， $g(s) = (4.5-5.37)^2 + (4.7-5.37)^2+ (4.9-5.37)^2+ (5.3-5.37)^2+ (5.8-5.37)^2+ (7-5.37)^2+ 0 = 4.2734$

$s$	1.5	2.5	3.5	4.5	5.5	6.5
$g(s)$	8.01	6.21	4.22	2.57	1.48	4.27

所以切分點 $s$ 應爲5.5，即

$T_1(x) = \begin{cases} 5.04, & x \le 5.5 \\ 7.45, & x \gt 5.5 \end{cases}$

於是 $f_1(x) = T_1(x)$

用 $f_1(x)$ 擬合訓練數據的平方損失爲

$Loss(y, f_1(x)) = \sum_{i = 1}^{7}{(y_i - f_1(x_i))}^2 = 1.477$

下面用 $f_1(x)$ 擬合殘差 $r_{2,i} = y_i - f_1(x_i)$ ：

$x_i$	1	2	3	4	5	6	7
$r_{2,i}$	-0.54	-0.34	-0.14	0.26	0.76	-0.45	0.45

接着求 $T_2(x)$ ，方法同求 $T_1(x)$ ，得

$s$	1.5	2.5	3.5	4.5	5.5	6.5
$g(s)$	1.14	0.92	0.87	1.14	1.48	1.24

所以切分點爲3.5，即

$T_2(x) = \begin{cases} -0.34, & x \le 3.5 \\ 0.255, & x \gt 3.5 \end{cases}$

於是 $f_2(x) = f_1(x) + T_2(x)$ ，有

$f_2(x) = \begin{cases} 4.7, & x \le 3.5 \\ 5.295, & 3.5 \lt x \le 5.5 \\ 7.705, & x \gt 5.5 \end{cases}$

用 $f_2(x)$ 擬合訓練數據的平方損失爲

$Loss(y, f_2(x)) = \sum_{i = 1}^{7}{(y_i - f_2(x_i))}^2 = 0.870$

假設我們提前設定決策樹的數量就是2，那麼 $f_2(x)$ 即爲所求，否則繼續重複前面的運算直到決策樹數量達到預先設定值 $M$ 或者 $f_m(x)$ 擬合訓練數據的平方誤差損失滿足要求

梯度提升

梯度提升（Gradient boosting）是一種用於迴歸、分類和排序任務的機器學習技術，屬於Boosting算法族的一部分。梯度提升通過集成多個弱學習器，通常是決策樹，來構建最終的預測模型，在提升樹的基礎上，當損失函數不能使用便於計算的平方損失或者指數損失時，利用損失函數的負梯度

$-{\left[\frac{\partial{Loss(y, f(x_i))}}{\partial{f(x_i}}\right]}_{f(x) = f_{m-1}(x)}$

作爲迴歸問題提升樹算法中的殘差的近似值，即

$r_{m,i} = -{\left[\frac{\partial{Loss(y, f(x_i))}}{\partial{f(x_i}}\right]}_{f(x) = f_{m-1}(x)}$

GBDT

參考：https://www.zybuluo.com/yxd/note/611571#gbdt算法

基於梯度提升算法的學習器叫做GBM(Gradient Boosting Machine)。理論上，GBM可以選擇各種不同的學習算法作爲基學習器。現實中，用得最多的基學習器是決策樹。爲什麼梯度提升方法傾向於選擇決策樹（通常是CART樹）作爲基學習器呢？這與決策樹算法自身的優點有很大的關係。決策樹可以認爲是if-then規則的集合，易於理解，可解釋性強，預測速度快。同時，決策樹算法相比於其他的算法需要更少的特徵工程，比如可以不用做特徵標準化，可以很好的處理字段缺失的數據，也可以不用關心特徵間是否相互依賴等。決策樹能夠自動組合多個特徵，學習特徵之間更高級別的相互關係，它可以毫無壓力地處理特徵間的交互關係並且是非參數化的，因此你不必擔心異常值或者數據是否線性可分（舉個例子，決策樹能輕鬆處理好類別A在某個特徵維度x的末端，類別B在中間，然後類別A又出現在特徵維度x前端的情況）。此外，模型很容易實現擴展。不過，單獨使用決策樹算法時，有容易過擬合缺點。所幸的是，通過各種方法，抑制決策樹的複雜性，降低單顆決策樹的擬合能力，再通過梯度提升的方法集成多個決策樹，最終能夠很好的解決過擬合的問題。由此可見，梯度提升方法和決策樹學習算法可以互相取長補短，是一對完美的搭檔。至於抑制單顆決策樹的複雜度的方法有很多，比如限制樹的最大深度、限制葉子節點的最少樣本數量、限制節點分裂時的最少樣本數量、吸收bagging的思想對訓練樣本採樣（subsample），在學習單顆決策樹時只使用一部分訓練樣本、借鑑隨機森林的思路在學習單顆決策樹時只採樣一部分特徵、在目標函數中添加正則項懲罰複雜的樹結構（XGboost）等。現在主流的GBDT算法實現中這些方法基本上都有實現，因此GBDT算法的超參數還是比較多的，應用過程中需要精心調參，並用交叉驗證的方法選擇最佳參數。

具體的算法如下：

輸入：訓練數據集 $T = \{(x_1,y_1),(x_2,y_2),\cdots,(x_N,y_N)\}$ ， $x_i \in \mathcal{X}\subseteq \mathbb{R}^n$ ， $y_i \in \mathcal{Y}\subseteq \mathbb{R}$ ；損失函數 $Loss(y,f(x))$ ；

輸出：迴歸樹 $\hat{f}(x)$

(1) 初始化

$f_0(x) = \mathop{\arg \min}_{c} \sum_{i = 1}^{N} Loss(y_i, c)$

(2) For $m = 1, 2, \cdots, M$ do:

(a) 計算

$r_{m,i} = -{\left[\frac{\partial{Loss(y, f(x_i))}}{\partial{f(x_i}}\right]}_{f(x) = f_{m-1}(x)}$

(b) 對(x_i, r_{m,i})擬合一個迴歸樹，得到第m棵樹的葉節點區域 $R_{m,j}$ ， $j = 1, 2, \cdots, J$

$c_{m,j} = \mathop{\arg\min}_{c} \sum_{x_i \in R_{m,j}}Loss(y_i, f_{m-1}+c)$

(d) 更新 $f_m(x) = f_{m-1}(x)+\sum_{j = 1}^{J}c_{m,j}I(x \in R_{m,j})$

(3) 得到迴歸樹

$\hat{f}(x) = f_M(x) = \sum_{m = 1}^{M}\sum_{j = 1}^{J}c_{m,j}I(x \in R_{m,j})$

GBDT實例

我們仍以前面的例子爲例，損失函數也仍然選擇平方誤差損失，我們看看結果是否會有不同

解第一步先初始化 $f_0(x)$ ，易知使平方誤差損失達到最小的 $c$ 值爲 $c_0 = \frac{1}{7}\sum_{i = 1}^{7}y_i = 5.73$ ，於是 $f_0(x) = c_0 = 5.73$

接着擬合第一輪的殘差

$r_{1,i} = -{\left[\frac{\partial{Loss(y_i, f(x_i))}}{\partial{f(x_i)}}\right]}_{f_0(x)} = 2{[y_i-f(x_i)]}_{f_0(x)}$

於是有

$x_i$	1	2	3	4	5	6	7
$r_{1,i}$	-2.45	-2.06	-1.66	-0.86	0.14	2.54	4.34

對 $r_{1,i}$ 擬合迴歸樹，易求得各切分點下的 $g(s)$ 爲

$s$	1.5	2.5	3.5	4.5	5.5	6.5
$g(s)$	32.05	24.80	16.87	10.26	5.89	17.06

於是確定切分點 $s = 5.5$ ，葉節點區域 $R_{1,1} = \{x | x \le 5.5\}$ ， $R_{1,2} = \{x | x \gt 5.5\}$ ，有迴歸樹

$T_1(x) = \begin{cases} -1.38, & x \le 5.5 \\ 3.44, & x \gt 5.5 \end{cases}$

進而有

$c_{1,1} = \frac{1}{5}\sum_{i = 1}^{5}[y_i - f_0(x_i)] = \frac{1}{5}\sum_{i = 1}^{5}y_i - f_0(x) = -0.69$

$c_{1,2} = \frac{1}{2}\sum_{i = 6,7}[y_i - f_0(x_i)] = \frac{1}{5}\sum_{i = 6,7}y_i - f_0(x) = 1.72$

最後學得

$f_1(x) = f_0(x) + \sum_{j = 1}^{2}c_{1,j}I(x \in R_{1,j}) = \begin{cases} 5.04, & x \le 5.5 \\ 7.45, & x \gt 5.5 \end{cases}$

用 $f_1(x)$ 擬合訓練數據的平方損失誤差：

$Loss(y,f_1(x)) = \sum_{i = 1}^{7}{[y_i-f_1(x_i)]}^2 = 1.477$

至此我們可以看到雖然學習過程不同，但是學得的 $f_1(x)$ 和提升樹是相同的，那麼繼續下一輪的學習，計算

$r_{2,i} = -{\left[\frac{\partial{Loss(y_i, f(x_i))}}{\partial{f(x_i)}}\right]}_{f_1(x)} = 2{[y_i-f(x_i)]}_{f_1(x)}$

於是有

$x_i$	1	2	3	4	5	6	7
$r_{2,i}$	-1.08	-0.68	-0.28	0.52	1.52	-0.9	0.9

對 $r_{2,i}$ 擬合迴歸樹，易求得各切分點下的 $g(s)$ 爲

$s$	1.5	2.5	3.5	4.5	5.5	6.5
$g(s)$	4.55	3.74	3.48	4.56	5.89	4.96

於是確定切分點 $s = 3.5$ ，葉節點區域 $R_{2,1} = \{x | x \le 3.5\}$ ， $R_{2,2} = \{x | x \gt 3.5\}$ ，有迴歸樹

$T_1(x) = \begin{cases} -0.68, & x \le 3.5 \\ 0.51, & x \gt 3.5 \end{cases}$

進而有

$c_{2,1} = \frac{1}{3}\sum_{i = 1}^{3}[y_i - f_1(x_i)] = \frac{1}{3}\sum_{i = 1}^{3}y_i - f_1(x) = -0.34$

$c_{2,2} = \frac{1}{4}\sum_{i = 4}^{7}[y_i - f_1(x_i)] = \frac{1}{4}\left[ \sum_{i = 4}^{7}y_i - \sum_{i = 4}^{7}f_1(x_i) \right] = 0.255$

最後學得

$f_2(x) = f_1(x) + \sum_{j = 1}^{2}c_{2,j}I(x \in R_{1,j}) = \begin{cases} 4.70, & x \le 3.5 \\ 5.30, & 3.5 \lt x \le 5.5 \\ 7.70, & x \gt 5.5 \end{cases}$

用 $f_2(x)$ 擬合訓練數據的平方損失誤差：

$Loss(y,f_2(x)) = \sum_{i = 1}^{7}{[y_i-f_2(x_i)]}^2 = 0.870$

最終結果仍和前面是一樣的，這說明提升樹用平方誤差損失實質上也是利用了損失函數的負梯度，提升樹用殘差進行學習本質上仍然是計算負梯度（因爲平方誤差損失的負梯度恰好是殘差），如果使用其他損失函數就不對了。也就是說殘差學習是負梯度學習的一個特例，或者說提升樹是梯度提升樹的一個特例，但是顯然梯度提升適用範圍更廣，在實際中梯度提升可以使用任意可微的損失函數

GBDT爲什麼要用負梯度來近似殘差？

首先我覺得這個問題本身就是錯誤的，我們並沒有這樣做過。GBDT的本來就是在每一輪迭代中對損失函數的負梯度構造迴歸樹，可以證明，這樣能夠使學習器收斂爲一個強學習器（損失函數最小）。而當損失函數爲平方誤差損失時，恰好損失函數的負梯度就是殘差。同時對於損失函數爲平方誤差的例子，直接用殘差來進行解釋反而更容易理解，並且由於提出的更早所以變成了負梯度近似殘差這一說法。就好像我們做一道數學選擇題，我們可能用特殊值法很容易得到了答案，但是並沒有接觸這道題本質的解法，如果題幹變化了（損失函數不再是平方誤差損失）或者題型改爲解答題，那麼特殊值法也不再適用。

GBDT編程（迴歸）

參考：https://scikit-learn.org/stable/auto_examples/ensemble/plot_gradient_boosting_regression.html#sphx-glr-auto-examples-ensemble-plot-gradient-boosting-regression-py

用波士頓房價數據集做GBDT迴歸

import numpy as np
import matplotlib.pyplot as plt

from sklearn.ensemble import GradientBoostingRegressor
from sklearn import datasets
from sklearn.utils import shuffle
from sklearn.metrics import mean_squared_error

加載數據

boston = datasets.load_boston()
X, y = shuffle(boston.data, boston.target, random_state=13)
X = X.astype(np.float32)
offset = int(X.shape[0] * 0.9)
X_train, y_train = X[:offset], y[:offset]
X_test, y_test = X[offset:], y[offset:]

訓練模型

params = {'n_estimators': 500, 'max_depth': 4, 'min_samples_split': 2,
          'learning_rate': 0.01, 'loss': 'ls'}
clf = GradientBoostingRegressor(**params)

clf.fit(X_train, y_train)
mse = mean_squared_error(y_test, clf.predict(X_test))
print("MSE: %.4f" % mse)

MSE: 6.5401

繪製偏差圖

# 計算測試集偏差
test_score = np.zeros((params['n_estimators'],), dtype=np.float64)

for i, y_pred in enumerate(clf.staged_predict(X_test)):
    test_score[i] = clf.loss_(y_test, y_pred)

plt.figure(figsize=(12, 6))
plt.subplot(1, 2, 1)
plt.title('Deviance')
plt.plot(np.arange(params['n_estimators']) + 1, clf.train_score_, 'b-',
         label='Training Set Deviance')
plt.plot(np.arange(params['n_estimators']) + 1, test_score, 'r-',
         label='Test Set Deviance')
plt.legend(loc='upper right')
plt.xlabel('Boosting Iterations')
plt.ylabel('Deviance')

Text(0, 0.5, ‘Deviance’)

繪製特徵重要性圖

feature_importance = clf.feature_importances_
# make importances relative to max importance
feature_importance = 100.0 * (feature_importance / feature_importance.max())
sorted_idx = np.argsort(feature_importance)
pos = np.arange(sorted_idx.shape[0]) + .5

plt.barh(pos, feature_importance[sorted_idx], align='center')
plt.yticks(pos, boston.feature_names[sorted_idx])
plt.xlabel('Relative Importance')
plt.title('Variable Importance')
plt.show()

梯度提升 Vs Adaboost

參考：https://www.zybuluo.com/yxd/note/611571#gbdt算法

與AdaBoost算法不同，梯度提升方法在迭代的每一步構建一個能夠沿着梯度最陡的方向降低損失（steepest-descent）的學習器來彌補已有模型的不足，而傳統的boosting算法只能是儘量向梯度方向減小（在Adaboost實例中能夠看到隨着學習器數量增加精度有起伏的情況），這是GBDT與傳統boosting算法最大的區別，這也是爲什麼GBDT相比傳統boosting算法可以用更少的樹個數與深度達到更好的效果。經典的AdaBoost算法只能處理採用指數損失函數的二分類學習任務等，而梯度提升方法可以使用更多種類的可微損失函數來處理各類學習任務（多分類、迴歸、排序等），應用範圍大大擴展。另一方面，AdaBoost算法對異常點（outlier）比較敏感（作爲錯分值而受到較大關注），而梯度提升算法通過引入bagging思想、加入正則項等方法能夠有效地抵禦訓練數據中的噪音，具有更好的健壯性。這也是爲什麼梯度提升算法（尤其是採用決策樹作爲弱學習器的GBDT算法）如此流行的原因。

XGboost

參考：https://arxiv.org/pdf/1603.02754v1.pdf

XGboost 的全稱是Extreme Gradient Boosting，由華盛頓大學的陳天奇博士提出

●準備知識

在推導算法之前，先回憶一下目標函數（objective function）的概念，目標函數通常定義爲如下形式：

$Obj(\Theta) = L(\Theta) + \Omega(\Theta)$

其中， $L(\Theta)$ 是損失函數，用來衡量模型擬合訓練數據的好壞程度； $\Omega(\Theta)$ 稱之爲正則項，用來衡量學習到的模型的複雜度，按照奧卡姆剃刀的原則，因此模型越簡單越好。常用的損失函數有平方損失（square loss）： $L(y,\hat{y}) = \sum_{i = 1}^{N}{(y_i - \hat{y}_i)}^2$ ；Logistic損失： $L(y,\hat{y}) = - \sum_{i = 1}^{N} [y_i \log sigmoid(\hat{y}_i) + (1 - y_i) \log sigmoid(1 - \hat{y}_i)]$ 。常用的正則項有L1範數： $\Omega(w) = \lambda ||w||_1$ 和L2範數： $\Omega(w) = \lambda ||w||^2$ 。Ridge regression就是指使用平方損失和L2範數正則項的線性迴歸模型；Lasso regression就是指使用平方損失和L1範數正則項的線性迴歸模型；邏輯迴歸（Logistic regression）指使用logistic損失和L2範數或L1範數正則項的線性模型。

再回顧一下泰勒公式的一些知識，根據泰勒公式把函數 $f(x + \triangle x)$ 在 $x$ 點處二階展開，可得到如下等式：

$f(x + \triangle x) \approx f(x) + f'(x)\triangle x + \frac{1}{2}f''(x){\triangle x}^2$

那麼對於一棵決策樹，我們如何描述它的複雜度呢，假設 $f_m(x)$ 是一棵已經建好的決策樹，它有 $T$ 個葉子結點， $w_j$ 是葉子結點的權重， $j = 1, 2, \cdots, T$ ，用 $q$ 表示樹的結構，它把特徵維度爲 $d$ 的數據映射到葉子結點的索引： $q: \mathbb{R}^d \to \{1, 2, \cdots, T\}$ ，即 $q(x)$ 表示數據 $x$ 在樹中的位置

於是 $f(x)$ 可以表示爲：

$f_m(x) = w_{q(x)}, \quad, w \in \mathbb{R}^T, \quad, q: \mathbb{R}^d \to \{1, 2, \cdots, T\}$

有了公式化的決策樹，我們便可以用樹的葉子數量和葉子權重的平滑程度（用權重的L2範數表示）來定義樹的複雜度

$\Omega(f_m(x)) = \gamma T + \frac{1}{2}\lambda \sum_{j = 1}^{T}w_j^2$

●模型

有了前面的準備，我們可以進入XGboost的模型推導，在前向分步算法中，第m步的目標函數爲

$\begin{aligned} Obj(\Theta)^{(m)} & = L(\Theta) + \Omega(\Theta) \\ & = \sum_{i = 1}^{N}l(y_i,{\hat{y}_i}^{(m)}) + \sum_{j = 1}^{m} \Omega (f_j) \\ & = \sum_{i = 1}^{N}l(y_i,{\hat{y}_i}^{(m-1)} + f_m(x_i)) + \Omega (f_m) +Constant1 \end{aligned}$

其中 $Constant1 = \sum_{j = 1}^{m-1} \Omega (f_j)$ ，因其不參與第 $m$ 輪的計算，所以可視爲常數。

由泰勒公式

$\begin{aligned} Obj(\Theta)^{(t)} & \approx \sum_{i = 1}^{N}[l(y_i,{\hat{y}_i}^{(m-1)}) + l'(y_i,{\hat{y}_i}^{(m-1)})f_m(x_i) + \frac{1}{2}l''(y_i,{\hat{y}_i}^{(m-1)})f_m^2(x_i)] + \Omega (f_m) +Constant1 \\ & \triangleq \sum_{i = 1}^{N}[l(y_i,{\hat{y}_i}^{(m-1)}) + g_i f_m(x_i) + \frac{1}{2}h_i f_m^2(x_i)] + \Omega (f_m) +Constant1 \\ & = \sum_{i = 1}^{N}[g_i f_m(x_i) + \frac{1}{2}h_i f_m^2(x_i)] + \Omega (f_m) + Constant1 + Constant2 \end{aligned}$

其中 $Constant2 = l(y_i,{\hat{y}_i}^{(m-1)})$ 同樣對於目標函數求最優值無影響，所以目標函數可以進一步簡化爲

$\begin{aligned} Obj(\Theta)^{(m)} & \approx \sum_{i = 1}^{N}[g_i f_m(x_i) + \frac{1}{2}h_i f_m^2(x_i)] + \Omega (f_m) \\ & = \sum_{i = 1}^{N}[g_i w_{q(x_i)} + \frac{1}{2}h_i w_{q(x_i)}^2] + \gamma T + \frac{1}{2}\lambda \sum_{j = 1}^{T}w_j^2 \end{aligned}$

定義 $I_j = \{i | q(x_i) = j\}$ 表示第 $j$ 個葉子結點裏的樣本序號的集合，於是目標函數又可表示爲

$\begin{aligned} Obj(\Theta)^{(m)} & = \sum_{j = 1}^{T}\left[\sum_{i \in I_j}\left(g_i w_j + \frac{1}{2}h_i w_j^2\right)\right] + \gamma T + \frac{1}{2}\lambda \sum_{j = 1}^{T}w_j^2 \\ & = \sum_{j = 1}^{T} \left[\left(\sum_{i \in I_j}g_i \right)w_j + \frac{1}{2}\left(\sum_{i \in I_j}h_i + \lambda\right)w_j^2\right] + \gamma T \\ & \triangleq \sum_{j = 1}^{T} [G_j w_j + \frac{1}{2}(H_j + \lambda)w_j^2] + \gamma T \end{aligned}$

假設我們的樹的結構 $q(x)$ 已經確定，爲求目標函數最小值，令目標函數對 $w_j$ 求偏導，令偏導數爲0得

$w_j^* = - \frac{G_j}{H_j + \lambda}$

於是

$\min_{w}Obj(\Theta)^{(m)} = - \frac{1}{2}\sum_{j = 1}^{T} \frac{G_j^2}{H_j + \lambda} + \gamma T$

那麼現在現在只要知道樹的結構，就能得到一個該結構下的最好分數（ $Obj$ ）。可是樹的結構應該怎麼確定呢？沒法用枚舉，畢竟可能的狀態基本屬於無窮種。這種情況，貪婪算法是個好方法。從樹的深度爲0開始，每一節點都遍歷所有的特徵。對於某個特徵，先按照該特徵裏的值進行排序（複雜度爲 $O(n\log n)$ ），然後線性掃描（複雜度爲 $O(n)$ ）該特徵來決定最好的分割點，最後在所有特徵（複雜度爲 $O(d(n \log n +n)))$ ）裏選擇分割後，分數增加（Gain）最高的那個特徵。

$Obj_{split} = − \frac{1}{2}[\frac{G^2_L}{H_L+\lambda}+\frac{G^2_R}{H_R+\lambda}]+\gamma T_{split}$

$Obj_{no-split} = − \frac{1}{2}\frac{{G_L + G_R}^2}{H_L+H_R +\lambda}+\gamma T_{no-split}$

$Gain=Obj_{no-split}−Obj_{split}=\frac{1}{2}[\frac{G^2_L}{H_L+\lambda}+\frac{G^2_R}{H_R+\lambda} - \frac{{G_L + G_R}^2}{H_L+H_R +\lambda}]−\gamma(T_{split}−T_{no-split})$

這時，就有兩種後續。一種是當最好的分割的情況下， $Gain$ 爲負時就停止樹的生長（類似前剪枝），這樣的話效率會比較高也簡單，但是這樣就放棄了未來可能會有更好的情況（因爲引入了正則項，所以增加分支並不能保證 $Gain$ 爲正）。另外一種就是一直分割到最大深度，然後進行修剪，遞歸得把劃分葉子得到的 $Gain$ 爲負的收回（類似後剪枝）。一般來說，後一種要好一些，於是我們採用後一種，完整的（無修剪）算法如下（若樹的深度爲 $k$ ，那麼總的複雜度爲 $O(kdn \log n)$ ）

此外還有近似算法，該算法更適合對數據量較大、難以計算的任務，算法描述如下，這裏不做過多介紹

●注意

對某個節點的分割時，是需要按某特徵的值排序，那麼對於無序的類別變量，就需要進行one-hot化。不然的話，假設某特徵有1，2，3三種變量。我們進行比較時，就會只比較左子樹爲1,2或者右子樹爲2,3，或者不分割，哪個更好。這樣的話就沒有考慮，左子樹爲1,3的分割。因爲 $Gain$ 於特徵的值範圍是無關的，它採用的是已經生成的樹的結構與權重來計算的。所以不需要對特徵進行歸一化處理。

●總述

我們回顧一下求解模型的整個過程來看看在具體應用XGboost時需要哪些操作

輸出：迴歸樹 $\hat{f}(x)$

(1) 初始化 $f_0(x) = 0$

(2) For $m = 1, 2, \cdots, M$ do:

a) 計算

$g_i = \frac{\partial{l(y_i,{\hat{y}_i}^{(m-1)})}}{\partial{{\hat{y}_i}^{(m-1)}}}$

$h_i = \frac{\partial^2{l(y_i,{\hat{y}_i}^{(m-1)})}}{{\left(\partial{\hat{y}_i}^{(m-1)}\right)}^2}$

b) 計算

$Obj = - \frac{1}{2}\sum_{j = 1}^{T} \frac{G_j^2}{H_j + \lambda} + \gamma T$

c) 採用貪婪算法生成樹 $f_m(x)$

d) 更新 $\hat{y}^{m} = \hat{y}^{m-1} + \epsilon_m f_m(x)$ ，其中 $\epsilon_m$ 是一個收縮係數，或者叫做步進。加上它的好處就是每一步我們都不是做一個完全的最優化，留下餘地給未來的循環，這樣能防止過擬合。

(3) 得到迴歸樹

$\hat{f}(x) = \sum_{m = 1}^{M}\epsilon_m f_m(x)$

XGboost編程

官網：https://xgboost.readthedocs.io/en/latest/python/python_intro.html

數據接口

xgboost可以把不同類型的數據轉化爲模塊內置的DMatrix數據類型：xgboost.DMatrix(data, label=None, missing=None, weight=None, silent=False, feature_names=None, feature_types=None, nthread=None)，各參數的說明如下：

參數名	參數說明	備註
data	數據來源：可以是文本文件/csv文件/二進制文件/np.array/scipy.sparse/pd.DataFrame/dt.Frame	1. 注意對於分類特徵，即非數值型特徵要轉化成one-hot編碼。2. 對於有標題的csv文件，不能直接導入，要先轉化爲pdDataFrame格式
label	標籤列
missing	缺失值	如果不設置，默認爲np.nan
weight	樣本權值	在排序任務中，權值是針對每個分組的，單個樣本的權值沒有意義
silent	是否打印信息
feature_names	特徵名
feature_types	特徵類型
nthread	使用線程數量，默認使用最大線程

Tips：把數據轉化爲二進制可以加快加載速度

dtrain = xgb.DMatrix('train.svm.txt')
dtrain.save_binary('train.buffer')

設置參數

可以通過字典或者列表的方式初始化參數

param = {'max_depth': 2, 'eta': 1, 'silent': 1, 'objective': 'binary:logistic'}
param['nthread'] = 4
param['eval_metric'] = 'auc'

通過字典的方法修改或者增加參數

param['eval_metric'] = ['auc', '[email protected]']

# alternatively:
# plst = param.items()
# plst += [('eval_metric', '[email protected]')]

設置驗證集

evallist = [(dtest, 'eval'), (dtrain, 'train')]

很明顯這裏的param非常重要，那麼它有哪些可以設置的參數呢

增加隨機性

eta 這個就是學習步進，也就是上面中的 $\epsilon$ 。
subsample 這個就是隨機森林的方式，每次不是取出全部樣本，而是有放回地取出部分樣本。有人把這個稱爲行抽取，subsample就表示抽取比例
colsample_bytree和colsample_bylevel 這個是模仿隨機森林的方式，這是特徵抽取。colsample_bytree是每次準備構造一棵新樹時，選取部分特徵來構造，colsample_bytree就是抽取比例。colsample_bylevel表示的是每次分割節點時，抽取特徵的比例。
max_delta_step 這個是構造樹時，允許得到 $f_m(x)$ 的最大值。如果爲0，表示無限制。也是爲了後續構造樹留出空間，和eta相似

控制模型複雜度

max_depth 樹的最大深度
min_child_weight 如果一個節點的權重和小於這玩意，那就不分了
gamma每次分開一個節點後，造成的最小下降的分數。類似於上面的 $Gain$
alpha和lambda就是目標函數裏的表示模型複雜度中的L1範數和L2範數前面的係數

其他參數

booster 表示用哪種模型，一共有gbtree, gbline, dart三種選擇。一般用gbtree。
nthread 並行線成數。如果不設置就是能採用的最大線程。
sketch_eps 這個就是近似算法裏的 $\epsilon$ 。
scale_pos_weight 這個是針對二分類問題時，正負樣例的數量差距過大。

訓練模型

訓練一個模型需要一個參數列表和一個數據集

bst = xgb.train(param, dtrain, num_round = 10, evallist)

保存模型

bst.save_model('0001.model')

模型及其特徵映射也可以轉儲到文本文件中

# dump model
bst.dump_model('dump.raw.txt')
# dump model with feature map
bst.dump_model('dump.raw.txt', 'featmap.txt')

加載模型

bst = xgb.Booster({'nthread': 4})  # init model
bst.load_model('model.bin')  # load data

提前結束

如果你有一個驗證集，你可以使用early_stopping_rounds來找到更合適的boost輪次。提前停止需要至少一組evals（估計指標）。如果有多個，它將使用最後一個。

train(..., evals=evals, early_stopping_rounds=10)

設置該參數後，訓練會持續進行直到驗證集的得分不再增加（評價標準是每early_stopping_rounds輪，驗證集的錯分率不再降低），當提前結束髮生時，模型會額外輸出三個字段：bst.best_score, bst.best_iteration 和 bst.best_ntree_limit，注意返回的模型是最後一次迭代的模型而不是最優的模型。

預測

# 7 entities, each contains 10 features
data = np.random.rand(7, 10)
dtest = xgb.DMatrix(data)
ypred = bst.predict(dtest)

如果生成模型的過程中出發了提前結束，那麼在預測的時候可以設置bst.best_ntree_limit參數來使用最優的模型

ypred = bst.predict(dtest, ntree_limit=bst.best_ntree_limit)

可視化

特徵重要性可視化，需要導入matplotlib模塊

xgb.plot_importance(bst)

指定序號的樹可視化，需要導入graphviz和matplotlib模塊

xgb.plot_tree(bst, num_trees=2)

如果通過ipython運行，可以把指定的樹轉化爲graphviz實例

xgb.to_graphviz(bst, num_trees=2)

實例

參考：https://github.com/dmlc/xgboost/tree/master/demo/guide-python

基礎練習

# import numpy as np
import xgboost as xgb

# 加載數據
dtrain = xgb.DMatrix('file/2/agaricus.txt.train')
dtest = xgb.DMatrix('file/2/agaricus.txt.test')

# 指定參數（樹的深度爲2，學習步進爲1，打印信息，目標函數爲二分類邏輯迴歸）
param = {'max_depth':2, 'eta':1, 'silent':1, 'objective':'binary:logistic'}

# 指定驗證集並開始訓練
watchlist = [(dtest, 'eval'), (dtrain, 'train')]
num_round = 2
bst = xgb.train(param, dtrain, num_round, watchlist)

# 預測並保存模型
preds = bst.predict(dtest)
labels = dtest.get_label()
print('error=%f' % (sum(1 for i in range(len(preds)) if int(preds[i] > 0.5) != labels[i]) / float(len(preds))))
bst.save_model('model/2/0001.model')

# 將模型轉儲爲文本文件
bst.dump_model('model/2/dump.raw.txt')
# 將模型和特徵映射轉儲爲文本文件
# bst.dump_model('model/2/dump.raw.txt', 'model/2/featmap.txt') # 這裏無法生成featmap文件，我們一會再來解決

[11:02:22] 6513x127 matrix with 143286 entries loaded from file/2/agaricus.txt.train
[11:02:22] 1611x127 matrix with 35442 entries loaded from file/2/agaricus.txt.test
[0] eval-error:0.042831 train-error:0.046522
[1] eval-error:0.021726 train-error:0.022263
error=0.021726

import numpy as np

# 把DMatrix格式數據轉化爲二進制格式數據
dtest.save_binary('file/2/dtest.buffer')

# 加載訓練好的模型和轉化好的數據
bst2 = xgb.Booster(model_file='model/2/0001.model')
dtest2 = xgb.DMatrix('file/2/dtest.buffer')

preds2 = bst2.predict(dtest2)

# 查看兩次預測結果是否相同
assert np.sum(np.abs(preds2 - preds)) == 0

[11:02:57] 1611x127 matrix with 35442 entries loaded from file/2/dtest.buffer

import pickle   # 可以序列化對象並保存到磁盤中，用於存儲訓練好的模型

# 序列化模型
pks = pickle.dumps(bst2)

# 加載序列化模型
bst3 = pickle.loads(pks)
preds3 = bst3.predict(dtest2)

assert np.sum(np.abs(preds3 - preds)) == 0

import scipy.sparse

# 通過CSR稀疏矩陣加載數據
print('start running example of build DMatrix from scipy.sparse CSR Matrix')
labels = []
row = []; col = []; dat = []
i = 0

for line in open('file/2/agaricus.txt.train'):
    arr = line.split()
    labels.append(int(arr[0]))
    for item in arr[1:]:
        k,v = item.split(':')
        row.append(i); col.append(int(k)); dat.append(float(v))
    i += 1
    
csr = scipy.sparse.csr_matrix((dat, (row, col)))    # 創建稀疏矩陣
dtrain = xgb.DMatrix(csr, label=labels)    # 加載稀疏矩陣
watchlist = [(dtest, 'eval'), (dtrain, 'train')]    # 指定訓練集，測試集
bst = xgb.train(param, dtrain, num_round, watchlist)    # 訓練模型


# 通過CSC稀疏矩陣加載數據
print('='*30, '\n', 'start running example of build DMatrix from scipy.sparse CSC Matrix')

csc = scipy.sparse.csc_matrix((dat, (row, col)))
dtrain = xgb.DMatrix(csc, label=labels)
watchlist = [(dtest, 'eval'), (dtrain, 'train')]
bst = xgb.train(param, dtrain, num_round, watchlist)

start running example of build DMatrix from scipy.sparse CSR Matrix
[0] eval-error:0.042831 train-error:0.046522
[1] eval-error:0.021726 train-error:0.022263
start running example of build DMatrix from scipy.sparse CSC Matrix
[0] eval-error:0.042831 train-error:0.046522
[1] eval-error:0.021726 train-error:0.022263

# 通過numpy矩陣加載數據
print('start running example of build DMatrix from numpy array')

npymat = csr.todense()
dtrain = xgb.DMatrix(npymat, label=labels)
watchlist = [(dtest, 'eval'), (dtrain, 'train')]
bst = xgb.train(param, dtrain, num_round, watchlist)

start running example of build DMatrix from numpy array
[0] eval-error:0.042831 train-error:0.046522
[1] eval-error:0.021726 train-error:0.022263

定義損失函數和評估指標

注意：

自定義的評估函數需要兩個返回值，分別是評估指標的名稱和評估結果，並且指標名稱中不能含有冒號和空格
注意模型給出的預測是margin，而在自定義損失函數有時候需要對預測進行處理，如用sigmoid函數

import numpy as np
import xgboost as xgb

# 自定義損失函數
print('start running example to used customized objective function')

dtrain = xgb.DMatrix('file/2/agaricus.txt.train')
dtest = xgb.DMatrix('file/2/agaricus.txt.test')

param = {'max_depth': 2, 'eta': 1, 'silent': 1}
watchlist = [(dtest, 'eval'), (dtrain, 'train')]
num_round = 2

start running example to used customized objective function
[13:18:20] 6513x127 matrix with 143286 entries loaded from file/2/agaricus.txt.train
[13:18:20] 1611x127 matrix with 35442 entries loaded from file/2/agaricus.txt.test

我們以Logistic損失爲例

$l(y_i, \hat{y}_i) = y_i\log (1 + e^{-\hat{y}_i}) +(1-y_i)\log (1 + e^{\hat{y}_i})$

於是

$\begin{aligned} g_i & = \frac{ \partial l(y_i,{\hat{y}_i}^{(m-1)}) } {\partial {\hat{y}_i}^{(m-1)} } \\ & = -y_i \frac{e^{-{\hat{y}_i}^{(m-1)}}} {1 + e^{-{\hat{y}_i}^{(m-1)}}} + (1 - y_i) \frac { e^{ {\hat{y}_i}^{(m-1)} } } {1 + e^{{\hat{y}_i}^{(m-1)}}} \\ & = -y_i \left( 1-\frac{1}{1 + e^{-{\hat{y}_i}^{(m-1)}}}\right) + (1 - y_i)\frac{1}{1 + e^{-{\hat{y}_i}^{(m-1)}}} \\ & = \frac{1}{1 + e^{-{\hat{y}_i}^{(m-1)}}} - y_i \\ & = Pred - Label \\ \\ h_i & = \frac{ \partial g_i } {\partial {\hat{y}_i}^{(m-1)} } \\ & = \frac{ e^{-{\hat{y}_i}^{(m-1)}} }{ {(1+e^{-{\hat{y}_i}^{(m-1)}})}^2 } \\ & = \frac{1}{1 + e^{-{\hat{y}_i}^{(m-1)}}} \left( \frac{1}{1 + e^{-{\hat{y}_i}^{(m-1)}}} -1 \right) \\ & = Pred * (Pred - 1) \end{aligned}$ gihi=∂y^i(m−1)∂l(yi,y^i(m−1))=−yi1+e−y^i(m−1)e−y^i(m−1)+(1−yi)1+ey^i(m−1)ey^i(m−1)=−yi(1−1+e−y^i(m−1)1)+(1−yi)1+e−y^i(m−1)1=1+e−y^i(m−1)1−yi=Pred−Label=∂y^i(m−1)∂gi=(1+e−y^i(m−1))2e−y^i(m−1)=1+e−y^i(m−1)1(1+e−y^i(m−1)1−1)=Pred∗(Pred−1)

# 自定義目標函數——對數似然損失，計算g_i,h_i
def logregobj(preds, dtrain):
    labels = dtrain.get_label()
    preds = 1.0 / (1.0 + np.exp(-preds))
    grad = preds - labels
    hess = preds * (1.0 - preds)
    return grad, hess


# 自定義評估函數
# 注意在自定義損失函數時，默認的預測值（Pred）是margin，如logistic迴歸的默認預測值是未經過logistic變換的值
# 所以在損失函數定義中要對Pred做一下siagmoid處理。而內置的評估函數是假定pred經過logistic變換的
def evalerror(preds, dtrain):
    labels = dtrain.get_label()
    # since preds are margin(before logistic transformation, cutoff at 0)
    return 'my-error', float(sum(labels != (preds > 0.0))) / len(labels)

bst = xgb.train(param, dtrain, num_round, watchlist, obj=logregobj, feval=evalerror)    # obj參數傳入目標函數，feval傳入評估函數

[0] eval-rmse:1.59229 train-rmse:1.59597 eval-my-error:0.042831 train-my-error:0.046522
[1] eval-rmse:2.40519 train-rmse:2.40977 eval-my-error:0.021726 train-my-error:0.022263

分步訓練

注意

因爲模型內置的評估函數是假定對預測進行處理的，因爲它需要給我們直接的結果，而在分步訓練中，我們要使用的上一輪未經處理的預測值，所以需要設置參數output_margin=True使我們獲得上一輪未經處理的結果。以logistic損失爲例，我們要的是上一輪的預測得分而不是預測標籤

import numpy as np
import xgboost as xgb

dtrain = xgb.DMatrix('file/2/agaricus.txt.train')
dtest = xgb.DMatrix('file/2/agaricus.txt.test')
watchlist = [(dtest, 'eval'), (dtrain, 'train')]

print ('='*30, '\n', 'start running example to start from a initial prediction')
param = {'max_depth':2, 'eta':1, 'silent':1, 'objective':'binary:logistic'}

# 第一輪訓練
bst = xgb.train(param, dtrain, 1, watchlist)

# 獲取第一輪訓練結果（margin），而不是預測標籤
ptrain = bst.predict(dtrain, output_margin=True)
ptest = bst.predict(dtest, output_margin=True)

# 更新訓練集和測試集的margin
dtrain.set_base_margin(ptrain)
dtest.set_base_margin(ptest)


print('='*30, '\n', 'this is result of running from initial prediction')
bst = xgb.train(param, dtrain, 1, watchlist)

[15:39:34] 6513x127 matrix with 143286 entries loaded from file/2/agaricus.txt.train
[15:39:34] 1611x127 matrix with 35442 entries loaded from file/2/agaricus.txt.test
start running example to start from a initial prediction
[0] eval-error:0.042831 train-error:0.046522
this is result of running from initial prediction
[0] eval-error:0.021726 train-error:0.022263

使用最初生成的n棵樹進行預測

import numpy as np
import xgboost as xgb

dtrain = xgb.DMatrix('file/2/agaricus.txt.train')
dtest = xgb.DMatrix('file/2/agaricus.txt.test')
param = {'max_depth':2, 'eta':1, 'silent':1, 'objective':'binary:logistic'}
watchlist = [(dtest, 'eval'), (dtrain, 'train')]
num_round = 3
bst = xgb.train(param, dtrain, num_round, watchlist)

print('='*30, '\n', 'start testing prediction from first n trees')

# 指定用1棵樹進行預測
label = dtest.get_label()
ypred1 = bst.predict(dtest, ntree_limit=1)

# 指定用2棵樹進行預測
ypred2 = bst.predict(dtest, ntree_limit=2)

# 默認情況下使用所有樹進行預測
ypred3 = bst.predict(dtest)
print('error of ypred from first 1 tree =%f' % (np.sum((ypred1 > 0.5) != label) / float(len(label))))
print('error of ypred from first 2 trees =%f' % (np.sum((ypred2 > 0.5) != label) / float(len(label))))
print('error of ypred from all trees =%f' % (np.sum((ypred3 > 0.5) != label) / float(len(label))))

[15:44:55] 6513x127 matrix with 143286 entries loaded from file/2/agaricus.txt.train
[15:44:55] 1611x127 matrix with 35442 entries loaded from file/2/agaricus.txt.test
[0] eval-error:0.042831 train-error:0.046522
[1] eval-error:0.021726 train-error:0.022263
[2] eval-error:0.006207 train-error:0.007063
start testing prediction from first n trees
error of ypred from first 1 tree =0.042831
error of ypred from first 2 trees =0.021726
error of ypred from all trees =0.006207

使用廣義線性模型作爲基學習器

import xgboost as xgb

dtrain = xgb.DMatrix('file/2/agaricus.txt.train')
dtest = xgb.DMatrix('file/2/agaricus.txt.test')

# change booster to gblinear, so that we are fitting a linear model
# alpha is the L1 regularizer
# lambda is the L2 regularizer
# you can also set lambda_bias which is L2 regularizer on the bias term
param = {'silent':1, 'objective':'binary:logistic', 'booster':'gblinear',
         'alpha': 0.0001, 'lambda': 1}

# normally, you do not need to set eta (step_size)
# XGBoost uses a parallel coordinate descent algorithm (shotgun)（平行座標下降算法）
# there could be affection on convergence（收斂） with parallelization on certain cases
# setting eta to be smaller value, e.g 0.5 can make the optimization more stable
# param['eta'] = 1

watchlist = [(dtest, 'eval'), (dtrain, 'train')]
num_round = 4
bst = xgb.train(param, dtrain, num_round, watchlist)
preds = bst.predict(dtest)
labels = dtest.get_label()
print('error=%f' % (sum(1 for i in range(len(preds)) if int(preds[i] > 0.5) != labels[i]) / float(len(preds))))

[15:56:57] 6513x127 matrix with 143286 entries loaded from file/2/agaricus.txt.train
[15:56:57] 1611x127 matrix with 35442 entries loaded from file/2/agaricus.txt.test
[0] eval-error:0.113594 train-error:0.104099
[1] eval-error:0.117939 train-error:0.105788
[2] eval-error:0.121043 train-error:0.107631
[3] eval-error:0.121043 train-error:0.107477
error=0.121043

交叉驗證

import numpy as np
import xgboost as xgb

dtrain = xgb.DMatrix('file/2/agaricus.txt.train')
param = {'max_depth':2, 'eta':1, 'silent':1, 'objective':'binary:logistic'}
num_round = 2

print('='*30, '\n', 'running cross validation')
# 結果：[輪次] 指標名:指標均值+指標標準差
# 設置5折交叉驗證
xgb.cv(param, dtrain, num_round, nfold=5,
       metrics={'error'}, seed=0,
       callbacks=[xgb.callback.print_evaluation(show_stdv=True)])

print('='*30, '\n', 'running cross validation, disable standard deviation display')
# num_boost_round是最大迭代次數，early_stopping_rounds是提前停止，若最近3輪內指標不再提高，停止並輸出最好輪次，結果不顯示指標的標準差
res = xgb.cv(param, dtrain, num_boost_round=10, nfold=5,
             metrics={'error'}, seed=0,
             callbacks=[xgb.callback.print_evaluation(show_stdv=False),
                        xgb.callback.early_stop(3)])

print(res)


print('='*30, '\n', 'running cross validation, with preprocessing function')
# 定義預處理函數對訓練集、測試集和參數進行處理，可以用來重塑權重
def fpreproc(dtrain, dtest, param):
    # 重塑權重
    label = dtrain.get_label()
    ratio = float(np.sum(label == 0)) / np.sum(label == 1)
    param['scale_pos_weight'] = ratio
    return (dtrain, dtest, param)

# 進行每折驗證前先對數據和參數進行預處理
res = xgb.cv(param, dtrain, num_round, nfold=5,
             metrics={'auc'}, seed=0, fpreproc=fpreproc)

print(res)


# 用自定義的損失函數和評估函數進行交叉驗證
print('='*30, '\n', 'running cross validation, with cutomsized loss function')
def logregobj(preds, dtrain):
    labels = dtrain.get_label()
    preds = 1.0 / (1.0 + np.exp(-preds))
    grad = preds - labels
    hess = preds * (1.0 - preds)
    return grad, hess
def evalerror(preds, dtrain):
    labels = dtrain.get_label()
    return 'error', float(sum(labels != (preds > 0.0))) / len(labels)

param = {'max_depth':2, 'eta':1, 'silent':1}

res = xgb.cv(param, dtrain, num_round, nfold=5, seed=0,
             obj=logregobj, feval=evalerror)

print(res)

[16:14:25] 6513x127 matrix with 143286 entries loaded from file/2/agaricus.txt.train

==============================
running cross validation
[0] train-error:0.0506682+0.009201 test-error:0.0557316+0.0158887
[1] train-error:0.0213034+0.00205561 test-error:0.0211884+0.00365323

==============================
running cross validation, disable standard deviation display
[0] train-error:0.0506682 test-error:0.0557316
Multiple eval metrics have been passed: ‘test-error’ will be used for early stopping.

Will train until test-error hasn’t improved in 3 rounds.
[1] train-error:0.0213034 test-error:0.0211884
[2] train-error:0.0099418 test-error:0.0099786
[3] train-error:0.0141256 test-error:0.0144336
[4] train-error:0.0059878 test-error:0.0062948
[5] train-error:0.0020344 test-error:0.0016886
[6] train-error:0.0012284 test-error:0.001228
[7] train-error:0.0012284 test-error:0.001228
[8] train-error:0.0009212 test-error:0.001228
[9] train-error:0.0006142 test-error:0.001228
Stopping. Best iteration:
[6] train-error:0.0012284+0.000260265 test-error:0.001228+0.00104094

train-error-mean train-error-std test-error-mean test-error-std
0 0.050668 0.009201 0.055732 0.015889
1 0.021303 0.002056 0.021188 0.003653
2 0.009942 0.006076 0.009979 0.004828
3 0.014126 0.001706 0.014434 0.003517
4 0.005988 0.001878 0.006295 0.003123
5 0.002034 0.001470 0.001689 0.000574
6 0.001228 0.000260 0.001228 0.001041

==============================
running cross validation, with preprocessing function
train-auc-mean train-auc-std test-auc-mean test-auc-std
0 0.958228 0.001442 0.958232 0.005778
1 0.981414 0.000647 0.981431 0.002595

==============================
running cross validation, with cutomsized loss function
train-error-mean train-error-std train-rmse-mean train-rmse-std
0 0.050668 0.009201 1.595072 0.003868
1 0.021303 0.002056 2.442600 0.076834

test-error-mean test-error-std test-rmse-mean test-rmse-std
0 0.055732 0.015889 1.598043 0.012826
1 0.021188 0.003653 2.449282 0.080900

預測葉子索引

import xgboost as xgb

dtrain = xgb.DMatrix('file/2/agaricus.txt.train')
dtest = xgb.DMatrix('file/2/agaricus.txt.test')
param = {'max_depth':2, 'eta':1, 'silent':1, 'objective':'binary:logistic'}
watchlist = [(dtest, 'eval'), (dtrain, 'train')]
num_round = 3
bst = xgb.train(param, dtrain, num_round, watchlist)

print ('start testing predict the leaf indices')

# 使用最初的2棵樹
print('='*30, '\n', 'using first 2 trees')
leafindex = bst.predict(dtest, ntree_limit=2, pred_leaf=True)
print(leafindex.shape)
print(leafindex)

# 使用所有樹
print('='*30, '\n', 'using all tree')
leafindex = bst.predict(dtest, pred_leaf=True)
print(leafindex.shape)
print(leafindex)

[16:49:55] 6513x127 matrix with 143286 entries loaded from file/2/agaricus.txt.train
[16:49:55] 1611x127 matrix with 35442 entries loaded from file/2/agaricus.txt.test
[0] eval-error:0.042831 train-error:0.046522
[1] eval-error:0.021726 train-error:0.022263
[2] eval-error:0.006207 train-error:0.007063
start testing predict the leaf indices
using first 2 trees
(1611, 2)
[[4 3]
[3 3]
[4 3]
…
[3 3]
[5 4]
[3 3]]
using all tree
(1611, 3)
[[4 3 5]
[3 3 5]
[4 3 5]
…
[3 3 3]
[5 4 5]
[3 3 3]]

結果顯示1611個測試樣本在每棵樹的葉子索引

利用sklearn中的數據學習xgboost

iris：鳶尾花數據集：多分類（3）數據集
digits：手寫數字數據集：多分類（10）數據集
boston：波士頓房價數據集：迴歸任務數據集

import pickle
import xgboost as xgb

import numpy as np
from sklearn.model_selection import KFold, train_test_split, GridSearchCV
from sklearn.metrics import confusion_matrix, mean_squared_error
from sklearn.datasets import load_iris, load_digits, load_boston

rng = np.random.RandomState(31337)    # 用於產生爲隨機數


# 使用手寫數字數據集中標籤爲0和1的數據做二分類任務
print("Zeros and Ones from the Digits dataset: binary classification")
digits = load_digits(2)
y = digits['target']
X = digits['data']
kf = KFold(n_splits=2, shuffle=True, random_state=rng)    # 交叉驗證
for train_index, test_index in kf.split(X):
    xgb_model = xgb.XGBClassifier().fit(X[train_index], y[train_index])
    predictions = xgb_model.predict(X[test_index])
    actuals = y[test_index]
    print(confusion_matrix(actuals, predictions))

    
# 使用鳶尾花數據集做多分類任務
print('='*30, '\n', "Iris: multiclass classification")
iris = load_iris()
y = iris['target']
X = iris['data']
kf = KFold(n_splits=2, shuffle=True, random_state=rng)
for train_index, test_index in kf.split(X):
    xgb_model = xgb.XGBClassifier().fit(X[train_index], y[train_index])
    predictions = xgb_model.predict(X[test_index])
    actuals = y[test_index]
    print(confusion_matrix(actuals, predictions))


# 使用波士頓房價數據集做迴歸任務
print('='*30, '\n', "Boston Housing: regression")
boston = load_boston()
y = boston['target']
X = boston['data']
kf = KFold(n_splits=2, shuffle=True, random_state=rng)
for train_index, test_index in kf.split(X):
    xgb_model = xgb.XGBRegressor().fit(X[train_index], y[train_index])
    predictions = xgb_model.predict(X[test_index])
    actuals = y[test_index]
    print(mean_squared_error(actuals, predictions))

    
# 參數優化
print('='*30, '\n', "Parameter optimization")
y = boston['target']
X = boston['data']
xgb_model = xgb.XGBRegressor()
clf = GridSearchCV(xgb_model,
                   {'max_depth': [2,4,6],
                    'n_estimators': [50,100,200]}, verbose=1, cv=5)    # 5折交叉驗證，需要對max_depth和n_estimators進行優化，偶爾輸出日誌
clf.fit(X,y)
print(clf.best_score_)
print(clf.best_params_)


# 對sklearn模型進行序列化
print('='*30, '\n', "Pickling sklearn API models")
# 必須用二進制格式打開
pickle.dump(clf, open("best_boston.pkl", "wb"))
clf2 = pickle.load(open("best_boston.pkl", "rb"))
print(np.allclose(clf.predict(X), clf2.predict(X)))


# 設置提前停止
print('='*30, '\n', "With early stopping")
X = digits['data']
y = digits['target']
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=0)
clf = xgb.XGBClassifier()
clf.fit(X_train, y_train, early_stopping_rounds=10, eval_metric="auc",
        eval_set=[(X_test, y_test)])

Zeros and Ones from the Digits dataset: binary classification
[[87 0]
[ 1 92]]
[[91 0]
[ 3 86]]

==============================
Iris: multiclass classification
[[19 0 0]
[ 0 31 3]
[ 0 1 21]]
[[31 0 0]
[ 0 16 0]
[ 0 3 25]]

==============================
Boston Housing: regression
9.860776812557337
15.942418468446029

==============================
Parameter optimization
Fitting 5 folds for each of 9 candidates, totalling 45 fits
[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.
0.6699572097100618
{‘max_depth’: 2, ‘n_estimators’: 100}

==============================
Pickling sklearn API models
True

==============================
With early stopping
[0] validation_0-auc:0.999497
Will train until validation_0-auc hasn’t improved in 10 rounds.
[1] validation_0-auc:0.999497
[2] validation_0-auc:0.999497
[3] validation_0-auc:0.999749
[4] validation_0-auc:0.999749
[5] validation_0-auc:0.999749
[6] validation_0-auc:0.999749
[7] validation_0-auc:0.999749
[8] validation_0-auc:0.999749
[9] validation_0-auc:0.999749
[10] validation_0-auc:1
[11] validation_0-auc:1
[12] validation_0-auc:1
[13] validation_0-auc:1
[14] validation_0-auc:1
[Parallel(n_jobs=1)]: Done 45 out of 45 | elapsed: 8.4s finished

[15] validation_0-auc:1
[16] validation_0-auc:1
[17] validation_0-auc:1
[18] validation_0-auc:1
[19] validation_0-auc:1
[20] validation_0-auc:1
Stopping. Best iteration:
[10] validation_0-auc:1

XGBClassifier(base_score=0.5, booster=‘gbtree’, colsample_bylevel=1,
colsample_bytree=1, gamma=0, learning_rate=0.1, max_delta_step=0,
max_depth=3, min_child_weight=1, missing=None, n_estimators=100,
n_jobs=1, nthread=None, objective=‘binary:logistic’, random_state=0,
reg_alpha=0, reg_lambda=1, scale_pos_weight=1, seed=None,
silent=True, subsample=1)

*使用klearn進行並行

import os

if __name__ == "__main__":
    # NOTE: on posix systems, this *has* to be here and in the
    # `__name__ == "__main__"` clause to run XGBoost in parallel processes
    # using fork, if XGBoost was built with OpenMP support. Otherwise, if you
    # build XGBoost without OpenMP support, you can use fork, which is the
    # default backend for joblib, and omit this.
    try:
        from multiprocessing import set_start_method
    except ImportError:
        raise ImportError("Unable to import multiprocessing.set_start_method."
                          " This example only runs on Python 3.4")
    set_start_method("forkserver")

    import numpy as np
    from sklearn.model_selection import GridSearchCV
    from sklearn.datasets import load_boston
    import xgboost as xgb

    rng = np.random.RandomState(31337)

    print("Parallel Parameter optimization")
    boston = load_boston()

    os.environ["OMP_NUM_THREADS"] = "2"  # or to whatever you want
    y = boston['target']
    X = boston['data']
    xgb_model = xgb.XGBRegressor()
    clf = GridSearchCV(xgb_model, {'max_depth': [2, 4, 6],
                                   'n_estimators': [50, 100, 200]}, verbose=1,
                       n_jobs=2)
    clf.fit(X, y)
    print(clf.best_score_)
    print(clf.best_params_)

使用sklearn評估結果

import xgboost as xgb
import numpy as np
from sklearn.datasets import make_hastie_10_2

X, y = make_hastie_10_2(n_samples=2000, random_state=42)

# Map labels from {-1, 1} to {0, 1}
labels, y = np.unique(y, return_inverse=True)

X_train, X_test = X[:1600], X[1600:]
y_train, y_test = y[:1600], y[1600:]

param_dist = {'objective':'binary:logistic', 'n_estimators':2}

clf = xgb.XGBModel(**param_dist)
# Or you can use: clf = xgb.XGBClassifier(**param_dist)

clf.fit(X_train, y_train,
        eval_set=[(X_train, y_train), (X_test, y_test)], 
        eval_metric='logloss',
        verbose=True)

# Load evals result by calling the evals_result() function
evals_result = clf.evals_result()

print('Access logloss metric directly from validation_0:')
print(evals_result['validation_0']['logloss'])

print('')
print('Access metrics through a loop:')
for e_name, e_mtrs in evals_result.items():
    print('- {}'.format(e_name))
    for e_mtr_name, e_mtr_vals in e_mtrs.items():
        print(' - {}'.format(e_mtr_name))
        print(' - {}'.format(e_mtr_vals))
 

print