《Tidy Modeling with R》这本书，作者是经验丰富的资深统计学家，有系统全面的理论和方法论介绍，也有作者丰富实践经验的介绍。该书从第10章《Resampling》开始的一些章节，对建模工具底下的统计学理论稍为深入介绍了一下，对已经离开学校近30年的我来说，数学推理已经看不懂了，但大致看懂会用还是没问题的。另外，这本书的源码，对于学习ggplot绘图，也是一本不错的实例教材，书中的统计图简明握要，很精当。
  总体上，该书的主题不是介绍各种模型的原理，那已经有很多书了，而是介绍用Posit的tidymodels工具套件完成分析建模及评估优化的完整流程。tidymodels套件的设计理念是易用、复用、管道和函数式编程，这让R语言的分析建模代码很简洁，也很容易阅读。确实是易于使用，所以在执行性能上就要稍微作出点牺牲。软件永远都是这样，人们需要作出平衡和选择，屏蔽了复杂性，易于使用就意味着后台干了更多的工作，就要有消耗。
  具体来说，recipes包等提供了丰富的预处理工具集，一个函数调用就可以完成一种复杂的预处理，所以它们是组特征工程的工具。parsnip包尝试统一各种算法模型的建模调用接口，让参数和调用更一致更易用易记。workflows及workflowsets包把特征工程和建模连接成完整的工作流。rsample包提供了数据集划分、重抽样、并行多重交叉验证等的实现。yarstick包提供了评价模型性能的指标集实现。tune及dials包等提供了贝叶斯优化等模型调参的支持。purrr包则提供了函数式编程的支持（类似于apply族函数的矢量运算等）。其它的包如tibble、tidyr包等是格式辅助包，dplyr是高效操作大数据集的非SQL语法工具包，比如稍后的filter()过滤数据及arrange()排序。最后，这些部件都通过管道操作连接，数据在管道内流动，这使得R语言程序非常简洁和易于阅读。不过tidymodels套件不包含模型解释的工具，它使用DALEXtra包支持模型的SHAP值解释。

> library(tidymodels)
── Attaching packages ─────────────────────── tidymodels 1.0.0 ──
✔ broom        1.0.1      ✔ recipes      1.0.3 
✔ dials        1.1.0      ✔ rsample      1.1.1 
✔ dplyr        1.0.10     ✔ tibble       3.1.8 
✔ ggplot2      3.4.0      ✔ tidyr        1.2.1 
✔ infer        1.0.4      ✔ tune         1.0.1 
✔ modeldata    1.0.1      ✔ workflows    1.1.2 
✔ parsnip      1.0.3      ✔ workflowsets 1.0.0 
✔ purrr        1.0.0      ✔ yardstick    1.1.0 
── Conflicts ────────────────────────── tidymodels_conflicts() ──
✖ purrr::discard() masks scales::discard()
✖ dplyr::filter()  masks stats::filter()
✖ dplyr::lag()     masks stats::lag()
✖ recipes::step()  masks stats::step()
• Use tidymodels_prefer() to resolve common conflicts.

  本篇先用随机森林算法演示一下。

一、加载经预处理的数据
  使用的数据与Python版相同，已经过预处理，稍后再研究一下用R语言做相同的数据探索与预处理过程。

library(tidymodels)
library(kableExtra)
library(tidyr)
library(ggplot2)
library(showtext)
# R Plot添加中文及其他字体[showtext]
# https://blog.csdn.net/weixin_46128755/article/details/125825935
showtext_auto()

# # Windows，注册并行处理
# library(doParallel)
# cl <- makePSOCKcluster(parallel::detectCores())
# registerDoParallel(cl)
# Unix like only
library(doMC)
registerDoMC(cores = parallel::detectCores())

# 优先使用tidymodels的同名函数。
tidymodels_prefer()

# 异常值阈值30，作图用。
threshold<- 30

# ----------------------------------------------------------------------------------------
# 加载经过预处理的数据
melbourne<- read.csv("/home/jean/data/Melbourne_housing_pre.csv")
# 过滤缺失值, 已经作过填零处理，会导致下面的错误。
# Error: Missing data in columns: BuildingArea.
# 47 obs.
missing <- filter(melbourne, BuildingArea==0)
melbourne <- filter(melbourne, BuildingArea!=0)

# 划分训练集与测试集，这些都是rsample包的函数。
set.seed(2023)
melbourne_split <- initial_split(melbourne, prop = 0.80)
melbourne_train <- training(melbourne_split)
melbourne_test  <-  testing(melbourne_split)

> melbourne_split
<Training/Testing/Total>
<7174/1794/8968>

二、随机森林模型
  Tidy Models数据预处理由recipes等包提供，由工作流连接特征工程预处理与模型，预处理与建模等主要通过管道操作符"%>%"连接操作，R-4.1后也可用R语言内生的管道操作符"|>"。
  R语言公式化建模是与Python的重要不同，它的公式表达式可以高效表达比较复杂的建模公式，比如广义线性回归时，y~I(x)^n就可以对x的n次多项式建模，还可以在公式中简洁的指定输入变量之间的交叉项等。这是它方便的地方。

