寫在前面

AdaBoost是機器學習領域一個很重要很流行的算法，而Julia是一門新興的發展迅速的科學計算語言。本文將從一個實際例子出發，展示如何用Julia語言實現AdaBoost算法。

什麽是AdaBoost

這方面的資料有很多，我將基於Hastie和Tibshirani的ESL(The Elements of Statistical Learning)有關章節的內容，從統計學習的角度簡單介紹一下。另外，我一直在進行ESL的翻譯工作，並試圖實現書中有關算法，歡迎訪問ESL-CN項目主頁，本節的相關翻譯內容參見這裏。

給定預報向量\(X\)，分類器\(G(X)\)在二值\(\\{-1,1\\}\)

\[ \overline{err}=\frac{1}{N}\sum\limits_{i=1}^NI(y_i\neq G(x_i)) \]

在未來預測值上的期望誤差率為\(E_{XY}I(Y\neq G(X))\)

弱分類器是誤差率僅僅比隨機猜測要好一點的分類器。boosting的目的是依次對反復修改的數據應用弱分類器算法，因此得到弱分類器序列\(G_m(x),m=1,2,\ldots,M\) 根據它們得到的預測再通過一個加權來得到最終的預測

\[ G(x)=\mathrm {sign}(\sum\limits_{m=1}^M\alpha_mG_m(x)) \]

用一個概念圖（圖來自ESL原書）表示如下：

具體來說，對每步boosting的數據修改是對每個訓練觀測\((x_i,y_i),i=1,2,\ldots,N\)賦予權重\(w_1,w_2,\ldots,w_N\)。初始化所有的權重設為\(w_i=1/N\)，使得第一步以通常的方式對數據進行訓練分類器。對每個接下來的叠代\(m=2,3,\ldots,M\)，單獨修改觀測的權重，然後將分類算法重新應用到加權觀測值上。在第\(m\)步，上一步中被分類器\(G_{m-1}(x)\)的誤分類的觀測值增大了權重，而正確分類的觀測值權重降低了。因此當叠代繼續，很難正確分類的觀測受到越來越大的影響。每個相繼的分類器因此被強制集中在上一步誤分類的訓練數據上。

算法10.1顯示了AdaBoost.M1算法的詳細細節。當前的分類器\(G_m(x)\)由第2(a)行的加權觀測值得到。在第2(b)行計算加權誤差率。第2(c)行計算賦予\(G_m(x)\)的權重\(\alpha_m\)來得到最終的分類器\(G(x)\)(第3行)。每個觀測的個體權重在第2(d)行進行更新。在導出序列中下一個分類器\(G_{m+1}(x)\)時，被分類器\(G(x)\)錯誤分類的觀測值的權重被因子\(exp(\alpha_m)\)進行縮放以提高它們的相對影響。

例子

