如何替换包 randomForest r 中的引导步骤

Posted 2023-03-12

【中文标题】如何替换包 randomForest r 中的引导步骤

【英文标题】：How do I replace the bootstrap step in the package randomForest r

【发布时间】：2015-11-10 12:20:10

【问题描述】：

首先是一些背景信息，这可能在 stats.stackexchange 上更有趣：

在我的数据分析中，我尝试比较不同机器学习方法在时间序列数据（回归，而不是分类）上的表现。例如，我训练了一个 Boosting 训练模型并将其与随机森林训练模型（R 包 randomForest）进行比较。

我使用时间序列数据，其中解释变量是其他数据和因变量的滞后值。

由于某种原因，随机森林的表现严重不佳。我能想到的问题之一是随机森林对每棵树的训练数据执行采样步骤。如果它对时间序列数据执行此操作，则该序列的自回归性质将被完全消除。

为了测试这个想法，我想用所谓的逐块引导步骤替换 randomForest() 函数中的（引导）采样步骤。这基本上意味着我将训练集分成 k 个部分，其中k<<N，其中每个第 k 个部分都按原始顺序排列。如果我对这 k 个部分进行采样，我仍然可以从随机森林中的“随机性”中受益，但时间序列性质基本上保持不变。

现在我的问题是这样的：

为此，我通常会复制现有函数并编辑所需的步骤/行。

randomForest2 <- randomForest()

但 randomForest() 函数似乎是另一个包装器的包装器，用于更深层次的底层函数。那么如何在 randomForest() 函数中编辑实际的引导步骤并仍然定期运行该函数的其余部分？

【问题讨论】：

您可能想要获取原始函数并从那里开始 - 类似于 randomForest <- randomForest:::randomForest.formula。数据有多大？也许您可以通过将滞后值包含为自变量来删除时间序列元素位，然后使用标准方法？

我不确定strata 参数是否能达到你想要的效果，可能不会。否则，重采样都发生在编译后的 C（回归）或 Fortran（分类）代码中，因此要修改它，您需要下载源代码，修改它并重新编译。

如果您尝试使用等于您的训练数据的replace==FALSE 和sampsize 怎么办？这将完全摆脱引导程序，您基本上会得到一组袋装树。

感谢您的建议，我明天会检查并相应更新。

【参考方案1】：

所以对我来说，解决方案不是编辑现有的 randomForest 函数。相反，我使用 Soren H. Welling 给出的 split2 函数自己编写了块级引导程序来创建块。一旦我的数据块被引导，我寻找一个包 (rpart)，它只执行一个回归树并自己聚合它（采取手段）。

就 RMSPE 而言，我的实际数据的结果是比正常随机森林性能略微但持续改进的版本。

对于下面的代码，性能似乎是抛硬币。

以Soren的代码为例，看起来有点像这样：

library(randomForest)
library(doParallel) #parallel package and mclapply is better for linux
library(rpart)

#parallel backend ftw
nCPU = detectCores()
cl = makeCluster(nCPU)
registerDoParallel(cl)

#simulated time series(y) with time roll and lag=1
timepoints=1000;var=6;noise.factor=.2

#past to present orientation    
y = sin((1:timepoints)*pi/30) * 1000 +
  sin((1:timepoints)*pi/40) * 1000 + 1:timepoints
y = y+rnorm(timepoints,sd=sd(y))*noise.factor
plot(y,type="l")

#convert to absolute change, with lag=1
dy = c(0,y[-1]-y[-length(y)]) # c(0,t2-t1,t3-t2,...)

#compute lag 
dy = dy + rnorm(timepoints)*sd(dy)*noise.factor #add noise
dy = c(0,y[-1]-y[-length(y)]) #convert to absolute change, with lag=1 
dX = sapply(1:40,function(i)
  getTheseLags = (1:timepoints) - i
  getTheseLags[getTheseLags<1] = NA #remove before start timePoints
  dx.lag.i = dy[getTheseLags]
)
dX[is.na(dX)]=-100 #quick fix of when lag exceed timeseries
pairs(data.frame(dy,dX[,1:5]),cex=.2)#data structure