# 定义预处理菜谱：回归公式与预处理，参与回归的变量与Python版相同，没有交叉项。
# LogPrice ~'Year','YearBuilt','Distance','Lattitude','Longtitude','Propertycount',
# 'Landsize','BuildingArea', 'Rooms','Bathroom', 'Car','Type_h','Type_t','Type_u'
melbourne_rec<-    # 定义回归公式，目标变量是LogPrice 
  recipe(LogPrice ~ Year + YearBuilt + Distance + Lattitude + Longtitude + BuildingArea
         + Rooms + Bathroom + Car + Type_h + Type_t + Type_u, data = melbourne_train) %>%
  # 标准化数值型变量，这些是recipes等包提供的预处理步骤函数。
  # 随机森林不需要标准化数据量纲，SVM等需要，这里是演示。
  #step_log(BuildingArea, base = 10) %>%
  step_normalize(all_numeric_predictors())

# 定义模型：随机森林  
rf_model <- 
  rand_forest(trees = 1000) %>%  # parsnip模型主参数
  set_engine("ranger") %>%       # 使用ranger实现的随机森林算法
  set_mode("regression")         # 用于回归分析

# 定义工作流，组合预处理与模型
rf_wflow <- 
  workflow() %>% 
  add_model(rf_model) %>% 
  add_recipe(melbourne_rec)

# 训练模型
t1<-proc.time()
rf_fit <- fit(rf_wflow, melbourne_train)
t2<-proc.time()
# 8.39 0.16 8.55 NA NA
cat(t2-t1)

# 测试集预测
# 预测值
melbourne_test_res <- predict(rf_fit, new_data = melbourne_test %>% select(-LogPrice))
# 合并真实值
melbourne_test_res <- bind_cols(melbourne_test_res, melbourne_test %>% select(LogPrice))
# 浏览
melbourne_test_res

# 画拟合效果，这个图效果较好
# 可以看到随机森林预测的效果是：异常值中低价区间预测偏高，高价区间预测偏低。
data <- arrange(melbourne_test_res, LogPrice)
data["index"]<-1:nrow(data)
data["upper"]<-data["LogPrice"]+log(1+threshold/100)
data["lower"]<-data["LogPrice"]+log(1-threshold/100)

ggplot(data = data)+
  geom_point(mapping = aes(x = index, y = .pred, color = "blue"), size=0.1)+
  geom_line(mapping = aes(x = index, y = LogPrice, color = "red"))+
  geom_ribbon(aes(index, ymax=upper, ymin=lower, fill="yellow"), alpha=0.25)+
  ggtitle("测试集拟合效果图")+
  xlab("样本")+
  ylab("ln(房价)")+
  scale_fill_identity(name = '异常值区间', guide = 'legend',labels = c('正常值')) +
  scale_colour_manual(name = '颜色', 
                      values =c('blue'='blue','red'='red'), labels = c('预测值','实际值'))

# 计算性能指标
# rmse(melbourne_test_res, truth = LogPrice, estimate = .pred)
# 0.852, 0.203, 0.144
melbourne_metrics <- metric_set(rsq, rmse,  mae)
melbourne_metrics(melbourne_test_res, truth = LogPrice, estimate = .pred)

测试集上准确率rsq=85.2%。rsq即统计量，即模型反映了数据中方差的百分比，在回归分析中相当于准确率。因为回归分析一般使用rmse，标准差，越小越好。

> melbourne_metrics(melbourne_test_res, truth = LogPrice, estimate = .pred)
# A tibble: 3 × 3
  .metric .estimator .estimate
  <chr>   <chr>          <dbl>
1 rsq     standard       0.852
2 rmse    standard       0.203
3 mae     standard       0.144

测试集拟合效果是异常值中低价位预测偏高及高价位预测偏低的相对多一点。

测试集拟合效果

三、多折交叉验证
  多折交叉验证使用多核CPU多进程并行处理，5折是Python的默认。

set.seed(1001)
# 定义5折交叉验证
melbourne_folds <- vfold_cv(melbourne_train, v = 5)
# 用以在交叉验证抽取训练的模型
get_model <- function(x) {
  extract_fit_parsnip(x)
}

# 并行训练，用fit_resamples()函数。
t1<-proc.time()
# 保留交叉验证预测结果，保留训练的工作流对象，抽取训练时拟合的模型
keep_pred <- control_resamples(save_pred = TRUE, save_workflow = TRUE, extract = get_model)
set.seed(1003)
rf_res <- 
  rf_wflow %>% 
  fit_resamples(resamples = melbourne_folds, control = keep_pred)
t2<-proc.time()
# 1.11 0.32 11.98 NA NA
cat(t2-t1)
# 看一下结果
rf_res