特征\(X_1,\ldots,X_{10}\)是標準獨立高斯分布，目標\(Y\)定義如下
\[ Y= \left\{ \begin{array}{ll} 1&\text{if } \sum_{j=1}^{10}X_j^2>\chi_{10}^2(0.5)\-1 & \text{otherwise} \end{array} \right. \qquad (10.2) \]

這裏\(\chi_{10}^2(0.5)=9.34\)是自由度為10的卡方隨機變量的中位數（10個標準的高斯分布的平方和）。有2000個訓練情形，每個類別大概有1000個情形，以及10000個測試觀測值。這裏我們取稱為“stump”的弱分類器：含兩個終止結點的分類樹。

實現

Julia的具體細節參見官方manual。

首先我們定義模型的結構，我們需要兩個參數，弱分類器的個數n_clf和存儲n_clf個弱分類器的n_clf\(\times 4\)的矩陣。因為對於每個弱分類器——兩個終止結點的stump，我們需要三個參數確定，分割變量的編號idx，該分割變量對應的cutpoint值val，以及分類的方向flag（當flag取1時則所有比cutpoint大的觀測值分到樹的右結點，而flag取0時分到左結點），另外算法中需要確定的alpha參數，所以一個stump需要四個參數。下面代碼默認弱分類器個數為10。

struct Adaboost
    n_clf::Int64
    clf::Matrix
end

function Adaboost(;n_clf::Int64 = 10)
    clf = zeros(n_clf, 4)
    return Adaboost(n_clf, clf)
end

訓練模型

function train!(model::Adaboost, X::Matrix, y::Vector)
    n_sample, n_feature = size(X)
    ## initialize weight
    w = ones(n_sample) / n_sample
    threshold = 0
    ## indicate the classification direction
    ## consider observation obs which is larger than cutpoint.val
    ## if flag = 1, then classify obs as 1
    ## else if flag = -1, classify obs as -1
    flag = 0
    feature_index = 0
    alpha = 0
    for i = 1:model.n_clf
        ## step 2(a): stump
        err_max = 1e10
        for feature_ind = 1:n_feature
            for threshold_ind = 1:n_sample
                flag_ = 1
                err = 0
                threshold_ = X[threshold_ind, feature_ind]

                for sample_ind = 1:n_sample
                    pred = 1
                    x = X[sample_ind, feature_ind]
                    if x < threshold_
                        pred = -1
                    end
                    err += w[sample_ind] * (y[sample_ind] != pred)
                end
                err = err / sum(w)
                if err > 0.5
                    err = 1 - err
                    flag_ = -1
                end

                if err < err_max
                    err_max = err
                    threshold = threshold_
                    flag = flag_
                    feature_index = feature_ind
                end
            end
        end
        ## step 2(c)
        #alpha = 1/2 * log((1-err_max)/(err_max))
        alpha = 1/2 * log((1.000001-err_max)/(err_max+0.000001))
        ## step 2(d)
        for j = 1:n_sample
            pred = 1
            x = X[j, feature_index]
            if flag * x < flag * threshold
                pred = -1
            end
            w[j] = w[j] * exp(-alpha * y[j] * pred)
        end
        model.clf[i, :] = [feature_index, threshold, flag, alpha]
    end
end

預測模型

function predict(model::Adaboost,
                 x::Matrix)
    n = size(x,1)
    res = zeros(n)
    for i = 1:n
        res[i] = predict(model, x[i,:])
    end
    return res
end

function predict(model::Adaboost,
                 x::Vector)
    s = 0
    for i = 1:model.n_clf
        pred = 1
        feature_index = trunc(Int64,model.clf[i, 1])
        threshold = model.clf[i, 2]
        flag = model.clf[i, 3]
        alpha = model.clf[i, 4]
        x_temp = x[feature_index]
        if flag * x_temp < flag * threshold
            pred = -1
        end
        s += alpha * pred
    end

    return sign(s)

end

接下來應用到模擬例子中

function generate_data(N)
    p = 10
    x = randn(N, p)
    x2 = x.*x
    c = 9.341818 #qchisq(0.5, 10)
    y = zeros(Int64,N)
    for i=1:N
        tmp = sum(x2[i,:])
        if tmp > c
            y[i] = 1
        else
            y[i] = -1
        end
    end
    return x,y
end

function test_Adaboost()
    x_train, y_train = generate_data(2000)
    x_test, y_test = generate_data(10000)
    m = 1:20:400
    res = zeros(size(m, 1))
    for i=1:size(m, 1)
        model = Adaboost(n_clf=m[i])
        train!(model, x_train, y_train)
        predictions = predict(model, x_test)
        println("The number of week classifiers ", m[i])
        res[i] = classification_error(y_test, predictions)
        println("classification error: ", res[i])
    end
    return hcat(m, res)
end

作出誤差隨叠代次數的圖象如下

完整代碼參見這裏覺得項目很好的話記得star鼓勵一下哦

AdaBoost對實際數據分類的Julia實現