#make train- and test-set
train=1:600
dy.train = dy[ train]
dy.test  = dy[-train]
dX.train  = dX[ train,]
dX.test   = dX[-train,]

#classic rf
rf = randomForest(dX.train,dy.train,ntree=500)
print(rf)

#like function split for a vector without mixing
split2 = function(aVector,splits=31) 
  lVector = length(aVector)
  mod = lVector %% splits
  lBlocks = rep(floor(lVector/splits),splits)
  if(mod!=0) lBlocks[1:mod] = lBlocks[1:mod] + 1
  lapply(1:splits,function(i) 
    Stop  = sum(lBlocks[1:i])
    Start = Stop - lBlocks[i] + 1
    aVector[Start:Stop]
  )
  


#create a list of block-wise bootstrapped samples
aBlock <- list()
numTrees <- 500
splits <- 40
for (ttt in 1:numTrees)

  aBlock[[ttt]] <- unlist(
    sample(
      split2(1:nrow(dX.train),splits=splits),
      splits,
      replace=T
    )
  )


#put data into a dataframe so rpart understands it
df1 <- data.frame(dy.train, dX.train)
#perform regression trees for Blocks
rfBlocks = foreach(aBlock = aBlock,
                   .packages=("rpart")) %dopar% 
                     dBlock = df1[aBlock,] 
                     rf = predict( rpart( dy.train ~., data = dBlock, method ="anova" ), newdata=data.frame(dX.test) ) 
                    

#predict test, make results table
#use rowMeans to aggregate the block-wise predictions
results = data.frame(predBlock   = rowMeans(do.call(cbind.data.frame, rfBlocks)),
                     true=dy.test,
                     predBootstrap = predict(rf,newdata=dX.test)
                     )
plot(results[,1:2],xlab="OOB-CV predicted change",
     ylab="trueChange",
     main="black bootstrap and blue block train")
points(results[,3:2],xlab="OOB-CV predicted change",
       ylab="trueChange",
       col="blue")

#prediction results
print(cor(results)^2)


stopCluster(cl)#close cluster

【参考方案2】：

直接改变randomForest的采样(type="regression")：学习基本的C编程，从cran源代码下载randomForest.4.6-10.tar.gz，（如果windows安装Rtools），（如果OSX安装Xcode） ，安装并打开Rstudio，启动新项目，选择包，解压...tar.gz到文件夹，查看src文件夹，打开regrf.c，检查第151和163行。编写新的采样策略，偶尔按Ctrl+Shift+ B 包来重建/编译和覆盖 randomForest 库，更正声明的编译错误，偶尔测试包是否仍然有效，花几个小时找出旧的无信息代码，也许更改描述文件、命名空间文件和一些其他参考，以便包将将名称更改为 randomForestMod，重建，瞧。

下面描述了一种不改变 randomForest 的更简单的方法。任何具有相同特征输入的树都可以使用函数 randomForest::combine 进行修补，因此您可以在纯 R 代码中设计采样机制。我认为这实际上是一个坏主意，但是对于这个非常幼稚的模拟，它实际上具有相似/稍微更好的性能！请记住不要预测绝对目标值，而是要预测相对变化、绝对变化等平稳导数。如果预测绝对值，RF 将退回到预测明天与今天非常接近。这是一个微不足道的无用信息。

编辑代码 [22:42 CEST]

library(randomForest)
library(doParallel) #parallel package and mclapply is better for linux

#parallel backend ftw
nCPU = detectCores()
cl = makeCluster(nCPU)
registerDoParallel(cl)

#simulated time series(y) with time roll and lag=1
timepoints=1000;var=6;noise.factor=.2

