機器學習及演算法

林嶔 (Lin, Chin)

Lesson 16 機器學習概論7(梯度提升機與集成模型)

第一節:梯度提升機(1)

– 當時我們有強調還有一種東西叫做「boosting」,你還記得是什麼嘛?

– 本週讓我們直接使用真實資料來進行實驗吧!請至這裡下載範例資料

library(mice)

dat <- read.csv("ECG_train.csv", header = TRUE, fileEncoding = 'CP950', stringsAsFactors = FALSE, na.strings = "")
subdat <- dat[!(dat[,"K"] %in% NA), c(-1, -3, -4, -5)]

subdat[,'GENDER'] <- as.factor(subdat[,'GENDER'])
for (i in 1:31) {subdat[,paste0('rhythm.', i)] <- as.factor(subdat[,paste0('rhythm.', i)])}

used_dat.x <- subdat[,-1]
mice_dat <- mice(used_dat.x, m = 1, maxit = 5, meth = 'cart', seed = 123, printFlag = FALSE)
impute_dat.x <- mice:::complete(mice_dat, action = 1)

set.seed(0)
all_idx <- 1:nrow(subdat)

train_idx <- sample(all_idx, nrow(subdat) * 0.6)
valid_idx <- sample(all_idx[!all_idx %in% train_idx], nrow(subdat) * 0.2)
test_idx <- all_idx[!all_idx %in% c(train_idx, valid_idx)]

train_X <- impute_dat.x[train_idx,]
valid_X <- impute_dat.x[valid_idx,]
test_X <- impute_dat.x[test_idx,]

train_Y <- subdat[train_idx,"K"]
valid_Y <- subdat[valid_idx,"K"]
test_Y <- subdat[test_idx,"K"]

第一節:梯度提升機(2)

– 有別於隨機森林集成眾多深且獨立的樹模型,梯度提升機則是集成諸多淺且弱連續的樹模型,每個樹模型會以之前的樹模型為基礎去學習和精進,結果通常是難以擊敗的。

– 想想決策樹,他對於變項切點的選擇是如此的武斷,這在其他模型中也有類似的現象。

– 我們已經體驗到隨機森林的威力的,這就是「集成模型」的厲害,而隨機森林其實就是用Bagging的方法結合多個樹。

F01

第一節:梯度提升機(3)

– 有趣的是,這樣模型的組合將可以讓我們擁有一個非線性模型。

  1. 先建立第一個模型:\(y = F_1(x) + \epsilon_1\)

  2. 根據第一個模型的殘差建立第二個模型:\(\epsilon_1 = F_2(x) + \epsilon_2\)

  3. 先將上面兩個式子合併,產生一個整合模型,並取得整合模型的殘差:\(y = F_1(x) + \eta F_2(x) + \tilde{\epsilon_2}\)

  4. 再根據這個殘差建立第三個模型:\(\tilde{\epsilon_2} = F_3(x) + \epsilon_3\)

  5. 再將三個式子合併,產生一個整合模型,並取得整合模型的殘差:\(y = F_1(x) + \eta F_2(x) + \eta F_3(x) + \tilde{\epsilon_3}\)

  6. 再根據這個殘差建立第四個模型:\(\tilde{\epsilon_3} = F_4(x) + \epsilon_4\)

  7. 再將三個式子合併,產生一個整合模型,並取得整合模型的殘差:\(y = F_1(x) + \eta F_2(x) + \eta F_3(x) + \eta F_4(x) + \tilde{\epsilon_4}\)

  8. 依此類推…

第一節:梯度提升機(4)

library(party)

eta <- 0.3
tree.model_list <- list()

for (i in 1:10) {
  
  if (i == 1) {res_y <- train_Y}
  
  tree.model_list[[i]] <- ctree(formula = res_y ~ ., data = train_X, controls = ctree_control(mincriterion = 0.95, maxdepth = 2))
  
  for (j in 1:length(tree.model_list)) {
    now_y <- predict(tree.model_list[[j]], train_X)
    if (j == 1) {pred_y <- now_y} else {pred_y <- pred_y + eta * now_y}
  }
  
  res_y <- train_Y - pred_y

}
plot(tree.model_list[[1]], type = 'simple', main = 'Tree 1')

plot(tree.model_list[[2]], type = 'simple', main = 'Tree 2')

plot(tree.model_list[[3]], type = 'simple', main = 'Tree 3')

第一節:梯度提升機(5)

for (j in 1:length(tree.model_list)) {
  now_y <- predict(tree.model_list[[j]], valid_X)
  if (j == 1) {pred_valid <- now_y} else {pred_valid <- pred_valid + eta * now_y}
}

