原文首發于簡書于[2018.06.12]

---

本文是我自己動手用R語言寫的實現分類樹的代碼，以及在此基礎上寫的袋裝法（bagging）和隨機森林（random forest）的算法實現。全文的結構是：

• 分類樹

  • 基本知識
  • pred
  • gini
  • splitrule
  • splitrule_best
  • splitrule_random
  • splitting
  • buildTree
  • predict

• 裝袋法與隨機森林

  • 基本知識
  • bagging
  • predict_ensemble
• 性能測試
• 寫在后面

全部的代碼如下：

## x和y為自變量（矩陣）和因變量（列向量） ylevels = levels(y) nlevels = length(ylevels) n = length(y) k = dim(x)[2] pred = function(suby) {# majority vote 所占比例最大的值為預測值vote = rep(0,nlevels)for (i in 1:nlevels) { vote[i] = sum(suby==ylevels[i]) }return(ylevels[which.max(vote)]) }gini = function(suby,summary = FALSE) { # 給定一列因變量，計算它的基尼系數if (!summary) suby = as.vector(summary(suby))obs = sum(suby)return( 1-(drop(crossprod(suby)))/(obs*obs) ) }splitrule = function(subx,suby) { #這里subx和suby都是向量，給定一個自變量，求出它的最優劃分條件subylen = length(suby)xvalues = sort(unique(subx)) if (length(xvalues)>1) { #只有在自變量的取值不是完全相同時（有不同的x），才能夠對x劃分cutpoint = ( xvalues[1:(length(xvalues)-1)] + xvalues[2:length(xvalues)] )/2minimpurity = 1 #初始啟動值for ( i in 1:length(cutpoint) ) {lefty = suby[subx>=cutpoint[i]]righty = suby[subx<cutpoint[i]]impurity = ( gini(lefty)*length(lefty) + gini(righty)*length(righty) )/subylen# 如果一個點是父節點，它的不純度是它的兩個子節點的基尼系數加權和；# 如果是葉節點，則不純度為它本身的基尼系數if (impurity<minimpurity) {minimpurity = impuritysplitpoint = cutpoint[i]}}}else {splitpoint = xvaluesminimpurity = gini(suby)}return(c(splitpoint,minimpurity)) }splitrule_best = function(subx,suby) {# 這是對上面splitrule原有版本的改進。自變量按從小到大排列，每次劃分時，把落入到左邊的觀測記錄下來# 隨著劃分節點不斷增大，落入左節點的數據越來越多，在上一次的基礎上進行累加就行。suby_length = length(suby)xvalues = sort(unique(subx))if (length(xvalues)>1) {intervals = cut(subx,xvalues,right = FALSE) #區間左閉右開suby_splited = matrix(unlist(by(suby,intervals,summary)),ncol = nlevels, byrow = TRUE)suby_splited_left = matrix(apply(suby_splited,2,cumsum),ncol = nlevels)suby_splited_right = sweep(-suby_splited_left,2,as.vector(summary(suby)),FUN = "+") suby_splited_left_obs = apply(suby_splited_left,1,sum)suby_splited_right_obs = suby_length - suby_splited_left_obsimpurity = NULLfor (i in 1:(length(xvalues)-1) ) {impurity = c(impurity,(gini(suby_splited_left[i,],summary = TRUE)*suby_splited_left_obs[i]+ gini(suby_splited_right[i,],summary = TRUE)*suby_splited_right_obs[i])/suby_length)}minimpurity = min(impurity)splitpoint = xvalues[which.min(impurity)]} else {splitpoint = xvaluesminimpurity = gini(suby)}return(c(splitpoint,minimpurity)) }splitrule_random = function(subx,suby) {# splitrule的一個變形。較為極端的方法，不排序，不取唯一值，任意抽取一個數作劃分節點。suby_length = length(suby)subx_withoutmax = subx[subx!=max(subx)]if (length(subx_withoutmax)>0) {splitpoint = subx_withoutmax[sample(length(subx_withoutmax),1)]suby_splited_left = suby[subx<=splitpoint]suby_splited_right = suby[subx>splitpoint]impurity = (gini(suby_splited_left)*length(suby_splited_left) + gini(suby_splited_right)*length(suby_splited_right))/suby_length} else{splitpoint = 0impurity = 1}return(c(splitpoint,impurity)) }splitting = function(subx,suby,split,rf) {# subx是一個矩陣，suby是列向量。給定自變量矩陣，返還最優劃分變量和相應最優劃分條件if (!rf) chosen_variable = 1:kif (rf == TRUE) chosen_variable = sample(1:k,round(sqrt(k)))if (split == "best") temp = apply(subx[,chosen_variable],2,splitrule_best,suby=suby) if (split == "random") temp = apply(subx[,chosen_variable],2,splitrule_random,suby=suby) splitpoint = temp[1,]minimpurity = temp[2,]splitvariable = chosen_variable[which.min(minimpurity)] #確定第幾個變量是最優劃分變量splitpoint = splitpoint[which.min(minimpurity)]minimpurity = min(minimpurity)return(c(splitvariable,splitpoint,minimpurity)) } buildTREE = function(x,y,split = "best",rf = FALSE) {TREE = NULLindex = 1:n tree = as.data.frame(list(leftnode=0,rightnode=0,splitvariable=0,splitpoint=0,obs=n,pred=levels(y)[1],leaf=TRUE,begin=1,end=n)) cur = 1 #當前正在分析的節點編號，它要追趕nodes，每次循環過后增加1。nodes = 1 #當前這棵樹的總節點數量while ((cur<=nodes) ) {beginC = tree$begin[cur]; endC = tree$end[cur]indexCurrent = index[beginC:endC]subx = x[indexCurrent,]suby = y[indexCurrent]impurityCurrent = gini(suby)if (impurityCurrent>0.1) {temp <- splitting(subx,suby,split) if (temp[3]<impurityCurrent) {#如果能進一步劃分為兩個子節點（即不純度能夠降低），nodes增加2# 總會出現無法劃分的情況，即nodes不會增大,這種情況下cur就有可能追趕上nodes了tree$splitvariable[cur] = temp[1]tree$splitpoint[cur] = temp[2]tree$leftnode[cur] = nodes+1tree$rightnode[cur] = nodes+2nodes = nodes +2 tree$leaf[cur] = FALSE #一個節點默認是葉節點，滿足劃分條件時才設為FALSE，這樣寫比較方便。indexL = indexCurrent[subx[,tree$splitvariable[cur]] <= tree$splitpoint[cur]]indexR = indexCurrent[subx[,tree$splitvariable[cur]] > tree$splitpoint[cur]]index[beginC:endC] <- c(indexL,indexR) #這是最為聰明的一步，index被一步一步地精煉更新直至最后成品，搭配begin和end使用。predL = pred(y[indexL])predR = pred(y[indexR]) nodeL = as.data.frame(list(leftnode=0,rightnode=0,splitvariable=0,splitpoint=0,obs=length(indexL),pred=predL,leaf=TRUE,begin = beginC, end = beginC+length(indexL)-1) )nodeR = as.data.frame(list(leftnode=0,rightnode=0,splitvariable=0,splitpoint=0,obs=length(indexR),pred=predR,leaf=TRUE,begin = beginC+length(indexL),end = beginC+length(indexCurrent)-1) )tree = as.data.frame(rbind(tree,nodeL,nodeR))}}cur = cur + 1 #迭代循環iterator}TREE <- treereturn(TREE) }predict = function(TREE,newx) { #newx是一個自變量矩陣，返還一列因變量的預測值。subpredict = function(subx) { #給定一行觀測的自變量，我預測它的因變量row = 1 #從第一個節點開始出發while(TRUE) { if (TREE$leaf[row]==TRUE) {predicted = TREE$pred[row]break}#給定一個觀測值，判斷它是往左走還是往右走，一直走下去，直到抵達葉節點leaf為止if (subx[TREE$splitvariable[row]] <= TREE$splitpoint[row]) {row = TREE$leftnode[row]}else {row = TREE$rightnode[row]}}return(predicted)}return(apply(newx,1,subpredict)) } bagging = function(x,y,trees = 150,split = "best",rf = FALSE) {TREES = list()for (i in 1:trees) {bootstrap_index = sample(1:n,n,replace = TRUE)x_bagging = x[bootstrap_index,]y_bagging = y[bootstrap_index]tree = buildTREE(x_bagging,y_bagging,split,rf)TREES = c(TREES,list(tree)) }return(TREES) }predict_ensemble = function(newx,TREES) { #此處newx為自變量矩陣，TREES為bagging樹林或是隨機森林subpredict_ensemble = function(subx,TREES){ #給定一行觀測的自變量和一片樹林，我預測它的因變量subpredict = function(TREE,subx) { #給定一行觀測的自變量和一棵樹的結構，我預測它的因變量row = 1 #從第一個節點開始出發while(TRUE) { if (TREE$leaf[row]==TRUE) {predicted = TREE$pred[row]break}#給定一個觀測值，判斷它是往左走還是往右走，一直走下去，直到抵達葉節點leaf為止if (subx[TREE$splitvariable[row]] <= TREE$splitpoint[row]) {row = TREE$leftnode[row]}else {row = TREE$rightnode[row]}}return(predicted)}sub_predicted = unlist(lapply(TREES,subpredict,subx)) #對每一行數據進行預測return( names(summary(sub_predicted))[which.max(summary(sub_predicted))] ) #每棵樹進行投票，確定最終預測值}apply(newx,1,subpredict_ensemble,TREES) }