#past to present orientation    
y = sin((1:timepoints)*pi/30) * 1000 +
    sin((1:timepoints)*pi/40) * 1000 + 1:timepoints
y = y+rnorm(timepoints,sd=sd(y))*noise.factor
plot(y,type="l")

#convert to absolute change, with lag=1
dy = c(0,y[-1]-y[-length(y)]) # c(0,t2-t1,t3-t2,...)

#compute lag 
dy = dy + rnorm(timepoints)*sd(dy)*noise.factor #add noise
dy = c(0,y[-1]-y[-length(y)]) #convert to absolute change, with lag=1 
dX = sapply(1:40,function(i)
  getTheseLags = (1:timepoints) - i
  getTheseLags[getTheseLags<1] = NA #remove before start timePoints
  dx.lag.i = dy[getTheseLags]
)
dX[is.na(dX)]=-100 #quick fix of when lag exceed timeseries
pairs(data.frame(dy,dX[,1:5]),cex=.2)#data structure

#make train- and test-set
train=1:600
dy.train = dy[ train]
dy.test  = dy[-train]
dX.train  = dX[ train,]
dX.test   = dX[-train,]

#classic rf
rf = randomForest(dX.train,dy.train,ntree=500)
print(rf)

#like function split for a vector without mixing
split2 = function(aVector,splits=31) 
  lVector = length(aVector)
  mod = lVector %% splits
  lBlocks = rep(floor(lVector/splits),splits)
  if(mod!=0) lBlocks[1:mod] = lBlocks[1:mod] + 1
  lapply(1:splits,function(i) 
    Stop  = sum(lBlocks[1:i])
    Start = Stop - lBlocks[i] + 1
    aVector[Start:Stop]
  )
  

nBlocks=10 #combine do not support block of unequal size
rfBlocks = foreach(aBlock = split2(train,splits=nBlocks),
                   .combine=randomForest::combine,
                   .packages=("randomForest")) %dopar% 
                     dXblock = dX.train[aBlock,] ; dyblock = dy.train[aBlock]
                     rf = randomForest(x=dXblock,y=dyblock,sampsize=length(dyblock),
                                       replace=T,ntree=50)
                   
print(rfBlocks)

#predict test, make results table
results = data.frame(predBlock   = predict(rfBlocks,newdata=dX.test),
                     true=dy.test,
                     predBootstrap = predict(rf,newdata=dX.test))
plot(results[,1:2],xlab="OOB-CV predicted change",
     ylab="trueChange",
     main="black bootstrap and blue block train")
points(results[,3:2],xlab="OOB-CV predicted change",
       ylab="trueChange",
       col="blue")

#prediction results
print(cor(results)^2)


stopCluster(cl)#close cluster

【讨论】：

感谢您的两个回答，尽管第二个听起来像是更可行的选择。这是对分块引导程序的不同解释，但我想我的描述不够清楚。实际的 blockwise bootstrap 将数据集划分为 k 个部分，并重新采样这些块，就好像它们是观察值一样。所以在特定的k部分，实际的观察还是按照原来的顺序。

函数 split2 做到这一点

25 个区块中的 600 个训练时间点。第一个块，第一个树，前 24 个时间点。如果您愿意，可以将 replace=FALSE I 设置为。

我不小心按回车发布了，你很快就回复了，所以这是我剩下的评论：你的实现实际上也有优点，所以我试了一下。不幸的是，就我的数据而言，它在 rmspe 方面的表现略差。我能想到的原因是，我的数据不能分成很好的偶数段，所以我必须切断观察。因此包含更少的信息。阅读您的回复后：我将再次尝试使用 replace = F 并回复您，谢谢！另外：感谢并行实现，它为我节省了 1/3 的运行时间。

哦，我从采样池中拉出一个块 dXblock = dX.train[-aBlock,] ; dyblock = dy.train[-aBlock] 我去掉了减号，这样采样池就是一个单独的块： dXblock = dX.train[aBlock,] ; dyblock = dy.train[aBlock]