# 查看性能指标
# 0.875, 0.191
collect_metrics(rf_res)

# 查看预测效果
assess_res <- collect_predictions(rf_res)
assess_res

# 画预测效果图
data <- arrange(assess_res, LogPrice)
data["index"]<-1:nrow(data)
data["upper"]<-data["LogPrice"]+log(1+threshold/100)
data["lower"]<-data["LogPrice"]+log(1-threshold/100)

ggplot(data = data)+
  geom_point(mapping = aes(x = index, y = .pred, color = "blue"), size=0.1)+
  geom_line(mapping = aes(x = index, y = LogPrice, color = "red"))+
  geom_ribbon(aes(index, ymax=upper, ymin=lower, fill="yellow"), alpha=0.25)+
  ggtitle("五折交叉验证拟合效果图")+
  xlab("样本")+
  ylab("ln(房价)")+
  scale_fill_identity(name = '异常值区间', guide = 'legend',labels = c('正常值')) +
  scale_colour_manual(name = '颜色', 
                      values =c('blue'='blue','red'='red'), labels = c('预测值','实际值'))

5折交叉验证把训练集划分成5等份，每次用占20%的一份作验证集，其它用于训练。验证集上的准确率是87.5%。这样就可以及早处理训练集中的过拟合问题。

> rf_res
# Resampling results
# 5-fold cross-validation 
# A tibble: 5 × 6
  splits              id    .metrics         .notes           .extracts        .predictions        
  <list>              <chr> <list>           <list>           <list>           <list>              
1 <split [5739/1435]> Fold1 <tibble [2 × 4]> <tibble [0 × 3]> <tibble [1 × 2]> <tibble [1,435 × 4]>
2 <split [5739/1435]> Fold2 <tibble [2 × 4]> <tibble [0 × 3]> <tibble [1 × 2]> <tibble [1,435 × 4]>
3 <split [5739/1435]> Fold3 <tibble [2 × 4]> <tibble [0 × 3]> <tibble [1 × 2]> <tibble [1,435 × 4]>
4 <split [5739/1435]> Fold4 <tibble [2 × 4]> <tibble [0 × 3]> <tibble [1 × 2]> <tibble [1,435 × 4]>
5 <split [5740/1434]> Fold5 <tibble [2 × 4]> <tibble [0 × 3]> <tibble [1 × 2]> <tibble [1,434 × 4]>
> collect_metrics(rf_res)
# A tibble: 2 × 6
  .metric .estimator  mean     n std_err .config             
  <chr>   <chr>      <dbl> <int>   <dbl> <chr>               
1 rmse    standard   0.191     5 0.00312 Preprocessor1_Model1
2 rsq     standard   0.875     5 0.00402 Preprocessor1_Model1

五折交叉验证拟合效果图

四、网格搜索
  这是比较简单的模型超参数优化方式，在参数空间中作有限次数的网格式搜索，稍为了解一下。

# 定义菜谱：回归公式与预处理，这里重新写一下是为了阅读清晰。
melbourne_rec<-
  recipe(LogPrice ~ Year + YearBuilt + Distance + Lattitude + Longtitude + BuildingArea
         + Rooms + Bathroom + Car + Type_h + Type_t + Type_u, data = melbourne_train) %>%
  # 标准化数值型变量
  step_normalize(all_numeric_predictors())

# 定义模型：随机森林， 定义要调整的参数，这是与上面不同的地方。
# tune()函数起占位符的作用，实际值由后面tune_grid()中的参数网格赋值。
# parsnip模型中只有这3个主参数可调。
rf_spec <- 
  rand_forest(mtry = tune(), trees = tune(), min_n = tune()) %>% 
  set_engine("ranger") %>% 
  set_mode("regression")

# 定义工作流
rf_wflow <- 
  workflow() %>% 
  add_model(rf_spec) %>% 
  add_recipe(melbourne_rec)

# mtry参数的边界未完全确定，用finalize()函数确定。
rf_param <- rf_wflow %>%
  extract_parameter_set_dials() %>%
  finalize(melbourne_train) # 要带上训练集作参数。
# 查看参数边界，都已确定
rf_param %>% extract_parameter_dials("trees")
rf_param %>% extract_parameter_dials("mtry")
rf_param %>% extract_parameter_dials("min_n")

# 执行网格搜索，每个参数搜索3个候选，每组参数用5折交叉验证。
set.seed(2023)
t1<-proc.time()
rf_res_grid <- 
  rf_wflow %>% 
  tune_grid(
    resamples = melbourne_folds,
    grid = rf_param %>% grid_regular(levels = 3),  # 用参数空间构造一个网格
    metrics = metric_set(rsq, rmse, mae)
  )
t2<-proc.time()
# 0.36 0 224.35 NA NA
cat(t2-t1)

# 画图查看网格搜索效果
autoplot(rf_res_grid)

# 查看准确率最高的模型
show_best(rf_res_grid, metric="rsq")

