分類樹/裝袋法/隨機森林算法的R語言實現

原文首發于簡書于[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)
`就能夠把所有的自變量取值進行劃分。

kyphosis數據集劃分情況

以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迭代示意圖

使用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多個小時，把各個變量整合起來尤其困難，系統報錯太多的時候心態都有些崩。雖然過程有些辛苦，但總歸是寫了出來，算是對自己的一次鍛煉，寫了這篇文章記錄一下。

本文來自互聯網用戶投稿，該文觀點僅代表作者本人，不代表本站立場。本站僅提供信息存儲空間服務，不擁有所有權，不承擔相關法律責任。
如若轉載，請注明出處：http://www.pswp.cn/news/389743.shtml
繁體地址，請注明出處：http://hk.pswp.cn/news/389743.shtml
英文地址，請注明出處：http://en.pswp.cn/news/389743.shtml

如若內容造成侵權/違法違規/事實不符，請聯系多彩編程網進行投訴反饋email:809451989@qq.com，一經查實，立即刪除！

相關文章

數據科學學習心得_學習數據科學時如何保持動力

數據科學學習心得_學習數據科學時如何保持動力

數據科學學習心得When trying to learn anything all by yourself, it is easy to lose motivation and get thrown off track.嘗試自己學習所有東西時&#xff0c;很容易失去動力并偏離軌道。 In this article, I will provide you with some tips that I used to stay focus…
閱讀更多...
用php當作cat使用

用php當作cat使用

今天&#xff0c;本來是想敲 node test.js 執行一下&#xff0c;test.js文件&#xff0c;結果 慣性的敲成了 php test.js, 原文輸出了 test.js的內容。 突然覺得&#xff0c;這東西 感覺好像是 cat 命令&#xff0c;嘿嘿&#xff0c;以后要是ubuntu 上沒裝 cat &#xff0c; …
閱讀更多...
建信01. 間隔刪除鏈表結點

建信01. 間隔刪除鏈表結點

建信01. 間隔刪除鏈表結點 給你一個鏈表的頭結點 head&#xff0c;每隔一個結點刪除另一個結點&#xff08;要求保留頭結點&#xff09;。 請返回最終鏈表的頭結點。 示例 1&#xff1a; 輸入&#xff1a;head [1,2,3,4] 輸出: [1,3] 解釋&#xff1a; 藍色結點為刪除的結點…
閱讀更多...
python安裝Crypto:NomodulenamedCrypto.Cipher

python安裝Crypto:NomodulenamedCrypto.Cipher

python 安裝Crypto時出現的錯誤:NomodulenamedCrypto.Cipher 首先直接pip install Crypto 這時會在lib/site-packages/ 文件夾下生成crypto文件夾&#xff0c;將其重命名為Crypto ...然而這個文件夾下沒有Cipher模塊&#xff0c;還需要pip安裝一個pycrypto&#xff0c;不過wind…
閱讀更多...
python多項式回歸_在python中實現多項式回歸

python多項式回歸_在python中實現多項式回歸

python多項式回歸Video Link影片連結 You can view the code used in this Episode here: SampleCode您可以在此處查看 此劇 集中使用的代碼&#xff1a; SampleCode 導入我們的數據 (Importing our Data) The first step is to import our data into python.第一步是將我們的…
閱讀更多...
Uboot 命令是如何被使用的？

Uboot 命令是如何被使用的？

有什么問題請 發郵件至syyxyoutlook.com, 歡迎交流~ 在uboot代碼中命令的模式是這個樣子&#xff1a; 這樣是如何和命令行交互的呢&#xff1f; 在command.h 中, 我們可以看到如下宏定義 將其拆分出來&#xff1a; #define U_BOOT_CMD(name,maxargs,rep,cmd,usage,help) \ U_…
閱讀更多...
2029. 石子游戲 IX

2029. 石子游戲 IX

2029. 石子游戲 IX Alice 和 Bob 再次設計了一款新的石子游戲。現有一行 n 個石子&#xff0c;每個石子都有一個關聯的數字表示它的價值。給你一個整數數組 stones &#xff0c;其中 stones[i] 是第 i 個石子的價值。 Alice 和 Bob 輪流進行自己的回合&#xff0c;Alice 先手…
閱讀更多...
大數據可視化應用_在數據可視化中應用種族平等意識

大數據可視化應用_在數據可視化中應用種族平等意識

大數據可視化應用The following post is a summarized version of the article accepted to the 2020 Visualization for Communication workshop as part of the 2020 IEEE VIS conference to be held in October 2020. The full paper has been published as an OSF Preprint…
閱讀更多...
Windows10電腦系統時間校準

Windows10電腦系統時間校準

有時候新安裝電腦系統&#xff0c;系統時間不對&#xff0c;需要主動去校準系統時間。1、點擊時間 2、日期和時間設置 3、其他日期、時間和區域設置 4、設置時間和日期 5、Internet 時間 6、點擊立即更新&#xff0c;如果更新失敗就查電腦是否已聯網&#xff0c;重試點擊立即更…
閱讀更多...
webpack問題

webpack問題

Cannot find module webpack/lib/node/NodeTemplatePlugin 全局&#xff1a;npm i webpack -g npm link webpack --save-dev 轉載于:https://www.cnblogs.com/doing123/p/8994269.html
閱讀更多...
397. 整數替換

397. 整數替換

397. 整數替換 給定一個正整數 n &#xff0c;你可以做如下操作&#xff1a; 如果 n 是偶數&#xff0c;則用 n / 2替換 n 。 如果 n 是奇數&#xff0c;則可以用 n 1或n - 1替換 n 。 n 變為 1 所需的最小替換次數是多少&#xff1f; 示例 1&#xff1a; 輸入&#xff1a;…
閱讀更多...
pd種知道每個數據的類型_每個數據科學家都應該知道的5個概念

pd種知道每個數據的類型_每個數據科學家都應該知道的5個概念

pd種知道每個數據的類型意見 (Opinion) 目錄 (Table of Contents) Introduction 介紹 Multicollinearity 多重共線性 One-Hot Encoding 一站式編碼 Sampling 采樣 Error Metrics 錯誤指標 Storytelling 評書 Summary 摘要 介紹 (Introduction) I have written about common ski…
閱讀更多...
td

td

單元格td設置padding&#xff0c;而不能設置margin。轉載于:https://www.cnblogs.com/fpcbk/p/9617629.html
閱讀更多...
清除浮動的幾大方法

清除浮動的幾大方法

對于剛接觸到html的一些人經常會用到浮動布局&#xff0c;但對于浮動的使用和清除浮動來說是大為頭痛的&#xff0c;在這里介紹幾個關于清除浮動的的方法。如果你說你要的就是浮動為什么要清除浮動的話&#xff0c;我就真的無言以對了&#xff0c;那你就當我沒說。 關于我們在布…
閱讀更多...
xgboost keras_用catboost lgbm xgboost和keras預測財務交易

xgboost keras_用catboost lgbm xgboost和keras預測財務交易

xgboost kerasThe goal of this challenge is to predict whether a customer will make a transaction (“target” 1) or not (“target” 0). For that, we get a data set of 200 incognito variables and our submission is judged based on the Area Under Receiver Op…
閱讀更多...
2017. 網格游戲

2017. 網格游戲

2017. 網格游戲 給你一個下標從 0 開始的二維數組 grid &#xff0c;數組大小為 2 x n &#xff0c;其中 grid[r][c] 表示矩陣中 (r, c) 位置上的點數。現在有兩個機器人正在矩陣上參與一場游戲。 兩個機器人初始位置都是 (0, 0) &#xff0c;目標位置是 (1, n-1) 。每個機器…
閱讀更多...
HUST軟工1506班第2周作業成績公布

HUST軟工1506班第2周作業成績公布

說明 本次公布的成績對應的作業為&#xff1a; 第2周個人作業&#xff1a;WordCount編碼和測試 如果同學對作業成績存在異議&#xff0c;在成績公布的72小時內&#xff08;截止日期4月26日0點&#xff09;可以進行申訴&#xff0c;方式如下&#xff1a; 畢博平臺的第二周在線答…
閱讀更多...
幣氪共識指數排行榜0910

幣氪共識指數排行榜0910

幣氪量化數據在今天的報告中給出DASH的近期買賣信號&#xff0c;可以看出從今年4月中旬起到目前為止&#xff0c;DASH_USDT的價格總體呈現出下降的趨勢。 轉載于:https://www.cnblogs.com/tokpick/p/9621821.html
閱讀更多...
走出囚徒困境的方法_囚徒困境的一種計算方法

走出囚徒困境的方法_囚徒困境的一種計算方法

走出囚徒困境的方法You and your friend have committed a murder. A few days later, the cops pick the two of you up and put you in two separate interrogation rooms such that you have no communication with each other. You think your life is over, but the polic…
閱讀更多...
2016. 增量元素之間的最大差值

2016. 增量元素之間的最大差值

2016. 增量元素之間的最大差值 給你一個下標從 0 開始的整數數組 nums &#xff0c;該數組的大小為 n &#xff0c;請你計算 nums[j] - nums[i] 能求得的 最大差值 &#xff0c;其中 0 < i < j < n 且 nums[i] < nums[j] 。 返回 最大差值 。如果不存在滿足要求的…
閱讀更多...
最新文章

Copyright @ 2022~2023