---

分類樹

基本知識

有監督的機器學習中有兩類主要問題，回歸問題（預測房價、預測銷售量）和分類問題（預測性別、預測是否信貸違約）。對應地，決策樹有回歸樹和分類樹兩種。決策樹算是最接近人類思考模式的算法之一，從圖像上來看就是一連串的流程圖，給定一個個體的信息，沿著流程往下走，判斷它最終的歸屬。下面放兩張馬老師的課件，對應的場景是預測某位用戶是否會拖欠貸款。

生成決策樹的規則是不斷地使得不純度下降，對于不純度的計算，可以采用基尼系數，公式為：

例如，十個用戶中有六個拖欠貸款，四個按時還貸，那么不純度為：1 - 0.4^2 - 0.6^2 = 1 - 0.16 - 0.36 = 0.48

給定一個樹節點，是否進行劃分，取決于劃分之后不純度是否會下降，如果下降，則繼續劃分，如果不能下降，該節點成為葉節點。如果能夠劃分，使得不純度下降幅度最大的那個自變量和劃分節點便是最優劃分變量和最優劃分條件。

講完基本概念，我們來看用R代碼實現分類樹，涉及到了下面幾個函數。

• pred：給定一列因變量，求出其中占比例最大的值。分類樹中，對于最終的結果的預測，是由落入到這個葉節點中的占比例最大的值確定的。
• gini：給定一列因變量，求出其基尼系數。
• splitrule：給定一列自變量和因變量，求出這個自變量的最優劃分條件和劃分之后的不純度。
• split_best：對上一版splitrule的改進，減少了計算量。
• split_random：對split_best的再改進，減少了計算量，同時有一定幾率能帶來預測精度的提高（沒有嚴謹的數學證明）。
• splitting：給定一個自變量矩陣和一列因變量，求出所有自變量中最優的劃分變量，以及對應的最優劃分條件和不純度。
• buildTREE：主函數，在上面各個基本函數的基礎上，給定一套數據，生成一棵分類樹。