# 选择准确率最高的模型
select_best(rf_res_grid, metric="rsq")

# 构造最佳参数
# 直接读取调参的最佳结果
rf_param_best<- select_best(rf_res_grid, metric="rsq")

# 最佳参数回填到工作流
rf_wflow_grid <-
  rf_wflow %>%
  finalize_workflow(rf_param_best)
rf_wflow_grid

# 用最佳参数在训练集全集上训练模型
t1<-proc.time()
rf_fit_grid<- rf_wflow_grid %>% fit(melbourne_train)
t2<-proc.time()
# 28.64 0.39 29.02 NA NA
cat(t2-t1)

# 预测值
melbourne_test_grid <- predict(rf_fit_grid, new_data = melbourne_test %>% select(-LogPrice))
# 合并真实值
melbourne_test_grid <- bind_cols(melbourne_test_grid, melbourne_test %>% select(LogPrice))
# 浏览
melbourne_test_grid
# 85.3, 0.201, 0.142
melbourne_metrics(melbourne_test_grid, truth = LogPrice, estimate = .pred)

finalize()之前的mtry参数，带？表示边界未定，因为它一般是变量数量的平方根，要根据具体的训练集确定。

> rf_param <- rf_wflow %>%
+     extract_parameter_set_dials()
> rf_param 
Collection of 3 parameters for tuning

 identifier  type    object
       mtry  mtry nparam[?]
      trees trees nparam[+]
      min_n min_n nparam[+]

Model parameters needing finalization:
   # Randomly Selected Predictors ('mtry')

See `?dials::finalize` or `?dials::update.parameters` for more information.

finalize()后mtry 的边界就确定了，然后就可以运行 tune_grid()调参。

> rf_param <- rf_wflow %>%
+     extract_parameter_set_dials() %>%
+     finalize(melbourne_train) # 要带上训练集作参数。
> 
> rf_param
Collection of 3 parameters for tuning

 identifier  type    object
       mtry  mtry nparam[+]
      trees trees nparam[+]
      min_n min_n nparam[+]

网格搜索各参数组合效果

显示与选择最佳模型。

> show_best(rf_res_grid, metric="rsq")
# A tibble: 5 × 9
   mtry trees min_n .metric .estimator  mean     n std_err .config              
  <int> <int> <int> <chr>   <chr>      <dbl> <int>   <dbl> <chr>                
1     8  2000     2 rsq     standard   0.874     3 0.00564 Preprocessor1_Model08
2     8  1000     2 rsq     standard   0.874     3 0.00549 Preprocessor1_Model05
3     8  1000    21 rsq     standard   0.871     3 0.00516 Preprocessor1_Model14
4     8  2000    21 rsq     standard   0.871     3 0.00522 Preprocessor1_Model17
5    16  2000     2 rsq     standard   0.868     3 0.00517 Preprocessor1_Model09
> select_best(rf_res_grid, metric="rsq")
# A tibble: 1 × 4
   mtry trees min_n .config              
  <int> <int> <int> <chr>                
1     8  2000     2 Preprocessor1_Model08

最佳模型在测试集上的准确率是85.3%，网格搜索有所提高（上面是85.2%）。

> melbourne_metrics(melbourne_test_grid, truth = LogPrice, estimate = .pred)
# A tibble: 3 × 3
  .metric .estimator .estimate
  <chr>   <chr>          <dbl>
1 rsq     standard       0.853
2 rmse    standard       0.201
3 mae     standard       0.142

五、贝叶斯优化
  贝叶斯优化可以调整菜谱参数、模型主参数及引擎相关参数，具体介绍请参阅《Tuning Parameters in tidymodels》一节，随机森林ranger实现可调的参数可以参阅turnRanger包文档，这里增加了引擎相关的参数regularization.factor及sample.fraction，本篇没有要调整的recipes预处理参数。贝叶斯优化函数tune_bayes()的参数no_improve表示经过该步数的搜索后，如果找不到更好的参数组合就停止搜索，这里设为Inf，表示不停止，执行完参数iter中定义的搜索次数。
  因为随机森林算法的实现不同，引擎相关的参数名字及作用也不同，与Python版不是一 一对应（Python版是调整7个超参数，此处是5个）。Tidy Models中的贝叶斯优化体系等用起来很方便，此处演示了在dials包框架内调整引擎相关参数的方法，但具体调整以达到Python版相同的性能，需要继续深入研究一下。
  dials包可调的参数有限，可以在该页面查询。定义新的dials包可调参数请参阅文章《How to create a tuning parameter function》。

# 定义菜谱：回归公式与预处理
melbourne_rec<-
  recipe(LogPrice ~ Year + YearBuilt + Distance + Lattitude + Longtitude + BuildingArea
         + Rooms + Bathroom + Car + Type_h + Type_t + Type_u, data = melbourne_train) %>%
  # 标准化数值型变量
  step_normalize(all_numeric_predictors())