cor(pred_valid, valid_Y)
##            [,1]
## res_y 0.2865046
set.seed(0)

tree.model_list <- list()

for (i in 1:3) {
  
  new_train_idx <- sample(1:nrow(train_X), 0.7 * nrow(train_X))
  sub_train_X <- train_X[new_train_idx,]
  sub_train_Y <- train_Y[new_train_idx]
  
  tree.model_list[[i]] <- ctree(formula = sub_train_Y ~ ., data = sub_train_X,
                                controls = ctree_control(mincriterion = 0.95, maxdepth = 4))

}

pred_valid <- 0

for (i in 1:length(tree.model_list)) {
  
  prob_y <- predict(tree.model_list[[i]], valid_X)
  pred_valid <- pred_valid + prob_y / length(tree.model_list)
  
}

cor(pred_valid, valid_Y)
##                  [,1]
## sub_train_Y 0.2734035

第二節:使用套件做出梯度提升機(1)

– 這個套件實現的是「eXtreme Gradient Boosting」,跟傳統的梯度提升機不完全一樣,但原理是類似的。

– 與套件「glmnet」類似,他一樣必須吃矩陣格式的資料,並且還要轉換成自己的格式。

library(xgboost)

train_X_mat <- model.matrix(~ ., data = train_X)
xgb.data_train <- xgb.DMatrix(data = train_X_mat[,-1], label = train_Y)
xgb_fit <- xgb.train(data = xgb.data_train, max.depth = 2, eta = 0.5, nthread = 2, nrounds = 100, objective = "reg:linear", verbose = FALSE)
valid_X_mat <- model.matrix(~ ., data = valid_X)
pred_valid <- predict(xgb_fit, valid_X_mat[,-1])

cor(pred_valid, valid_Y)
## [1] 0.3696494
test_X_mat <- model.matrix(~ ., data = test_X)
pred_test <- predict(xgb_fit, test_X_mat[,-1])

cor(pred_test, test_Y)
## [1] 0.3582173

第二節:使用套件做出梯度提升機(2)

– 這個函數最有趣的是可以設置提早停止,以避免過度擬合。

train_X_mat <- model.matrix(~ ., data = train_X)
xgb.data_train <- xgb.DMatrix(data = train_X_mat[,-1], label = train_Y)

valid_X_mat <- model.matrix(~ ., data = valid_X)
xgb.data_valid <- xgb.DMatrix(data = valid_X_mat[,-1], label = valid_Y)

xgb_fit <-  xgb.train(data = xgb.data_train, watchlist = list(eval = xgb.data_valid),
                      early_stopping_rounds = 10, eval_metric = 'rmse', verbose = FALSE,
                      max.depth = 2, eta = 0.5, nthread = 2, nrounds = 100, objective = "reg:linear")
valid_X_mat <- model.matrix(~ ., data = valid_X)
pred_valid <- predict(xgb_fit, valid_X_mat[,-1])

cor(pred_valid, valid_Y)
## [1] 0.3894055
test_X_mat <- model.matrix(~ ., data = test_X)
pred_test <- predict(xgb_fit, test_X_mat[,-1])

cor(pred_test, test_Y)
## [1] 0.3795727

第二節:使用套件做出梯度提升機(3)

– 如果你很想看每棵樹是怎樣分析的,你可以用函數「xgb.dump」來看看。

importance_matrix <- xgb.importance(model = xgb_fit)
print(importance_matrix)
##        Feature        Gain       Cover  Frequency
##  1:        AGE 0.224172044 0.107929678 0.11578947
##  2:       QRSd 0.135477011 0.120020943 0.12631579
##  3:        QTc 0.131120408 0.187512783 0.15789474
##  4:     Axes_T 0.118617802 0.126606510 0.10526316
##  5:       Rate 0.093157932 0.062387002 0.10526316
##  6:         QT 0.078398695 0.080425730 0.07368421
##  7:   Axes_QRS 0.055801122 0.078012386 0.07368421
##  8:     Axes_P 0.032946354 0.005432070 0.04210526
##  9: rhythm.201 0.028381946 0.011175012 0.01052632
## 10: rhythm.121 0.016130145 0.031250767 0.02105263
## 11:         PR 0.015242987 0.047473351 0.05263158
## 12: rhythm.111 0.012960160 0.031250767 0.02105263
## 13: rhythm.191 0.011586085 0.029262825 0.02105263
## 14:  rhythm.81 0.011450748 0.003272332 0.01052632
## 15:  rhythm.21 0.010607743 0.013596538 0.01052632
## 16:  rhythm.71 0.006287813 0.002830567 0.01052632
## 17: rhythm.291 0.004935006 0.015625383 0.01052632
## 18: rhythm.181 0.004451274 0.015625383 0.01052632
## 19:  rhythm.41 0.004338230 0.015625383 0.01052632
## 20:  rhythm.11 0.003936495 0.014684588 0.01052632
xgb.plot.importance(importance_matrix = importance_matrix)

