随机森林置信区间和预测

机器算法验证 随机森林 预测区间 覆盖概率
2022-03-17 13:11:32

这是一个简短的模拟,用于检查论文中介绍的随机森林置信区间的覆盖率,当用作预测区间时

S. Wager、T. Hastie 和 B. Efron。随机森林的置信区间:折刀和无穷小折刀机器学习研究杂志,15,第 1625-1651 页(2014 年)

数据是根据第 4.3 节中描述的过程模拟的

J.弗里德曼。多元自适应回归样条统计年鉴。19(1), pp 1-67 (1991)

有六个独立的预测变量,每个都具有分布,以及一个具有分布的独立随机变量响应变量定义为 X1,,X6U[0,1]ϵN(0,1)

Y=10sin(πX1X2)+20(X31/2)2+10X4+5X5+ϵ.

friedman <- function(n, p = 6) {
    X <- matrix(runif(n*p), nrow = n, ncol = p, dimnames = list(1:n, paste0("x_", 1:p)))
    y <- 10*sin(pi*X[, 1]*X[, 2]) + 20*(X[, 3] - 0.5)^2 + 10*X[, 4] + 5*X[, 5] + rnorm(n)
    data.frame(cbind(y, X))
}

我们生成大小为 1.000 的训练样本和大小为 100.000 的测试样本。

set.seed(42)

n <- 10^3
training <- friedman(n)

n_tst <- 10^5
test <- friedman(n_tst)

library(ranger)

rf <- ranger(y ~ ., data = training, num.trees = 10^3, keep.inbag = TRUE)

pred <- predict(rf, data = test, type = "se", se.method = "infjack")

y_hat <- pred$predictions

Lower <- y_hat - 1.96 * pred$se
Upper <- y_hat + 1.96 * pred$se

mean((Lower <= test$y) & (test$y <= Upper))

我们设定了95%的标称覆盖率,但模拟得出的覆盖率约为20.2%问题在于该模拟中获得的低有效预测覆盖率。

笔记:

我们正在计算级别的置信区间,如脚注(第 3 页)中所述:1αy^±zασ^

正如 usεr11852 在下面的评论中指出的那样,有效覆盖率随着我们增加森林中树木的数量而减少。例子:

num.trees =  50 => effective coverage = 0.94855
num.trees = 100 => effective coverage = 0.76876 
num.trees = 150 => effective coverage = 0.68959 
num.trees = 200 => effective coverage = 0.56038 
num.trees = 250 => effective coverage = 0.55393 
num.trees = 300 => effective coverage = 0.32304 
num.trees = 350 => effective coverage = 0.55413 
num.trees = 400 => effective coverage = 0.26372 
num.trees = 450 => effective coverage = 0.26232 
num.trees = 500 => effective coverage = 0.23139 
2个回答

OP 的编写方式和结果的评估方式,这个问题似乎混淆了预测区间和置信区间。

预测区间给出随机变量的区间估计,而置信区间提供参数的区间估计。两者之间的差异在概念上与变量的标准偏差与其平均值的标准误差之间的差异一样大。

编辑(关闭逻辑差距):混淆通常与模型预测有关,因为它们既可以被视为单个随机变量的猜测,也可以被视为条件均值的估计。结束编辑。

Jackknife 估计为估计的参数提供了不确定性,因此您尝试的方法为条件均值提供了置信区间,而您的 OP 要求提供预测区间。

这可以在一定程度上解释为什么更多的树(-> 条件均值的更稳健估计)导致更短的间隔。

我能想到的两种方法:

分位数随机森林

我想给你一个例子,如何使用分位数随机森林来产生(概念上有点太窄)预测区间,但我最终得到了 90% 的覆盖率,而不是 80% 的覆盖率,另见@Andy W 的回答和@Zen 的评论. 不同的参数化也会发生类似的情况。我稍微清理了代码,所以无论如何,这就是它。

library(ranger)

friedman <- function(n) {
  X <- matrix(runif(n * 6), nrow = n)
  data.frame(
    X,
    y = 10 * sin(pi * X[, 1] * X[, 2]) + 20 * (X[, 3] - 0.5)^2 +
        10 * X[, 4] + 5 * X[, 5] + rnorm(n)
  )
}

set.seed(42)

n <- 10^5
training <- friedman(n)

n_tst <- 10^5
test <- friedman(n_tst)

rf <- ranger(y ~ ., data = training, quantreg = TRUE)

pred <- predict(rf, data = test, quantiles = c(0.1, 0.9), type = "quantiles")

y_hat <- pred$predictions

Lower <- y_hat[, 1]
Upper <- y_hat[, 2]

mean((Lower <= test$y) & (test$y <= Upper))  # 0.91616

通用预测区间

我不时使用一种通用方法,它适用于任何回归方法。它基于两个模型:第一个是条件均值模型。第二个模型对模型的绝对残差进行建模,并因此为条件响应分布的条件标准差提供模型。然后以与您在示例中相同的方式构建预测区间。与上面的分位数方法类似,它忽略了模型的不确定性。第二个问题是第二个模型中使用的损失函数相当随意。

library(ranger)

friedman <- function(n) {
  X <- matrix(runif(n * 6), nrow = n)
  data.frame(
    X,
    y = 10 * sin(pi * X[, 1] * X[, 2]) + 20 * (X[, 3] - 0.5)^2 +
      10 * X[, 4] + 5 * X[, 5] + rnorm(n)
  )
}

set.seed(42)

n <- 10^5
training <- friedman(n)

n_tst <- 10^5
test <- friedman(n_tst)

rf_mean <- ranger(y ~ ., data = training)
rf_sd <- ranger(abs(rf_mean$predictions - training$y) ~ ., data = training)

pred_mean <- predict(rf_mean, data = test)$predictions
pred_sd <- predict(rf_sd, data = test)$predictions

Lower <- pred_mean - 1.96 * pred_sd
Upper <- pred_mean + 1.96 * pred_sd

mean((Lower <= test$y) & (test$y <= Upper))  # 89%

有趣的是,这种方法的覆盖范围与分位数随机森林的一种非常相似,这可能在某种程度上表明弗里德曼的数据很容易预测?

不是一个直接的答案(除非这里混合了置信区间和预测区间),但是从您的示例中的各个树中提取分位数确实会为个人级别预测产生近似正确的覆盖率(实际上比预期的覆盖率要好一些)间隔。

library(ranger)

friedman <- function(n, p = 6) {
    X <- matrix(runif(n*p), nrow = n, ncol = p)
    y <- 10*sin(pi*X[, 1]*X[, 2]) + 20*(X[, 3] - 0.5)^2 + 10*X[, 4] + 5*X[, 5] + rnorm(n)
    as.data.frame(cbind(y, X), col.names = c("y", sapply(1:p, function(i) paste0("x_", i))))
}

set.seed(10)

# Setting smaller to be a smidge
# more memory safe on my machine
n <- 10^3
training <- friedman(n)

n_tst <- 10^4
test <- friedman(n_tst)

rf <- ranger(y ~ ., data = training, num.trees = 500, keep.inbag = FALSE)

# Generating prediction intervals from trees directly seems to 
# Work just fine
pred <- predict(rf, data = test, predict.all = TRUE)
lower <- apply(pred$predictions,1,quantile,0.05) #90% interval
upper <- apply(pred$predictions,1,quantile,0.95)
mean((lower <= test$y) & (test$y <= upper)) 
# [1] 0.9655 on my machine, expected 90%

您可以通过随机森林免费获得这些。对于非常极端的分位数,需要增加树的数量,但在我跨项目的经验中,这些通常很难很好地覆盖真实数据示例中的极端百分位数。