# 定义模型：随机森林， 定义要调整的参数 
rf_spec <- 
  rand_forest(mtry = tune(), trees = tune(), min_n = tune()) %>%  # 模型主参数
  set_engine("ranger", 
             regularization.factor = tune("regularization"),      # 引擎相关参数
             sample.fraction = tune("max_samples")) %>%           # 引擎相关参数
  set_mode("regression")

# 定义工作流
rf_wflow <- 
  workflow() %>% 
  add_model(rf_spec) %>% 
  add_recipe(melbourne_rec)

# 定义一个引擎相关的新dials参数函数
max_samples <- function(range = c(0, 1), trans = NULL) {
  new_quant_param(
    type = "double",
    range = range,
    inclusive = c(FALSE, TRUE),
    trans = trans,
    label = c(sample.fraction = "# Fraction of observations to sample"),
    finalize = NULL
  )
}

# mtry参数的边界未完全确定，用finalize()函数确定。
rf_param <- rf_wflow %>%
  extract_parameter_set_dials() %>%
  # 注意用法：上面sample.fraction = tune("max_samples"))定义了参数的ID
  # 引用新定义的dials参数函数max_samples给出参数的取值范围。
  update(max_samples = max_samples()) %>%              # 暂时这样用着
  finalize(melbourne_train)                            # 要带上训练集作参数。
# 查看参数边界，都已确定
rf_param %>% extract_parameter_dials("trees")
rf_param %>% extract_parameter_dials("mtry")
rf_param %>% extract_parameter_dials("min_n")
rf_param %>% extract_parameter_dials("regularization")
rf_param %>% extract_parameter_dials("max_samples")

# 执行贝叶斯优化
ctrl <- control_bayes(verbose = TRUE, no_improve = Inf)
set.seed(2023)
t1<-proc.time()
rf_res_bo <-
  rf_wflow %>%
  tune_bayes(
    resamples = melbourne_folds,
    metrics = metric_set(rsq, rmse, mae),
    initial = 10,
    param_info = rf_param,
    iter = 100,
    control = ctrl
  )
t2<-proc.time()
cat(t2-t1)

# 画图查看贝叶斯优化效果
autoplot(rf_res_bo, type = "performance")
autoplot(rf_res_bo, type = "performance", metric = "rsq")

# 查看准确率最高的模型
show_best(rf_res_bo, metric="rsq")

# 选择准确率最高的模型
select_best(rf_res_bo, metric="rsq")

# 直接读取调参的最佳结果
rf_param_best<- select_best(rf_res_bo, metric="rsq")

# 最佳参数回填到工作流
rf_wflow_bo <-
  rf_wflow %>%
  finalize_workflow(rf_param_best)
rf_wflow_bo

# 用最佳参数在训练集全集上训练模型
t1<-proc.time()
rf_fit_bo<- rf_wflow_bo %>% fit(melbourne_train)
t2<-proc.time()
# 28.64 0.39 29.02 NA NA
cat(t2-t1)

# 测试集
# 预测值
melbourne_test_bo <- predict(rf_fit_bo, new_data = melbourne_test %>% select(-LogPrice))
# 合并真实值
melbourne_test_bo <- bind_cols(melbourne_test_bo, melbourne_test %>% select(LogPrice))
# 浏览
melbourne_test_bo
# 0.853, 0.201, 0.142
melbourne_metrics <- metric_set(rsq, rmse,  mae)
melbourne_metrics(melbourne_test_bo, truth = LogPrice, estimate = .pred)

# 训练集
# 预测值
melbourne_train_bo <- predict(rf_fit_bo, new_data = melbourne_train %>% select(-LogPrice))
# 合并真实值
melbourne_train_bo <- bind_cols(melbourne_train_bo, melbourne_train %>% select(LogPrice))
# 浏览
melbourne_train_bo
# 0.984, 0.0712, 0.0524
melbourne_metrics(melbourne_train_bo, truth = LogPrice, estimate = .pred)

开始时参数mtry是?状态，参数sample.fraction是missing，都需要指定边界。

> rf_param <- rf_wflow %>%
+   extract_parameter_set_dials() 
> rf_param
Collection of 5 parameters for tuning

     identifier                  type    object
           mtry                  mtry nparam[?]
          trees                 trees nparam[+]
          min_n                 min_n nparam[+]
 regularization regularization.factor nparam[+]
    max_samples       sample.fraction   missing

The parameter `max_samples` needs a `param` object. 
See `vignette('dials')` to learn more.
Model parameters needing finalization:
   # Randomly Selected Predictors ('mtry')

See `?dials::finalize` or `?dials::update.parameters` for more information.

现在它们的边界都确定了，可以运行tune_bayes()了。

> rf_param <- rf_wflow %>%
+   extract_parameter_set_dials() %>%
+   update(max_samples = max_samples()) %>%              # 暂时这样用着
+   finalize(melbourne_train)                                # 要带上训练集作参数。
> rf_param
Collection of 5 parameters for tuning

     identifier                  type    object
           mtry                  mtry nparam[+]
          trees                 trees nparam[+]
          min_n                 min_n nparam[+]
 regularization regularization.factor nparam[+]
    max_samples       sample.fraction nparam[+]