練習1:與隨機森林和SVM比較

– 我們就以二元分類預測LVD為例,請你分別在3個模型中都試著調整一下參數,你必須找出一組參數讓驗證組準確度最高。

library(mice)

subdat <- dat[!(dat[,"LVD"] %in% NA), c(-1, -2, -4, -5)]

subdat[,'LVD'] <- as.factor(subdat[,'LVD'])
subdat[,'GENDER'] <- as.factor(subdat[,'GENDER'])
for (i in 1:31) {subdat[,paste0('rhythm.', i)] <- as.factor(subdat[,paste0('rhythm.', i)])}

used_dat.x <- subdat[,-1]
mice_dat <- mice(used_dat.x, m = 1, maxit = 5, meth = 'cart', seed = 123, printFlag = FALSE)
impute_dat.x <- mice:::complete(mice_dat, action = 1)

set.seed(0)
all_idx <- 1:nrow(subdat)

train_idx <- sample(all_idx, nrow(subdat) * 0.6)
valid_idx <- sample(all_idx[!all_idx %in% train_idx], nrow(subdat) * 0.2)
test_idx <- all_idx[!all_idx %in% c(train_idx, valid_idx)]

train_X <- impute_dat.x[train_idx,]
valid_X <- impute_dat.x[valid_idx,]
test_X <- impute_dat.x[test_idx,]

train_Y <- subdat[train_idx,"LVD"]
valid_Y <- subdat[valid_idx,"LVD"]
test_Y <- subdat[test_idx,"LVD"]

練習1答案(1)

library(pROC)
library(randomForest)

result.rf <- data.frame(model = 'randomForest',
                        replace = rep(c(TRUE, FALSE), each = 4),
                        sampsize = rep(c(600, 800), each = 2),
                        nodesize = c(20, 50),
                        valid_auc = NA)

for (i in 1:nrow(result.rf)) {
  
  rf.model <- randomForest(formula = train_Y ~ ., data = train_X, ntree = 50,
                           replace = result.rf[i,'replace'],
                           sampsize = result.rf[i,'sampsize'],
                           nodesize =  result.rf[i,'nodesize'])
  
  pred.y <- predict(rf.model, valid_X, type = 'prob')[,2]
  roc_valid <- roc(valid_Y ~ pred.y)
  result.rf[i,'valid_auc'] <- roc_valid[['auc']]
  
}

result.rf
##          model replace sampsize nodesize valid_auc
## 1 randomForest    TRUE      600       20 0.7846612
## 2 randomForest    TRUE      600       50 0.7913805
## 3 randomForest    TRUE      800       20 0.7850618
## 4 randomForest    TRUE      800       50 0.7948031
## 5 randomForest   FALSE      600       20 0.7936813
## 6 randomForest   FALSE      600       50 0.7901557
## 7 randomForest   FALSE      800       20 0.7793613
## 8 randomForest   FALSE      800       50 0.7864812

練習1答案(2)

library(pROC)
library(e1071)

result.svm <- data.frame(model = 'svm',
                         cost = rep(2^(-3:-1), each = 3),
                         gamma = 2^(-2:0),
                         valid_auc = NA)

for (i in 1:nrow(result.svm)) {
  
  svm.model <- svm( train_Y ~ ., data = train_X, kernel = "radial", type = 'C-classification', cost = result.svm[i,'cost'], gamma = result.svm[i,'gamma'])
  pred_valid <- predict(svm.model, valid_X, decision.values = TRUE)
  pred.y <- attr(pred_valid, "decision.values")
  roc_valid <- roc(valid_Y ~ pred.y)
  result.svm[i,'valid_auc'] <- roc_valid[['auc']]
  
}

result.svm
##   model  cost gamma valid_auc
## 1   svm 0.125  0.25 0.7896978
## 2   svm 0.125  0.50 0.7935897
## 3   svm 0.125  1.00 0.7662546
## 4   svm 0.250  0.25 0.7893773
## 5   svm 0.250  0.50 0.7935897
## 6   svm 0.250  1.00 0.7664148
## 7   svm 0.500  0.25 0.7904304
## 8   svm 0.500  0.50 0.7831273
## 9   svm 0.500  1.00 0.7663462

