tidymodels交叉验证教程
library("tidyverse")
library("tidymodels")
library("parsnip")
library("brotools")
library("mlbench")
#随机森林预测房价
data("BostonHousing2")
head(BostonHousing2)
boston <- BostonHousing2 %>%
select(-medv, -town, -lon, -lat) %>%
rename(price = cmedv)
#删除了town ， lat和lon因为列栏中包含的信息就足够了。
#为了训练和评估模型的性能，数据分成两部分。 训练集分为两个训练集,测试集最后使用。
train_test_split <- initial_split(boston, prop = 0.9)
housing_train <- training(train_test_split)
housing_test <- testing(train_test_split)
#我想训练一个随机森林来预测房屋的价格，但随机森林有所谓的超参数，这些参数是无法从数据中估算或学习的参数。
#相反，这些参数必须由分析师选择。 为了选择它们，您可以使用似乎运行良好的文献中的值（如在宏计量经济学中完成）或者您可以进一步将列车集拆分为两个，创建超参数网格，在一个上训练模型网格的所有值的部分数据，并比较数据的第二部分上的模型的预测。
#然后，您将坚持使用最佳性能的模型，例如，具有最低RMSE的模型。 问题是，您无法仅使用一个值来估计RMSE的真实值。 这就像你想通过从人口中抽取一个单一的观察来估计人口的高度。
#为了给出一组超参数的RMSE的真实值，而不是做一次拆分，我会做30.然后计算平均RMSE，这意味着为超参数的每个值的组合训练30个模型我有兴趣。
#先让我们使用{rsample}包中的mc_cv()函数再次拆分训练数据。 此函数实现蒙特卡罗交叉验证：
validation_data <- mc_cv(housing_train, prop = 0.9, times = 30)
validation_data$splits[[1]]
#数据预处理
#我将分三步完成。 第一步和第二步用于标准化，第三步将字符和因子变量转换为虚拟变量。 这是必要的，因为随机森林不能直接处理因子变量。 让我们定义一个配方来做到这一点，并从预处理测试集开始。
#包装函数，可应用于各种数据集：
simple_recipe <- function(dataset){
recipe(price ~ ., data = dataset) %>%
step_center(all_numeric()) %>%
step_scale(all_numeric()) %>%
step_dummy(all_nominal())
}
#使用prep()函数，该函数从数据中估算处理数据所需的参数。 例如，对于居中， prep()估计平均值，然后从变量中减去平均值。
#使用bake() ，将估算值应用于数据：
testing_rec <- prep(simple_recipe(housing_test), testing = housing_test)
test_data <- bake(testing_rec, new_data = housing_test)
#在使用prep()和bake()之前分割数据很重要，因为如果没有，您将使用prep()步骤中测试集的观察结果，从而将测试集中的知识引入训练数据中。 这称为数据泄漏，必须避免。 这就是为什么有必要首先将训练数据分成分析和评估集，然后分别预处理这些集。
#但是， validation_data对象现在不能与recipe()一起使用，因为它不是数据帧。 不用担心，我只需编写一个函数，从validation_data对象中提取分析和评估集，应用预处理，训练模型，并返回RMSE。 这将是一项重要功能，是分析的核心。
#在此之前，让我们运行一个简单的线性回归作为基准。 对于线性回归，我不会使用任何CV，所以让我们预处理训练集：
trainlm_rec <- prep(simple_recipe(housing_train), testing = housing_train)
trainlm_data <- bake(trainlm_rec, new_data = housing_train)
linreg_model <- lm(price ~ ., data = trainlm_data)
broom::augment(linreg_model, newdata = test_data) %>%
rmse(price, .fitted)
#broom::augment()将预测添加到新列的.fitted 。 我不会在随机森林中使用这个技巧，因为我将使用{ranger}随机森林没有augment()方法。 我会自己将预测添加到数据中。
#现在让我们回到随机森林并编写一个大函数
my_rf <- function(mtry, trees, split, id){
analysis_set <- analysis(split)
analysis_prep <- prep(simple_recipe(analysis_set), training = analysis_set)
analysis_processed <- bake(analysis_prep, new_data = analysis_set)
model <- rand_forest(mtry = mtry, trees = trees) %>%
set_engine("ranger", importance = 'impurity') %>%
fit(price ~ ., data = analysis_processed)
assessment_set <- assessment(split)
assessment_prep <- prep(simple_recipe(assessment_set), testing = assessment_set)
assessment_processed <- bake(assessment_prep, new_data = assessment_set)
tibble::tibble("id" = id,
"truth" = assessment_processed$price,
"prediction" = unlist(predict(model, new_data = assessment_processed)))
}
#rand_forest()函数可从{parsnip}包中获得。 该软件包为许多其他机器学习包提供了统一的接口。
#这意味着不必学习range()和randomForest()的语法，而且......你可以简单地使用rand_forest()函数并将engine参数更改为你想要的randomForest（ranger， randomForest等）。
#试试这个功能
results_example <- map2_df(.x = validation_data$splits,
.y = validation_data$id,
~my_rf(mtry = 3, trees = 200, split = .x, id = .y))
head(results_example)
#我现在可以在mtry = 3和trees = 200时计算RMSE：
results_example %>%
group_by(id) %>%
rmse(truth, prediction) %>%
summarise(mean_rmse = mean(.estimate)) %>%
pull
#随机森林的RMSE已经低于线性回归。 现在的目标是通过调整mtry和trees超参数来降低此RMSE。 为此，我将使用{mlrMBO}包中实现的贝叶斯优化方法。
#贝叶斯超参数优化
#重新使用上面的代码，并定义一个函数，该函数执行从预处理到返回我希望通过调整超参数RMSE来最小化的度量的所有内容：
tuning <- function(param, validation_data){
mtry <- param[1]
trees <- param[2]
results <- purrr::map2_df(.x = validation_data$splits,
.y = validation_data$id,
~my_rf(mtry = mtry, trees = trees, split = .x, id = .y))
results %>%
group_by(id) %>%
rmse(truth, prediction) %>%
summarise(mean_rmse = mean(.estimate)) %>%
pull
}
tuning(c(3, 200), validation_data)
plot_points <- crossing("mtry" = 3, "trees" = seq(200, 300))
#绘制mtry = 3的RMSE值和200到300的trees 。因为我需要100次评估这个昂贵的函数。
# 如果评估函数，我可以通过改变mtry值来制作3D图，但是如果评估函数的价格便宜，我会进行详尽的网格搜索以找到超参数而不是使用贝叶斯优化。
plot_data <- plot_points %>%
mutate(value = map_dbl(seq(200, 300), ~tuning(c(3, .), validation_data)))
plot_data %>%
ggplot(aes(y = value, x = trees)) +
geom_line(colour = "#82518c") +
theme_blog() +
ggtitle("RMSE for mtry = 3")
#对于mtry = 3 ，最小值似乎在255左右。最小化的函数根本不是平滑的。
#我现在按照arxiv文章中的代码来运行优化。 但据我所知，一个更简单的模型，称为代理模型，用于寻找有希望的点，并在这些点评估函数的价值。 这似乎（在精神上）与Gourieroux，Monfort，Renault中描述的间接推理方法有些相似。
#让我们首先加载包并创建要优化的函数：
library("mlrMBO")
fn <- makeSingleObjectiveFunction(name = "tuning",
fn = tuning,
par.set = makeParamSet(makeIntegerParam("x1", lower = 3, upper = 8),
makeIntegerParam("x2", lower = 50, upper = 500)))
#此功能基于我之前定义的功能。 要优化的参数也被定义为它们的界限。 我会在3到8之间寻找mtry ，在50到500之间寻找trees 。
#我认为这意味着这10个点是用于启动整个过程的点。 我不明白为什么他们必须从超立方体中采样，但还可以。 然后我选择代理模型，随机森林，并预测标准误差。
library(lhs)# for randomLHS
des <- generateDesign(n = 5L * 2L, getParamSet(fn), fun = randomLHS)
surrogate <- makeLearner("regr.ranger", predict.type = "se", keep.inbag = TRUE)
# Set general controls 在这里我定义了一些选项
ctrl <- makeMBOControl()
ctrl <- setMBOControlTermination(ctrl, iters = 10L)
ctrl <- setMBOControlInfill(ctrl, crit = makeMBOInfillCritEI())
# Start optimization 这是优化部分：
result <- mbo(fn, des, surrogate, ctrl, more.args = list("validation_data" = validation_data))
#因此推荐的参数是6为mtry和381为trees 。 RMSE的值低于之前的值，等于0.393。 现在让我们用这个值训练训练数据上的随机森林。 首先，我预处理训练数据：
training_rec <- prep(simple_recipe(housing_train), testing = housing_train)
train_data <- bake(training_rec, new_data = housing_train)
#现在让我们训练我们的最终模型并预测价格：
final_model <- rand_forest(mtry = 6, trees = 381) %>%
set_engine("ranger", importance = 'impurity') %>%
fit(price ~ ., data = train_data)
price_predict <- predict(final_model, new_data = select(test_data, -price))
#让我们将数据转换回来，并将预测的价格与真实的价格进行直观比较：
cbind(price_predict * sd(housing_train$price) + mean(housing_train$price),
housing_test$price)
#我们现在计算RMSE
tibble::tibble("truth" = test_data$price,
"prediction" = unlist(price_predict)) %>%
rmse(truth, prediction)
library("tidymodels")
library("parsnip")
library("brotools")
library("mlbench")
#随机森林预测房价
data("BostonHousing2")
head(BostonHousing2)
boston <- BostonHousing2 %>%
select(-medv, -town, -lon, -lat) %>%
rename(price = cmedv)
#删除了town ， lat和lon因为列栏中包含的信息就足够了。
#为了训练和评估模型的性能，数据分成两部分。 训练集分为两个训练集,测试集最后使用。
train_test_split <- initial_split(boston, prop = 0.9)
housing_train <- training(train_test_split)
housing_test <- testing(train_test_split)
#我想训练一个随机森林来预测房屋的价格，但随机森林有所谓的超参数，这些参数是无法从数据中估算或学习的参数。
#相反，这些参数必须由分析师选择。 为了选择它们，您可以使用似乎运行良好的文献中的值（如在宏计量经济学中完成）或者您可以进一步将列车集拆分为两个，创建超参数网格，在一个上训练模型网格的所有值的部分数据，并比较数据的第二部分上的模型的预测。
#然后，您将坚持使用最佳性能的模型，例如，具有最低RMSE的模型。 问题是，您无法仅使用一个值来估计RMSE的真实值。 这就像你想通过从人口中抽取一个单一的观察来估计人口的高度。
#为了给出一组超参数的RMSE的真实值，而不是做一次拆分，我会做30.然后计算平均RMSE，这意味着为超参数的每个值的组合训练30个模型我有兴趣。
#先让我们使用{rsample}包中的mc_cv()函数再次拆分训练数据。 此函数实现蒙特卡罗交叉验证：
validation_data <- mc_cv(housing_train, prop = 0.9, times = 30)
validation_data$splits[[1]]
#数据预处理
#我将分三步完成。 第一步和第二步用于标准化，第三步将字符和因子变量转换为虚拟变量。 这是必要的，因为随机森林不能直接处理因子变量。 让我们定义一个配方来做到这一点，并从预处理测试集开始。
#包装函数，可应用于各种数据集：
simple_recipe <- function(dataset){
recipe(price ~ ., data = dataset) %>%
step_center(all_numeric()) %>%
step_scale(all_numeric()) %>%
step_dummy(all_nominal())
}
#使用prep()函数，该函数从数据中估算处理数据所需的参数。 例如，对于居中， prep()估计平均值，然后从变量中减去平均值。
#使用bake() ，将估算值应用于数据：
testing_rec <- prep(simple_recipe(housing_test), testing = housing_test)
test_data <- bake(testing_rec, new_data = housing_test)
#在使用prep()和bake()之前分割数据很重要，因为如果没有，您将使用prep()步骤中测试集的观察结果，从而将测试集中的知识引入训练数据中。 这称为数据泄漏，必须避免。 这就是为什么有必要首先将训练数据分成分析和评估集，然后分别预处理这些集。
#但是， validation_data对象现在不能与recipe()一起使用，因为它不是数据帧。 不用担心，我只需编写一个函数，从validation_data对象中提取分析和评估集，应用预处理，训练模型，并返回RMSE。 这将是一项重要功能，是分析的核心。
#在此之前，让我们运行一个简单的线性回归作为基准。 对于线性回归，我不会使用任何CV，所以让我们预处理训练集：
trainlm_rec <- prep(simple_recipe(housing_train), testing = housing_train)
trainlm_data <- bake(trainlm_rec, new_data = housing_train)
linreg_model <- lm(price ~ ., data = trainlm_data)
broom::augment(linreg_model, newdata = test_data) %>%
rmse(price, .fitted)
#broom::augment()将预测添加到新列的.fitted 。 我不会在随机森林中使用这个技巧，因为我将使用{ranger}随机森林没有augment()方法。 我会自己将预测添加到数据中。
#现在让我们回到随机森林并编写一个大函数
my_rf <- function(mtry, trees, split, id){
analysis_set <- analysis(split)
analysis_prep <- prep(simple_recipe(analysis_set), training = analysis_set)
analysis_processed <- bake(analysis_prep, new_data = analysis_set)
model <- rand_forest(mtry = mtry, trees = trees) %>%
set_engine("ranger", importance = 'impurity') %>%
fit(price ~ ., data = analysis_processed)
assessment_set <- assessment(split)
assessment_prep <- prep(simple_recipe(assessment_set), testing = assessment_set)
assessment_processed <- bake(assessment_prep, new_data = assessment_set)
tibble::tibble("id" = id,
"truth" = assessment_processed$price,
"prediction" = unlist(predict(model, new_data = assessment_processed)))
}
#rand_forest()函数可从{parsnip}包中获得。 该软件包为许多其他机器学习包提供了统一的接口。
#这意味着不必学习range()和randomForest()的语法，而且......你可以简单地使用rand_forest()函数并将engine参数更改为你想要的randomForest（ranger， randomForest等）。
#试试这个功能
results_example <- map2_df(.x = validation_data$splits,
.y = validation_data$id,
~my_rf(mtry = 3, trees = 200, split = .x, id = .y))
head(results_example)
#我现在可以在mtry = 3和trees = 200时计算RMSE：
results_example %>%
group_by(id) %>%
rmse(truth, prediction) %>%
summarise(mean_rmse = mean(.estimate)) %>%
pull
#随机森林的RMSE已经低于线性回归。 现在的目标是通过调整mtry和trees超参数来降低此RMSE。 为此，我将使用{mlrMBO}包中实现的贝叶斯优化方法。
#贝叶斯超参数优化
#重新使用上面的代码，并定义一个函数，该函数执行从预处理到返回我希望通过调整超参数RMSE来最小化的度量的所有内容：
tuning <- function(param, validation_data){
mtry <- param[1]
trees <- param[2]
results <- purrr::map2_df(.x = validation_data$splits,
.y = validation_data$id,
~my_rf(mtry = mtry, trees = trees, split = .x, id = .y))
results %>%
group_by(id) %>%
rmse(truth, prediction) %>%
summarise(mean_rmse = mean(.estimate)) %>%
pull
}
tuning(c(3, 200), validation_data)
plot_points <- crossing("mtry" = 3, "trees" = seq(200, 300))
#绘制mtry = 3的RMSE值和200到300的trees 。因为我需要100次评估这个昂贵的函数。
# 如果评估函数，我可以通过改变mtry值来制作3D图，但是如果评估函数的价格便宜，我会进行详尽的网格搜索以找到超参数而不是使用贝叶斯优化。
plot_data <- plot_points %>%
mutate(value = map_dbl(seq(200, 300), ~tuning(c(3, .), validation_data)))
plot_data %>%
ggplot(aes(y = value, x = trees)) +
geom_line(colour = "#82518c") +
theme_blog() +
ggtitle("RMSE for mtry = 3")
#对于mtry = 3 ，最小值似乎在255左右。最小化的函数根本不是平滑的。
#我现在按照arxiv文章中的代码来运行优化。 但据我所知，一个更简单的模型，称为代理模型，用于寻找有希望的点，并在这些点评估函数的价值。 这似乎（在精神上）与Gourieroux，Monfort，Renault中描述的间接推理方法有些相似。
#让我们首先加载包并创建要优化的函数：
library("mlrMBO")
fn <- makeSingleObjectiveFunction(name = "tuning",
fn = tuning,
par.set = makeParamSet(makeIntegerParam("x1", lower = 3, upper = 8),
makeIntegerParam("x2", lower = 50, upper = 500)))
#此功能基于我之前定义的功能。 要优化的参数也被定义为它们的界限。 我会在3到8之间寻找mtry ，在50到500之间寻找trees 。
#我认为这意味着这10个点是用于启动整个过程的点。 我不明白为什么他们必须从超立方体中采样，但还可以。 然后我选择代理模型，随机森林，并预测标准误差。
library(lhs)# for randomLHS
des <- generateDesign(n = 5L * 2L, getParamSet(fn), fun = randomLHS)
surrogate <- makeLearner("regr.ranger", predict.type = "se", keep.inbag = TRUE)
# Set general controls 在这里我定义了一些选项
ctrl <- makeMBOControl()
ctrl <- setMBOControlTermination(ctrl, iters = 10L)
ctrl <- setMBOControlInfill(ctrl, crit = makeMBOInfillCritEI())
# Start optimization 这是优化部分：
result <- mbo(fn, des, surrogate, ctrl, more.args = list("validation_data" = validation_data))
#因此推荐的参数是6为mtry和381为trees 。 RMSE的值低于之前的值，等于0.393。 现在让我们用这个值训练训练数据上的随机森林。 首先，我预处理训练数据：
training_rec <- prep(simple_recipe(housing_train), testing = housing_train)
train_data <- bake(training_rec, new_data = housing_train)
#现在让我们训练我们的最终模型并预测价格：
final_model <- rand_forest(mtry = 6, trees = 381) %>%
set_engine("ranger", importance = 'impurity') %>%
fit(price ~ ., data = train_data)
price_predict <- predict(final_model, new_data = select(test_data, -price))
#让我们将数据转换回来，并将预测的价格与真实的价格进行直观比较：
cbind(price_predict * sd(housing_train$price) + mean(housing_train$price),
housing_test$price)
#我们现在计算RMSE
tibble::tibble("truth" = test_data$price,
"prediction" = unlist(price_predict)) %>%
rmse(truth, prediction)
评论
发表评论