> rf_param %>% extract_parameter_dials("mtry")
# Randomly Selected Predictors (quantitative)
Range: [1, 16]
> rf_param %>% extract_parameter_dials("max_samples")
# Fraction of observations to sample (quantitative)
Range: (0, 1]

迭代优化的过程。

── Iteration 56 ──────────────────────────────────────────

i Current best:     rsq=0.8819 (@iter 17)
i Gaussian process model
✓ Gaussian process model
i Generating 5000 candidates
i Predicted candidates
i mtry=5, trees=1254, min_n=10, regularization=0.0577, max_samples=0.671
i Estimating performance
✓ Estimating performance
ⓧ Newest results:   rsq=0.8754 (+/-0.0038)

── Iteration 57 ──────────────────────────────────────────

i Current best:     rsq=0.8819 (@iter 17)
i Gaussian process model
✓ Gaussian process model
i Generating 5000 candidates
i Predicted candidates
i mtry=4, trees=1905, min_n=3, regularization=0.472, max_samples=0.776
i Estimating performance
✓ Estimating performance
♥ Newest results:   rsq=0.8823 (+/-0.00432)

── Iteration 58 ──────────────────────────────────────────

i Current best:     rsq=0.8823 (@iter 57)
i Gaussian process model
✓ Gaussian process model
i Generating 5000 candidates
i Predicted candidates
i mtry=3, trees=640, min_n=6, regularization=0.447, max_samples=0.98
i Estimating performance
✓ Estimating performance
ⓧ Newest results:   rsq=0.8755 (+/-0.00405)

> # 查看准确率最高的模型
> show_best(rf_res_bo, metric="rsq")
# A tibble: 5 × 12
   mtry trees min_n regularization max_samples .metric .estimator  mean     n std_err .config .iter
  <int> <int> <int>          <dbl>       <dbl> <chr>   <chr>      <dbl> <int>   <dbl> <chr>   <int>
1     4  1905     3          0.472       0.776 rsq     standard   0.882     3 0.00432 Iter57     57
2     7  1918     2          0.955       0.992 rsq     standard   0.882     3 0.00391 Iter17     17
3     7  1737     2          0.494       0.824 rsq     standard   0.882     3 0.00406 Iter11     11
4     8  1906     3          0.453       0.981 rsq     standard   0.881     3 0.00403 Iter94     94
5     8  1843     3          0.639       0.845 rsq     standard   0.880     3 0.00395 Iter80     80
> # 选择准确率最高的模型
> select_best(rf_res_bo, metric="rsq")
# A tibble: 1 × 6
   mtry trees min_n regularization max_samples .config
  <int> <int> <int>          <dbl>       <dbl> <chr>  
1     4  1905     3          0.472       0.776 Iter57 
> 

参数多于3个时，不画参数的组合效果图，因为组合太多了，可视化效果不好，这里画性能指标的效果图看看。

迭代过程中性能的变化

另一次搜索迭代中rsq的优化过程

测试集上的准确率提高到85.4%，又改进了0.1%。

> melbourne_metrics(melbourne_test_bo, truth = LogPrice, estimate = .pred)
# A tibble: 3 × 3
  .metric .estimator .estimate
  <chr>   <chr>          <dbl>
1 rsq     standard       0.854
2 rmse    standard       0.201
3 mae     standard       0.142

训练集上的准确率是97.1%，相差较大，训练集上过拟合比较严重。

> melbourne_metrics(melbourne_train_bo, truth = LogPrice, estimate = .pred)
# A tibble: 3 × 3
  .metric .estimator .estimate
  <chr>   <chr>          <dbl>
1 rsq     standard      0.971 
2 rmse    standard      0.0954
3 mae     standard      0.0702

检查了一下原来Python随机森林调参的代码，当时是用了极致的方法，没有多折交叉验证，训练集全集训练然后直接取测试集上的最高准确率，所以准确率虽然高（88.3%），但是不稳定的。此处调参则用了5折交叉验证，要稳定一点，也慢一点，然后重新在训练集全集上训练，训练数据与5折交叉验证又不同了（多了20%）。结果就是上面5折交叉验证有超过88%的准确率，训练集全集重新训练后，测试集的准确率就降到85.4%了，但这样稳定性和泛化性能应该是更好的。Python版也可以用同样的方法处理一次再比较一下，要方法相同才可比较。
不过，本篇到目前为止对ranger随机森林调参的效果有限，暂时还找不到本例中随机森林算法抗过拟合的有效方法。
Python版调参程序：

from sklearn.model_selection import train_test_split, KFold, cross_val_score
from sklearn.ensemble import RandomForestRegressor
from hyperopt import fmin, tpe, hp, Trials
from hyperopt.fmin import generate_trials_to_calculate