練習1答案(3)

library(pROC)
library(xgboost)

result.xgb <- data.frame(model = 'xgb',
                         max.depth = rep(2:4, each = 3),
                         eta = c(0.3, 0.5, 0.7),
                         valid_auc = NA)

train_X_mat <- model.matrix(~ ., data = train_X)
xgb.data_train <- xgb.DMatrix(data = train_X_mat[,-1], label = as.integer(train_Y) - 1L)

valid_X_mat <- model.matrix(~ ., data = valid_X)
xgb.data_valid <- xgb.DMatrix(data = valid_X_mat[,-1], label = as.integer(valid_Y) - 1L)

for (i in 1:nrow(result.xgb)) {
  
  xgb_fit <-  xgb.train(data = xgb.data_train, watchlist = list(eval = xgb.data_valid),
                        early_stopping_rounds = 10, eval_metric = 'auc', verbose = FALSE,
                        max.depth = result.xgb[i,'max.depth'], eta = result.xgb[i,'eta'],
                        nthread = 2, nrounds = 100, objective = "binary:logistic")
  
  pred_valid <- predict(xgb_fit, valid_X_mat[,-1])
  roc_valid <- roc(valid_Y ~ pred_valid)
  result.xgb[i,'valid_auc'] <- roc_valid[['auc']]
  
}

result.xgb
##   model max.depth eta valid_auc
## 1   xgb         2 0.3 0.8063416
## 2   xgb         2 0.5 0.8026213
## 3   xgb         2 0.7 0.7923191
## 4   xgb         3 0.3 0.7964515
## 5   xgb         3 0.5 0.8004922
## 6   xgb         3 0.7 0.7951007
## 7   xgb         4 0.3 0.7904304
## 8   xgb         4 0.5 0.7788690
## 9   xgb         4 0.7 0.7612179

第三節:集成模型(1)

– 讓我們在測試組上看看它的準確度:

best_pos.xgb <- which.max(result.xgb[,'valid_auc'])

best_xgb_fit <-  xgb.train(data = xgb.data_train, watchlist = list(eval = xgb.data_valid),
                           early_stopping_rounds = 10, eval_metric = 'auc', verbose = FALSE,
                           max.depth = result.xgb[best_pos.xgb,'max.depth'], eta = result.xgb[best_pos.xgb,'eta'],
                           nthread = 2, nrounds = 100, objective = "binary:logistic")

test_X_mat <- model.matrix(~ ., data = test_X)
pred_test <- predict(best_xgb_fit, test_X_mat[,-1])
roc_test <- roc(test_Y ~ pred_test)

plot(roc_test)
text(0.5, 0.5, paste0('AUC = ', formatC(roc_test[['auc']], 4, format = 'f')), col = 'red')

– 我們想想為什麼隨機森林能夠比決策樹好這麼多?原因是他綜合了多個模型的預測結果,從而讓我們獲得更準確的預測!

第三節:集成模型(2)

– 讓我們先拿到最佳的極限梯度提升機預測值:

best_pred.xgb <- predict(best_xgb_fit, valid_X_mat[,-1])

– 讓我們再拿到最佳的隨機森林預測值:

best_pos.rf <- which.max(result.rf[,'valid_auc'])

best_rf.model <- randomForest(formula = train_Y ~ ., data = train_X, ntree = 50,
                              replace = result.rf[best_pos.rf,'replace'],
                              sampsize = result.rf[best_pos.rf,'sampsize'],
                              nodesize =  result.rf[best_pos.rf,'nodesize'])
  
best_pred.rf <- predict(best_rf.model, valid_X, type = 'prob')[,2]

– 讓我們再拿到最佳的支持向量機預測值:

best_pos.svm <- which.max(result.svm[,'valid_auc'])

best_svm.model <-  svm(train_Y ~ ., data = train_X, kernel = "radial",
                       type = 'C-classification', cost = result.svm[best_pos.svm,'cost'], gamma = result.svm[best_pos.svm,'gamma'])


pred_valid <- predict(best_svm.model, valid_X, decision.values = TRUE)
best_pred.svm <- attr(pred_valid, "decision.values")

第三節:集成模型(3)

– 我們就不再重新搜索最佳參數了(這樣還需要一個獨立的組別),直接使用剛剛找到的最佳參數。

new_valid_X <- valid_X
new_valid_X[,'pred.rf'] <- best_pred.rf
new_valid_X[,'pred.svm'] <- best_pred.svm
new_valid_X[,'pred.xgb'] <- best_pred.xgb