下面逐個說明。

---

pred()

pred = function(suby) {# majority vote 所占比例最大的值為預測值vote = rep(0,nlevels)for (i in 1:nlevels) { vote[i] = sum(suby==ylevels[i]) }return(ylevels[which.max(vote)]) }

給定一列因變量，求出其中占比例最大的值。分類樹中，對于最終的結果的預測，是由落入到這個葉節點中的占比例最大的值確定的。

---

gini()

gini = function(suby,summary = FALSE) { # 給定一列因變量，計算它的基尼系數if (!summary) suby = as.vector(summary(suby))obs = sum(suby)return( 1-(drop(crossprod(suby)))/(obs*obs) ) }

這個函數就是照著數學公式來寫的，沒什么特別。

summary是R里面一個相當好用的函數，它能夠統計一個factor的頻數：

> example[1] absent absent present absent absent absent absent [8] absent absent present present absent absent absent [15] absent absent absent absent absent absent Levels: absent present > summary(example)absent present 17 3

因為不知道你要傳入的因變量是以像example這樣逐個列出的形式，還是以像summary(example)這樣進行頻數統計的形式，所以設置了一個summary參數，默認FALSE，即默認以前者的形式傳入因變量。

---

splitrule()

splitrule = function(subx,suby) { #這里subx和suby都是向量，給定一個自變量，求出它的最優劃分條件subylen = length(suby)xvalues = sort(unique(subx)) if (length(xvalues)>1) { #只有在自變量的取值不是完全相同時（有不同的x），才能夠對x劃分cutpoint = ( xvalues[1:(length(xvalues)-1)] + xvalues[2:length(xvalues)] )/2minimpurity = 1 #初始啟動值for ( i in 1:length(cutpoint) ) {lefty = suby[subx>=cutpoint[i]]righty = suby[subx<cutpoint[i]]impurity = ( gini(lefty)*length(lefty) + gini(righty)*length(righty) )/subylen# 如果一個點是父節點，它的不純度是它的兩個子節點的基尼系數加權和；# 如果是葉節點，則不純度為它本身的基尼系數if (impurity<minimpurity) {minimpurity = impuritysplitpoint = cutpoint[i]}}}else {splitpoint = xvaluesminimpurity = gini(suby)}return(c(splitpoint,minimpurity)) }

寫得很直觀的一個函數，完全就是按照上面的決策樹的流程圖來寫的。傳入一列因變量和一列自變量，對自變量進行從小到大排序，給定一個劃分點cutpoint，自變量小于等于cutpoint的觀測落入左邊的子節點，自變量大于cutpoint的觀測落入到右邊的子節點，劃分之后，計算不純度并記錄一下。遍歷所有可能的cutpoint，記錄最小的不純度，以及對應的cutpoint值。

---

splitrule_best()

splitrule_best = function(subx,suby) {# 這是對上面splitrule原有版本的改進。自變量按從小到大排列，每次劃分時，把落入到左邊的觀測記錄下來# 隨著劃分節點不斷增大，落入左節點的數據越來越多，在上一次的基礎上進行累加就行。suby_length = length(suby)xvalues = sort(unique(subx))if (length(xvalues)>1) {intervals = cut(subx,xvalues,right = FALSE) #區間左閉右開suby_splited = matrix(unlist(by(suby,intervals,summary)),ncol = nlevels, byrow = TRUE)suby_splited_left = matrix(apply(suby_splited,2,cumsum),ncol = nlevels)suby_splited_right = sweep(-suby_splited_left,2,as.vector(summary(suby)),FUN = "+") suby_splited_left_obs = apply(suby_splited_left,1,sum)suby_splited_right_obs = suby_length - suby_splited_left_obsimpurity = NULLfor (i in 1:(length(xvalues)-1) ) {impurity = c(impurity,(gini(suby_splited_left[i,],summary = TRUE)*suby_splited_left_obs[i]+ gini(suby_splited_right[i,],summary = TRUE)*suby_splited_right_obs[i])/suby_length)}minimpurity = min(impurity)splitpoint = xvalues[which.min(impurity)]} else {splitpoint = xvaluesminimpurity = gini(suby)}return(c(splitpoint,minimpurity)) }

這個是對剛才的splitrule的改進，最終實現的效果是相同的。上一版的缺點在于，每給一個cutpoint，所有的自變量值都要和它進行比較，從而判斷是要落入左邊還是落入右邊。這其實做了很多輪的大小比較，浪費了時間。改進的想法是，把所有的cutpoint都擺出來，把自變量值相應地劃入到不同的區間中，進行記錄。這個時候，R自帶的cut函數就派上用場了。

> example = 1:15 > example[1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 > cut(example,breaks = c(0,3,7,15))[1](0,3] (0,3] (0,3] (3,7] (3,7] (3,7] (3,7] (7,15](7,15](7,15](7,15](7,15](7,15](7,15](7,15] Levels: (0,3](3,7](7,15]

如上，自變量有15個取值，若以3和7作為劃分點，可以把所有的自變量歸入到三個區間內。這樣一來，通過命令`
intervals = cut(subx,xvalues,right = FALSE)
`就能夠把所有的自變量取值進行劃分。

以rpart包中的kyphosis數據集為例，該數據共有81個觀測，其中有64個因變量取值"absent"，有17個因變量取值"present"。

左邊是每個區間內的因變量情況，中間是累積的落入到左邊子節點的情況，右邊是累積的落入到右邊子節點的情況。先看第一行，能夠看到，第一次劃分時，有5個觀測落入到左節點，剩余的76個觀測落入右節點；同樣地，第二行可以看出，有8個觀測落入到左節點，剩下的73個觀測落入到右節點。這樣，只通過一輪的比較大小，就能夠完成任務，成功減少了運算量。

---

splitrule_random()

splitrule_random = function(subx,suby) {# splitrule的一個變形。較為極端的方法，不排序，不取唯一值，任意抽取一個數作劃分節點。suby_length = length(suby)subx_withoutmax = subx[subx!=max(subx)]if (length(subx_withoutmax)>0) {splitpoint = subx_withoutmax[sample(length(subx_withoutmax),1)]suby_splited_left = suby[subx<=splitpoint]suby_splited_right = suby[subx>splitpoint]impurity = (gini(suby_splited_left)*length(suby_splited_left) + gini(suby_splited_right)*length(suby_splited_right))/suby_length} else{splitpoint = 0impurity = 1}return(c(splitpoint,impurity)) }

即便進行了改進，splitrule_best的運算量還是大，因為每次調用這個函數時都要對所有的自變量取值進行排序。splitrule_random函數的做法是，不進行排序，從自變量中隨便挑出一個值作為cutpoint，并記錄不純度。例如，自變量取值是1:20，隨機挑選一個值（例如，選出7），那么自變量取1:7的觀測落入左節點，取8:20的落入右節點。

這是一種相當簡單粗暴的做法，因為計算得到的不純度有可能非常高，而上面的splitrule_best函數計算得到的，是所有可能的不純度里最低的那個，一定是優于splitrule_random得到的不純度的。但是，這種做法有它的道理，具體的會在下面講隨機森林的時候談到。

splitting()

splitting = function(subx,suby,split,rf) {# subx是一個矩陣，suby是列向量。給定自變量矩陣，返還最優劃分變量和相應最優劃分條件if (!rf) chosen_variable = 1:kif (rf == TRUE) chosen_variable = sample(1:k,round(sqrt(k)))if (split == "best") temp = apply(subx[,chosen_variable],2,splitrule_best,suby=suby) if (split == "random") temp = apply(subx[,chosen_variable],2,splitrule_random,suby=suby) splitpoint = temp[1,]minimpurity = temp[2,]splitvariable = chosen_variable[which.min(minimpurity)] #確定第幾個變量是最優劃分變量splitpoint = splitpoint[which.min(minimpurity)]minimpurity = min(minimpurity)return(c(splitvariable,splitpoint,minimpurity)) }

之前的splitrule函數是把一列自變量的最優劃分點給求出來，現在的這個splitting函數是把所有自變量的最優劃分點和相應的不純度都求出來。不純度最低的劃分變量，便是我們要的最優劃分變量。

---

buildTREE()

buildTREE = function(x,y,split = "best",rf = FALSE ) {TREE = NULLindex = 1:n tree = as.data.frame(list(leftnode=0,rightnode=0,splitvariable=0,splitpoint=0,obs=n,pred=levels(y)[1],leaf=TRUE,begin=1,end=n)) cur = 1 #當前正在分析的節點編號，它要追趕nodes，每次循環過后增加1。nodes = 1 #當前這棵樹的總節點數量while ((cur<=nodes) ) {beginC = tree$begin[cur]; endC = tree$end[cur]indexCurrent = index[beginC:endC]subx = x[indexCurrent,]suby = y[indexCurrent]impurityCurrent = gini(suby)if (impurityCurrent>0.1) {temp <- splitting(subx,suby,split) if (temp[3]<impurityCurrent) {#如果能進一步劃分為兩個子節點（即不純度能夠降低），nodes增加2# 總會出現無法劃分的情況，即nodes不會增大,這種情況下cur就有可能追趕上nodes了tree$splitvariable[cur] = temp[1]tree$splitpoint[cur] = temp[2]tree$leftnode[cur] = nodes+1tree$rightnode[cur] = nodes+2nodes = nodes +2 tree$leaf[cur] = FALSE #一個節點默認是葉節點，滿足劃分條件時才設為FALSE，這樣寫比較方便。indexL = indexCurrent[subx[,tree$splitvariable[cur]] <= tree$splitpoint[cur]]indexR = indexCurrent[subx[,tree$splitvariable[cur]] > tree$splitpoint[cur]]index[beginC:endC] <- c(indexL,indexR) #這是最為聰明的一步，index被一步一步地精煉更新直至最后成品，搭配begin和end使用。predL = pred(y[indexL])predR = pred(y[indexR]) nodeL = as.data.frame(list(leftnode=0,rightnode=0,splitvariable=0,splitpoint=0,obs=length(indexL),pred=predL,leaf=TRUE,begin = beginC, end = beginC+length(indexL)-1) )nodeR = as.data.frame(list(leftnode=0,rightnode=0,splitvariable=0,splitpoint=0,obs=length(indexR),pred=predR,leaf=TRUE,begin = beginC+length(indexL),end = beginC+length(indexCurrent)-1) )tree = as.data.frame(rbind(tree,nodeL,nodeR))}}cur = cur + 1 #迭代循環iterator}TREE <- treereturn(TREE) }

這是分類樹算法的主函數，也是最難寫的一個函數。代碼不多，但是寫了很長時間。它的作用是，在上面幾個基本函數的基礎上，把一棵分類樹給生成出來。寫這個函數時，一個首要的問題是：怎么表達一棵樹？

首先映入腦海的是這樣的圖像，一棵樹應該用這樣一幅圖來表達，把劃分變量和劃分條件講得很清楚。可是，R里面怎么用代碼把一幅圖給寫出來？

上了馬老師的課之后，我認識到，分類樹從表面上看是一幅圖，實際上它的結構是一張表（data.frame），圖像只不過是直觀地表達這張表格罷了。buildTREE這個函數，最終生成的便是這樣一個data.frame。

每一行表示一個節點的相關信息，節點的編號由最左端的1,2,3,4...確定。leftnode和rightnode表示這個節點劃分之后的子節點的編號，如果都是0，則表明這個節點是葉節點。splitvariable表示最優劃分變量是第幾個自變量，splitpoint表示劃分點，小于等于這個值的觀測落入左節點，大于這個值的觀測落入右節點。obs表示有多少個觀測落入到這個節點中，leaf表示該節點是否為葉節點，如果是，pred表示這個葉節點的預測值。

這個函數中最關鍵的三行代碼是

indexL = indexCurrent[subx[,tree$splitvariable[cur]] <= tree$splitpoint[cur]] indexR = indexCurrent[subx[,tree$splitvariable[cur]] > tree$splitpoint[cur]] index[beginC:endC] <- c(indexL,indexR)

index是一個十分重要的變量，用來標記落入每一個節點中的觀測編號。還是舉剛才的kyphosis的例子，觀測量有81個，index最初是1:81，在第一次劃分后，有62個觀測落入到左邊，19個觀測落入到右邊，此時對index進行順序調整，使其前62個數字表示左節點的觀測編號，后19個數字表示右節點的觀測編號。按照這種方法進行下去，每輪迭代都對index進行更新，在完成了所有的劃分之后，最終的index記錄了所有的葉節點的觀測編號。

使用index這個向量的意義在于，你不必每輪迭代都使用一個list來記錄落入到某個節點的觀測（我第一次實現分類樹就是這么做的），不必要浪費內存。能夠想到使用向量來記錄觀測，用data.frame來表達一棵決策樹，需要對數據結構有足夠深的理解，這一點我在上一篇博客里提到了。

---

predict()

predict = function(TREE,newx) { #newx是一個自變量矩陣，返還一列因變量的預測值。subpredict = function(subx) { #給定一行觀測的自變量，我預測它的因變量row = 1 #從第一個節點開始出發while(TRUE) { if (TREE$leaf[row]==TRUE) {predicted = TREE$pred[row]break}#給定一個觀測值，判斷它是往左走還是往右走，一直走下去，直到抵達葉節點leaf為止if (subx[TREE$splitvariable[row]] <= TREE$splitpoint[row]) {row = TREE$leftnode[row]}else {row = TREE$rightnode[row]}}return(predicted)}return(apply(newx,1,subpredict)) }

這個函數比較簡單。上一步已經把一棵樹（data.frame）給訓練出來了，現在給定一套數據，通過劃分變量和劃分條件，判斷觀測是往左節點走還是往右節點走，沿著樹一直走下去，直到走到葉節點為止。至此，預測就完成了。

---

裝袋法與隨機森林

基本知識

bagging和random forest都是基于bootstrap的集成學習方法（ensemble learning method），所謂集成學習，指的是把多個機器學習器給整合起來，綜合考慮它們的結果。例如，你訓練了100個學習器去做考試題，有7臺機器選擇A，4臺機器選擇B，80臺機器選擇C，9臺機器選擇D，那么你最后應該選擇C選項，這就是所謂的集體智慧（Wisedom of the Crowd）。

給定原始數據，通過bootstrap生成一套新數據，在這基礎上訓練出一棵樹，再bootstrap得到又一套數據，再訓練出一棵樹，持續進行下去，便得到了一片樹林。bagging和隨機森林便是這樣得到的樹林。之所以要種植多棵分類樹，是因為單棵分類樹有著嚴重的過度擬合問題，在訓練集上表現良好，在測試集上卻做得一塌糊涂。通過綜合多棵樹，能夠大幅度降低variance，而只是小幅度地增加bias，這樣一來測試集的MSE也就能夠顯著下降了。

隨機森林和bagging的區別在于，bagging每一輪迭代都會遍歷所有的自變量，從而找到最優的劃分變量，而隨機森林限定了搜索的自變量數量，在這自變量的子集里尋找最優劃分變量。例如，現有25個自變量，bagging在每一個樹節點處都要進行25輪循環，找到最優的那個，而隨機森林在每個節點處隨機地挑選出5個自變量，在這5個里找到最優的。至于這種做法究竟有什么道理，Hastie在中講得非常生動和清楚。

In other words, in building a random forest, at each split in the tree, the algorithm is not even allowed to consider a majority of the available predictors. This may sound crazy, but it has a clever rationale. Suppose that there is one very strong predictor in the data set, along with a number of other moderately strong predictors. Then in the collection of bagged trees, most or all of the trees will use this strong predictor in the top split. Consequently, all of the bagged trees will look quite similar to each other. Hence the predictions from the bagged trees will be highly correlated. In particular, this means that bagging will not lead to a substantial reduction in variance over a single tree in this setting.

決策樹的決策規則是，每一次劃分節點，都要使用不純度下降最多的劃分變量，一個變量能使不純度下降0.15，另一個變量使不純度下降0.1，決策樹便會選擇前一個變量。這種做法的問題在于，后者雖然效果不好，但是可能在劃分以后，再次劃分能夠使不純度再下降0.2，而前一個變量可能再次劃分時只能使不純度下降0.05。決策樹的問題在于，它能夠保證每一個節點處的劃分都是最有效的，但沒法保證這樣生成的一棵樹就是最好的。這有些像宏觀經濟學里講的動態不一致，每一步的最優和整體上的最優不一致。所以Hastie說，做決策樹的時候，目光要放得長遠一些。

The splitting rule is too short-sighted since a seemingly worthless split early on in the tree might be followed by a very good split — that is, a split that leads to a large reduction in RSS later on.

決策樹會按照最優劃分準則一根筋地生成下去，而隨機森林這一算法的意義在于改變這個最優劃分準側，使得每棵樹之間的相似度下降，隨機森林的隨機便體現在隨機挑選自變量上。

下面介紹一下bagging和隨機森林的實現。

bagging()

bagging = function(x,y,trees = 150,split = "best",rf = FALSE) {TREES = list()for (i in 1:trees) {bootstrap_index = sample(1:n,n,replace = TRUE)x_bagging = x[bootstrap_index,]y_bagging = y[bootstrap_index]tree = buildTREE(x_bagging,y_bagging,split,rf)TREES = c(TREES,list(tree)) }return(TREES) }

這個函數能夠得到一片樹林。做法非常簡單，從原本的數據中隨機抽樣出新數據，生成新的一棵樹就行了，默認是生成150棵樹。bagging和隨機森林都是使用這個函數，通過參數rf就能夠在兩種方法之間進行切換了，之前的splitting函數里也相應地設置了參數rf，其中k是原始數據的自變量個數，bagging遍歷全部k個自變量，而隨機森林只隨機挑選其中sqrt(k)個。

if (!rf) chosen_variable = 1:k if (rf == TRUE) chosen_variable = sample(1:k,round(sqrt(k)))

上面提到的split_random函數，實際上是受到了隨機森林算法的啟發。既然你能隨機地選取劃分變量，那么我為什么不能夠隨機地選取劃分條件呢？同樣是為了降低樹與樹之間的相似性，我這樣做應該也是合理的。只要保證總體上不純度是在下降就行，不必要追求每一次劃分都能帶來大幅度的下降。

---

predict_ensemble()

predict_ensemble = function(newx,TREES) { #此處newx為自變量矩陣，TREES為bagging樹林或是隨機森林subpredict_ensemble = function(subx,TREES){ #給定一行觀測的自變量和一片樹林，我預測它的因變量subpredict = function(TREE,subx) { #給定一行觀測的自變量和一棵樹的結構，我預測它的因變量row = 1 #從第一個節點開始出發while(TRUE) { if (TREE$leaf[row]==TRUE) {predicted = TREE$pred[row]break}#給定一個觀測值，判斷它是往左走還是往右走，一直走下去，直到抵達葉節點leaf為止if (subx[TREE$splitvariable[row]] <= TREE$splitpoint[row]) {row = TREE$leftnode[row]}else {row = TREE$rightnode[row]}}return(predicted)}sub_predicted = unlist(lapply(TREES,subpredict,subx)) #對每一行數據進行預測return( names(summary(sub_predicted))[which.max(summary(sub_predicted))] ) #每棵樹進行投票，確定最終預測值}apply(newx,1,subpredict_ensemble,TREES) }

這個函數和之前的單棵樹的predict函數基本上一樣的，無非是給定數據判斷它向左走還是向右走的問題罷了。唯一的區別是加了個多數投票的函數，對應的是上面舉的例子，100棵樹里，80棵樹選了C選項，所以最后我選C。

---

性能測試

現在試著用mlbench包里LetterRecognition數據集來做一下預測性能測試，數據的預處理過程就不寫了，一個簡單的分層抽樣。

### 單棵分類樹 > tree = buildTREE(train_x,train_y) > predicted1 = predict(tree,test_x) > result1 = mean(predicted1!=test_y)### Bagging > system.time(baggingtrees1 <- bagging(train_x,train_y,split = "best",rf = FALSE))user system elapsed 349.58 0.38 354.38 > system.time(predicted2 <- predict_ensemble(test_x,baggingtrees1))user system elapsed 7.72 0.00 7.80 > result2 = mean(predicted2!=test_y)> system.time(baggingtrees2 <- bagging(train_x,train_y,split = "random",rf = FALSE))user system elapsed 126.17 0.25 129.95 > system.time(predicted3 <- predict_ensemble(test_x,baggingtrees2))user system elapsed 7.15 0.00 7.18 > result3 = mean(predicted3!=test_y)### 隨機森林 > system.time(randomforest1 <- bagging(train_x,train_y,split = "best",rf = TRUE))user system elapsed 155.83 0.06 157.02 > system.time(predicted4 <- predict_ensemble(test_x,randomforest1))user system elapsed 7.84 0.00 7.99 > result4 = mean(predicted4!=test_y)> system.time(randomforest2 <- bagging(train_x,train_y,split = "random",rf = TRUE))user system elapsed 108.47 0.01 109.62 > system.time(predicted5 <- predict_ensemble(test_x,randomforest2))user system elapsed 7.64 0.02 7.66 > result5 = mean(predicted5!=test_y)## ### 結果展示 > result1 [1] 0.38 > result2 [1] 0.215 > result3 [1] 0.235 > result4 [1] 0.23 > result5 [1] 0.215###

從錯誤率來看，單棵分類樹是最差的，bagging和隨機森林差不多，可能是數據量小沒有體現出隨機森林的優勢。使用random的劃分和使用best的劃分在精度上差別不大，但是運算時間明顯減少了，這可以看做是算法改進成功了。

再用R里的randomforest包來試試，發現和我寫的函數性能差別并不大（也許是因為數據量太小了）。

### library(randomForest) bag = randomForest(x = train_x,y = train_y,ntree = 150, mtry = k) predicted6 = predict(bag,newdata = test_x,type = "response") result6 = mean(predicted6!=test_y)rf = randomForest(x = train_x,y = train_y,ntree = 150, mtry = sqrt(k)) predicted7 = predict(rf,newdata = test_x,type = "response") result7 = mean(predicted7!=test_y)> result6 [1] 0.23 > result7 [1] 0.21

---

寫在后面

這份代碼是我入門機器學習以來寫得比較完整的代碼，源代碼來自于馬老師的數據挖掘課，我是全部讀懂了之后憑著理解再自己寫出來的。相當吃力，前前后后大概花了20多個小時，把各個變量整合起來尤其困難，系統報錯太多的時候心態都有些崩。雖然過程有些辛苦，但總歸是寫了出來，算是對自己的一次鍛煉，寫了這篇文章記錄一下。

總結

以上是生活随笔為你收集整理的分类树/装袋法/随机森林算法的R语言实现的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。