# Auto search for better hyper parameters with hyperopt, only need to give a range
'''
1: max_depth，树深，3~50性能较好，过了10后测试集性能迅速下降，过拟合明显。有些问题可能需要较大的深度，如本例。
2: min_sample_split，分裂阈值，树的结点要分裂生长需要具有超过阈值的样本数，增大可以克服过拟合，过大会欠拟合。
3：max_leaf_nodes，树的最大叶子数，超过该值，树不再分裂生长。小值欠拟合，过大过拟合。本例要设置在1000~2000之间。
4: min_samples_leaf，叶子结点最少样本数，少于该阈值则树不再分裂生长。增大可以防止过拟合。
   <100，容易过拟合，100~400性能较好，>500欠拟合明显。本例需要设置较小的值，如=1。
5: n_estimators，树数，越大越好，但计算开销增大， >300无明显改进。
6: max_features，树采用的特征数，3~10，经过一定的值后开始过拟合，默认特征数的平方根或6最佳。
7: max_samples，数据采样率，(0,1],越低速度越快，但精度越低。超过一定的采样率，精度改善不明显。0.2以上合适。
'''
# 缩小参数取值范围，搜索会快很多  
space_rf = {
    'max_depth': hp.choice('max_depth', range(3, 51)), 
    'min_samples_split': hp.choice('min_samples_split', range(2, 11)),  
    'max_leaf_nodes': hp.choice('max_leaf_nodes', range(2, 2001)),      
    'min_samples_leaf': hp.choice('min_samples_leaf', range(1, 501)), 
    'n_estimators': hp.choice('n_estimators', range(50, 501)),    
    'max_features': hp.choice('max_features', range(3, 10)),
    'max_samples': hp.uniform('max_samples', 0.2, 1.0),                         
}

def f_rf(params):
    rf = RandomForestRegressor( random_state = 0,verbose=0, n_jobs = -1, **params)  
    rf_model = rf.fit(train_X, train_y)
    acc = rf_model.score(valid_X,valid_y)
    #acc = cross_val_score(rf, train_X, train_y).mean()              
    return -acc
# trials = Trials()
# Set initial values, start searching from the best point of GridSearchCV(), and default values
trials = generate_trials_to_calculate([{
                                        'max_depth':29-3,                  # default None
                                        'min_samples_split':2-2,            # default 2 , >1                                       
                                        'max_leaf_nodes':1783-2,            # default None
                                        'min_samples_leaf':1-1,             # default 1
                                        'n_estimators':431-50,              # default 100                                        
                                        'max_features':7-3,                 # default auto, sqrt(n)
                                        'max_samples': 0.9703002612349153,  # default None                                        
                                        }])

t1 = time.time()  
# 1000trial [28:56,  1.74s/trial, best loss: -0.8875152912944326]
best_params = fmin(f_rf, space_rf, algo=tpe.suggest, max_evals=999, trials=trials)
t2 = time.time()
print("Time elapsed: ", t2-t1)

print('best:')
print(best_params)

六、Tidy Models支持的算法模型
Rstudio->Tools->Addins->parsnip->Execute，可以列出Tidy Models现在支持的算法模型，并且可以辅助生成使用模型的代码模板到代码编辑器当前位置。比如下面是支持的回归模型，墨尔本房价回归模型实例使用的SVM、随机森林与GBDT(xgboost)算法都有，后面将继续完成R语言上的这个实例：

Tidy Models支持的算法模型

Tidy Models支持的算法模型（续）

经过测试，XGboost的boost_tree算法在本数据集上的性能要比随机森林好，经过100次迭代的贝叶斯优化，找到的参数组合，训练集上的准确率是95.1%，5折交叉验证的准确率是88.82%，最终测试集上的准确率是87%，抗过拟合能力不错，并且5折交叉验证训练的速度很快。说明R语言的实现中，至少GBDT类算法是可以在生产中实际使用的。
XGB算法贝叶斯优化。

# 定义菜谱：回归公式与预处理
melbourne_rec<-
  recipe(LogPrice ~ Year + YearBuilt + Distance + Lattitude + Longtitude + BuildingArea
         + Rooms + Bathroom + Car + Type_h + Type_t + Type_u, data = melbourne_train) %>%
  # 标准化数值型变量
  step_normalize(all_numeric_predictors())

# 定义模型：XGB， 定义要调整的参数 
xgb_spec <-
  boost_tree(tree_depth = tune(), trees = tune(), learn_rate = tune(), min_n = tune(), loss_reduction = tune(), sample_size = tune(), stop_iter = tune()) %>%
  set_engine('xgboost') %>%
  set_mode('regression')

# 定义工作流
xgb_wflow <- 
  workflow() %>% 
  add_model(xgb_spec) %>% 
  add_recipe(melbourne_rec)

# 全部参数的边界都已确定。
xgb_param <- xgb_wflow %>%
  extract_parameter_set_dials()