new_test_X <- test_X
new_test_X[,'pred.rf'] <- predict(best_rf.model, test_X, type = 'prob')[,2]
new_test_X[,'pred.svm'] <- attr(predict(best_svm.model, test_X, decision.values = TRUE), "decision.values")
new_test_X[,'pred.xgb'] <- predict(xgb_fit, test_X_mat[,-1])

valid_rf.model <- randomForest(formula = valid_Y ~ ., data = new_valid_X, ntree = 50,
                               replace = result.rf[best_pos.rf,'replace'],
                               nodesize =  result.rf[best_pos.rf,'nodesize'])
par(mfrow = c(2, 2))

pred_test <- predict(best_rf.model, test_X, type = 'prob')[,2]
roc_test <- roc(test_Y ~ pred_test)
plot(roc_test, main = 'Random forest')
text(0.5, 0.5, paste0('AUC = ', formatC(roc_test[['auc']], 4, format = 'f')), col = 'red')

pred_test <- attr(predict(best_svm.model, test_X, decision.values = TRUE), "decision.values")
roc_test <- roc(test_Y ~ pred_test)
plot(roc_test, main = 'Support vector machine')
text(0.5, 0.5, paste0('AUC = ', formatC(roc_test[['auc']], 4, format = 'f')), col = 'red')

pred_test <- predict(best_xgb_fit, test_X_mat[,-1])
roc_test <- roc(test_Y ~ pred_test)
plot(roc_test, main = 'EXtreme gradient boosting')
text(0.5, 0.5, paste0('AUC = ', formatC(roc_test[['auc']], 4, format = 'f')), col = 'red')

pred_test <- predict(valid_rf.model, new_test_X, type = 'prob')[,2]
roc_test <- roc(test_Y ~ pred_test)
plot(roc_test, main = 'Ensemble learning')
text(0.5, 0.5, paste0('AUC = ', formatC(roc_test[['auc']], 4, format = 'f')), col = 'red')

期末報告:親手完成一次的完整實驗

– 讓我們複習一下整個分析流程:

  1. 使用多重插補法進行資料插補,而這裡可以插補出M個資料集

  2. 將資料分成「訓練組」、「驗證組」、「測試組」,其中「訓練組」與「驗證組」是可以合併的,並且用K-倍交叉驗證,這樣可以產生K個分割法

  3. 你可以選擇是不是要篩選變項,如果你想篩選變項的話,你可以使用「統計分析」、「隨機森林」、「梯度提升機」等方法篩選變項。

  4. 你也可以選擇要不要篩選樣本,我們這門課沒有特別教這個技術,但原則上我們是可以做資料篩選的(但不要在測試組上篩選,測試組要盡可能符合未來的使用情境)

  5. 接著你能使用N個算法進行分析,光我們教過的機器學習算法就有8種(線性判別分析、樸素貝葉斯分類、k-近鄰演算法、彈性網路、多層感知機、決策樹、隨機森林、梯度提升機),而每種模型都有一堆參數等著你調整。

  6. 另外,其實不同套件在實現同一個算法時結果多少有點不同,舉例來說你也可以用多個套件同時實現梯度提升機,從而得到更多模型。

  7. 最後,統合這麼多個模型後,你將能夠獲得一個最終的預測結果。

– 這裡有一些小技巧,假設你的樣本實在太大,那在調參數的時候其實可以隨機抽一個較小的樣本,並將最佳參數組合應用到大樣本之上。

課程總結

– 這些瘋狂的調參數技術是必要的技能,你會發現調參數帶來的效能提升是非常明顯的,這比起重新發展一個演算法來的容易得多。

– 但要特別注意的是,一定要小心測試組的資料,原則上在一個研究中,我們還是只能夠在測試組上進行一次運算,我們不能因為測試組的準確度而選擇不同的模型,這樣會「高估」未來應用的準確度。

– 不過發展演算法仍然是很重要的,當大家都有充足的調參數能力後,下一步就是看看在算法上能不能突破。這門課所介紹每個方法的原理都能帶給我們未來在發展新算法時的一些靈感,充分了解過去所有算法是發展新算法的重要支持。

– 在數值分析領域上,你會發現人工神經網路(多層感知機)的能力似乎遠不如這些機器學習算法,但它在2012年後被發現於非結構化數據的分析上有非常強大的能力,從而帶領這第三次人工智為革命的到來。

– 雖然我們已經介紹過多層感知機了,但深度學習的相關知識尚未有完整且深入的介紹,下學期開始我們將開始介紹人工神經網路的課程,我們將逐步介紹人工神經網路的原理、特殊型態、近代發展,從而了解為什麼人工神經網路打敗了支持向量機