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

Question

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

rfunctioneditrandom-foreststatistics-bootstrap

5

首先提供一些背景信息，这可能更有趣的统计堆栈交换：

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

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

出于某种原因，随机森林的性能严重不足。我能想到的问题之一是，随机森林对每个树的训练数据执行抽样步骤。如果它对时间序列数据执行此操作，则系列的自回归性质完全被移除。

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

现在我的问题是：

要实现这一点，我通常会复制现有函数并编辑所需的步骤/行。

randomForest2 <- randomForest()

但是randomForest()函数似乎是另一个更深层次底层函数的包装器。那么我如何编辑randomForest()函数中实际的bootstrap步骤并仍然正常运行函数的其余部分？

- DaReal

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

1

我不确定strata参数是否能实现你想要的，很可能不行。否则，所有重抽样都发生在编译后的C代码（回归）或Fortran代码（分类）中，因此如果要修改它，您需要下载源代码，对其进行修改并重新编译。 - joran

如果您尝试使用replace==FALSE和sampsize等于您的训练数据，会怎样呢？这将完全消除自助法，并且您基本上会得到一组装袋树。 - Tchotchke

谢谢您的建议，我明天会查看并相应更新。 - DaReal

2个回答

1

直接修改randomForest(type="reggression")的抽样方法：学习基本的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

- Soren Havelund Welling

谢谢你们两个的答案，虽然第二个听起来更可行一些。这是块状自助法的不同解释，但我想我的描述不够清晰。实际的块状自助法是将数据集分成k个部分，并重新采样这些块，就好像它们是观察值一样。因此，在特定的第k部分，实际观测值仍按原始顺序排列。 - DaReal

函数split2实现了这个功能。 - Soren Havelund Welling

600个火车时刻点分为25个块。第一个块，第一棵树，前24个时间点。如果您喜欢，可以将replace=FALSE I替换。 - Soren Havelund Welling

哦，我正在从样本池中取出一个块。dXblock = dX.train[-aBlock,] ; dyblock = dy.train[-aBlock]我去掉了负号，这样采样池就只有一个块：dXblock = dX.train[aBlock,] ; dyblock = dy.train[aBlock] - Soren Havelund Welling

我会更正我的答案。但是无论如何，我认为这种方法不会大大节省您的模型时间。祝你好运 :) - Soren Havelund Welling

显示剩余2条评论

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- DaReal · Accepted Answer

因此，对我来说解决方案并不是编辑现有的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