AI 助理
备案控制台
开发者社区 人工智能 文章 正文

分类树/装袋法/随机森林算法的R语言实现

2018-09-08 2114
版权
版权声明：
本文内容由阿里云实名注册用户自发贡献，版权归原作者所有，阿里云开发者社区不拥有其著作权，亦不承担相应法律责任。具体规则请查看《 阿里云开发者社区用户服务协议》和 《阿里云开发者社区知识产权保护指引》。如果您发现本社区中有涉嫌抄袭的内容，填写 侵权投诉表单进行举报，一经查实，本社区将立刻删除涉嫌侵权内容。
简介： 介绍了分类树和随机森林的基本知识，并详细介绍了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)] )/2
    minimpurity = 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 = impurity
        splitpoint = cutpoint[i]
      }
    }
  }else {
    splitpoint = xvalues
    minimpurity = 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_obs
    impurity = NULL
    for (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 = xvalues
    minimpurity = 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 = 0
    impurity = 1
  }
  return(c(splitpoint,impurity))
}

splitting = function(subx,suby,split,rf) {
  # subx是一个矩阵，suby是列向量。给定自变量矩阵，返还最优划分变量和相应最优划分条件
  if (!rf) chosen_variable = 1:k
  if (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 = NULL
  index = 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+1
        tree$rightnode[cur] = nodes+2
        nodes = 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 <- tree
  return(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)] )/2
    minimpurity = 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 = impurity
        splitpoint = cutpoint[i]
      }
    }
  }else {
    splitpoint = xvalues
    minimpurity = 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_obs
    impurity = NULL
    for (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 = xvalues
    minimpurity = 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 = 0
    impurity = 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:k
  if (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 = NULL
  index = 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+1
        tree$rightnode[cur] = nodes+2
        nodes = 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 <- tree
  return(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多个小时，把各个变量整合起来尤其困难，系统报错太多的时候心态都有些崩。虽然过程有些辛苦，但总归是写了出来，算是对自己的一次锻炼，写了这篇文章记录一下。

文章标签：
算法
机器学习/深度学习
关键词：
算法分类
算法语言
r语言算法
森林算法
随机森林算法
廖致君
目录
相关文章
java冯坚持
|
27天前
|
存储 人工智能 算法
数据结构与算法细节篇之最短路径问题：Dijkstra和Floyd算法详细描述，java语言实现。
这篇文章详细介绍了Dijkstra和Floyd算法，这两种算法分别用于解决单源和多源最短路径问题，并且提供了Java语言的实现代码。
java冯坚持
62 3
数据结构与算法细节篇之最短路径问题：Dijkstra和Floyd算法详细描述，java语言实现。
土木林森
|
3月前
|
存储 算法 C语言
"揭秘C语言中的王者之树——红黑树：一场数据结构与算法的华丽舞蹈，让你的程序效率飙升，直击性能巅峰！"
【8月更文挑战第20天】红黑树是自平衡二叉查找树，通过旋转和重着色保持平衡，确保高效执行插入、删除和查找操作，时间复杂度为O(log n)。本文介绍红黑树的基本属性、存储结构及其C语言实现。红黑树遵循五项基本规则以保持平衡状态。在C语言中，节点包含数据、颜色、父节点和子节点指针。文章提供了一个示例代码框架，用于创建节点、插入节点并执行必要的修复操作以维护红黑树的特性。
土木林森
99 1
dufadayang
|
3天前
|
算法
树的遍历算法有哪些？
不同的遍历算法适用于不同的应用场景。深度优先搜索常用于搜索、路径查找等问题；广度优先搜索则在图的最短路径、层次相关的问题中较为常用；而二叉搜索树的遍历在数据排序、查找等方面有重要应用。
dufadayang
9 2
java冯坚持
|
30天前
|
存储 算法 关系型数据库
数据结构与算法学习二一：多路查找树、二叉树与B树、2-3树、B+树、B*树。（本章为了解基本知识即可，不做代码学习）
这篇文章主要介绍了多路查找树的基本概念，包括二叉树的局限性、多叉树的优化、B树及其变体（如2-3树、B+树、B*树）的特点和应用，旨在帮助读者理解这些数据结构在文件系统和数据库系统中的重要性和效率。
java冯坚持
18 0
数据结构与算法学习二一：多路查找树、二叉树与B树、2-3树、B+树、B*树。（本章为了解基本知识即可，不做代码学习）
东方睿赢
|
2月前
|
大数据 UED 开发者
实战演练：利用Python的Trie树优化搜索算法，性能飙升不是梦！
在数据密集型应用中，高效搜索算法至关重要。Trie树（前缀树/字典树）通过优化字符串处理和搜索效率成为理想选择。本文通过Python实战演示Trie树构建与应用，显著提升搜索性能。Trie树利用公共前缀减少查询时间，支持快速插入、删除和搜索。以下为简单示例代码，展示如何构建及使用Trie树进行搜索与前缀匹配，适用于自动补全、拼写检查等场景，助力提升应用性能与用户体验。
东方睿赢
50 2
java冯坚持
|
30天前
|
编解码 算法 Java
数据结构与算法学习十九：赫夫曼树树（图文很详细）、赫夫曼编码、应用实践（数据压缩、数据解压）、这一章自我感觉看懂就好。。。
赫夫曼树是一种带权路径长度最短的二叉树，常用于数据压缩和编码。
java冯坚持
31 0
java冯坚持
|
30天前
|
存储 算法
数据结构与算法学习十六：树的知识、二叉树、二叉树的遍历(前序、中序、后序、层次)、二叉树的查找(前序、中序、后序、层次)、二叉树的删除
这篇文章主要介绍了树和二叉树的基础知识，包括树的存储方式、二叉树的定义、遍历方法（前序、中序、后序、层次遍历），以及二叉树的查找和删除操作。
java冯坚持
22 0
何雨晨
|
1月前
|
存储 算法 Java
数据结构和算法--分段树
数据结构和算法--分段树
何雨晨
14 0
糜终
|
1月前
|
算法 安全 Go
Python与Go语言中的哈希算法实现及对比分析
Python与Go语言中的哈希算法实现及对比分析
糜终
35 0
蓝易云
|
3月前
|
机器学习/深度学习 数据采集
R语言逻辑回归、GAM、LDA、KNN、PCA主成分分类分析预测房价及交叉验证
上述介绍仅为简要概述，每个模型在实施时都需要仔细调整与优化。为了实现高度精确的预测，模型选择与调参是至关重要的步骤，并且交叉验证是提升模型稳健性的有效途径。在真实世界的房价预测问题中，可能还需要结合地域经济、市场趋势等宏观因素进行综合分析。
蓝易云
67 3

热门文章

最新文章

  • 1
    算法训练 筛选号码
  • 2
    算法学习之C语言基础
  • 3
    18 回归算法 - Softmax回归 - 信贷审批案例
  • 4
    MapReduce的模式，算法以及用例
  • 5
    数据分析案例-基于随机森林算法的商品评价情感分析
  • 6
    A* + 弗洛伊德 算法 打造 360 向寻路 - [ActionScript3.0]
  • 7
    【编程练习】收集的一些c++代码片，算法排序，读文件，写日志，快速求积分等等
  • 8
    算法学习之顺序结构
  • 9
    字符串相似度算法 递归与动态规划求解分析
  • 10
    【算法研究】Bresenham画线算法
  • 1
    从零到精通：学习这些R语言必学包成为数据分析高手！
    185
  • 2
    R语言读取大型NetCDF文件
    135
  • 3
    R语言——AVOCADO“（异常植被变化检测）算法（1990-2015数据分析）监测森林干扰和再生（含GEE影像下载代码）
    99
  • 4
    R语言Ternary包绘制三元图、RGB三色空间分布图的方法
    146
  • 5
    R语言GD包地理探测器分析时报错、得不到结果等情况的解决方案
    165
  • 6
    R语言raster包批量拼接、融合大量栅格图像
    150
  • 7
    R语言raster包遍历多个文件夹并批量计算每一个文件夹下全部遥感影像的平均值
    85
  • 8
    R语言raster包计算多个栅格图像平均值、标准差的方法
    182
  • 9
    R语言GD包基于栅格图像实现地理探测器与连续参数的自动离散化
    133
  • 10
    R语言geodetector包基于栅格图像实现地理探测器操作
    154

    • 相关课程

    更多
  • 智能创作赛（复赛）：相册应用中的视频故事生成算法介绍
  • 智能创作赛（初赛）：相册应用中的故事生成算法介绍
  • 相册服务中的故事生成算法介绍
  • Go语言核心编程 - 数据结构和算法
  • 大数据之R语言速成与实战
  • 神经网络概览及算法详解

    • 相关电子书

    更多
  • 数据+算法定义新世界
  • 袋鼠云基于实时计算的反黄牛算法
  • Alink：基于Apache Flink的算法平台

    • 相关实验场景

    更多
  • 如何快速训练大模型
  • 使用Swing算法实现商品推荐
  • RSA密码算法设计与实现
  • 欧拉图的构造性证明与算法实现
  • 推荐系统入门之使用ALS算法实现打分预测
    • 下一篇
    无影云桌面