xgb_param
# 查看参数边界，都已确定
xgb_param %>% extract_parameter_dials("trees")
xgb_param %>% extract_parameter_dials("min_n")
xgb_param %>% extract_parameter_dials("tree_depth")
xgb_param %>% extract_parameter_dials("learn_rate")
xgb_param %>% extract_parameter_dials("loss_reduction")
xgb_param %>% extract_parameter_dials("sample_size")
xgb_param %>% extract_parameter_dials("stop_iter")

melbourne_folds <- vfold_cv(melbourne, v = 5)

# 执行贝叶斯优化
ctrl <- control_bayes(verbose = TRUE, no_improve = Inf)
# set.seed(2023)
t1<-proc.time()
xgb_res_bo <-
  xgb_wflow %>%
  tune_bayes(
    resamples = melbourne_folds,
    metrics = metric_set(rsq, rmse, mae),
    initial = 10,
    param_info = xgb_param,
    iter = 100,
    control = ctrl
  )
t2<-proc.time()
# 3014 1.99 3989.22 NA NA
cat(t2-t1)

# 画图查看贝叶斯优化效果
autoplot(xgb_res_bo, type = "performance", metric="rsq")

# 查看准确率最高的模型
show_best(xgb_res_bo, metric="rsq")

# 选择准确率最高的模型
select_best(xgb_res_bo, metric="rsq")
# # A tibble: 1 × 8   88.82%
# trees min_n tree_depth learn_rate loss_reduction sample_size stop_iter .config
# <int> <int>      <int>      <dbl>          <dbl>       <dbl>     <int> <chr>  
# 1792    13          8     0.0138         0.0115       0.334        11 Iter84 
# 直接读取调参的最佳结果
xgb_param_best<- select_best(xgb_res_bo, metric="rsq")

# 最佳参数回填到工作流
xgb_wflow_bo <-
  xgb_wflow %>%
  finalize_workflow(xgb_param_best)
xgb_wflow_bo

# 用最佳参数在训练集全集上训练模型
t1<-proc.time()
xgb_fit_bo<- xgb_wflow_bo %>% fit(melbourne_train)
t2<-proc.time()
# 28.64 0.39 29.02 NA NA
cat(t2-t1)

# 测试集
# 预测值
melbourne_test_bo <- predict(xgb_fit_bo, new_data = melbourne_test %>% select(-LogPrice))
# 合并真实值
melbourne_test_bo <- bind_cols(melbourne_test_bo, melbourne_test %>% select(LogPrice))
# 浏览
melbourne_test_bo
melbourne_metrics(melbourne_test_bo, truth = LogPrice, estimate = .pred)

# 训练集
# 预测值
melbourne_train_bo <- predict(xgb_fit_bo, new_data = melbourne_train %>% select(-LogPrice))
# 合并真实值
melbourne_train_bo <- bind_cols(melbourne_train_bo, melbourne_train %>% select(LogPrice))
# 浏览
melbourne_train_bo
melbourne_metrics(melbourne_train_bo, truth = LogPrice, estimate = .pred)

> show_best(xgb_res_bo, metric="rsq")
# A tibble: 5 × 14
  trees min_n tree_depth learn_rate loss_reduction sampl…¹ stop_…² .metric .esti…³  mean     n std_err .config .iter
  <int> <int>      <int>      <dbl>          <dbl>   <dbl>   <int> <chr>   <chr>   <dbl> <int>   <dbl> <chr>   <int>
1  1792    13          8    0.0138    0.0115         0.334      11 rsq     standa… 0.888     5 0.00445 Iter84     84
2  1660    15          9    0.00874   0.000000457    0.727       3 rsq     standa… 0.888     5 0.00459 Iter99     99
3  1728    36          6    0.0268    0.0000000931   0.690       5 rsq     standa… 0.888     5 0.00440 Iter69     69
4  1545    23         11    0.0155    0.0380         0.496      17 rsq     standa… 0.888     5 0.00441 Iter54     54
5  1144     6          8    0.0247    0.000000250    0.743       8 rsq     standa… 0.887     5 0.00419 Iter72     72
# … with abbreviated variable names ¹sample_size, ²stop_iter, ³.estimator
> 

XGB算法在测试集上的准确率。

> melbourne_metrics(melbourne_test_bo, truth = LogPrice, estimate = .pred)
# A tibble: 3 × 3
  .metric .estimator .estimate
  <chr>   <chr>          <dbl>
1 rsq     standard       0.870
2 rmse    standard       0.190
3 mae     standard       0.136

XGB算法在训练集上的准确率。

> melbourne_metrics(melbourne_train_bo, truth = LogPrice, estimate = .pred)
# A tibble: 3 × 3
  .metric .estimator .estimate
  <chr>   <chr>          <dbl>
1 rsq     standard      0.951 
2 rmse    standard      0.119 
3 mae     standard      0.0882

本篇的篇幅已经比较大了，先到此为止，其它内容另开新篇。