r语言随机森林回归预测_从零实现回归随机森林
一、前言
回歸隨機森林作為一種機器學(xué)習(xí)和數(shù)據(jù)分析領(lǐng)域常用且有效的算法，對其原理和代碼實現(xiàn)過程的掌握是非常有必要的。為此，本文將著重介紹從零開始實現(xiàn)回歸隨機森林的過程，對于隨機森林和決策樹的相關(guān)理論原理將不做太深入的描述。本文的目的只是為了演示回歸隨機森林主要功能的具體實現(xiàn)過程，在實現(xiàn)過程中不會考慮代碼性能，會更加注重代碼可讀性。
實現(xiàn)語言：Python
依賴：pandas, numpy
二、原理介紹
隨機森林屬于Bagging類算法，而Bagging 又屬于集成學(xué)習(xí)一種方法（集成學(xué)習(xí)方法大致分為Boosting和Bagging方法，兩個方法的不同請參考[10]），集成學(xué)習(xí)的大致思路是訓(xùn)練多個弱模型打包起來組成一個強模型，強模型的性能要比單個弱模型好很多（三個臭皮匠頂一個諸葛亮。注意：這里的弱和強是相對的），其中的弱模型可以是決策樹、SVM等模型，在隨機森林中，弱模型選用決策樹。
在訓(xùn)練階段，隨機森林使用bootstrap采樣從輸入訓(xùn)練數(shù)據(jù)集中采集多個不同的子訓(xùn)練數(shù)據(jù)集來依次訓(xùn)練多個不同決策樹；在預(yù)測階段，隨機森林將內(nèi)部多個決策樹的預(yù)測結(jié)果取平均得到最終的結(jié)果。本文主要介紹回歸隨機森林從零實現(xiàn)的過程，實現(xiàn)的RFR(回歸隨機森林)具有以下功能：
- 模型訓(xùn)練
- 模型數(shù)據(jù)預(yù)測
- 計算feature importance
2.1 模型訓(xùn)練
2.1.1 基礎(chǔ)理論
本文實現(xiàn)的RFR是將多個二叉決策樹（即CART，這也是sklearn,spark內(nèi)部實現(xiàn)的模型）打包組合而成的，訓(xùn)練RFR便是訓(xùn)練多個二叉決策樹。在訓(xùn)練二叉決策樹模型的時候需要考慮怎樣選擇切分變量(特征）、切分點以及怎樣衡量一個切分變量、切分點的好壞。針對于切分變量和切分點的選擇，本實現(xiàn)采用窮舉法，即遍歷每個特征和每個特征的所有取值，最后從中找出最好的切分變量和切分點；針對于切分變量和切分點的好壞，一般以切分后節(jié)點的不純度來衡量，即各個子節(jié)點不純度的加權(quán)和
，其計算公式如下： (2-1)其中，
為某一個切分變量， 為切分變量的一個切分值， 分別為切分后左子節(jié)點的訓(xùn)練樣本個數(shù)、右子節(jié)點的訓(xùn)練樣本個數(shù)以及當(dāng)前節(jié)點所有訓(xùn)練樣本個數(shù)， 分為左右子節(jié)點的訓(xùn)練樣本集合， 為衡量節(jié)點不純度的函數(shù)(impurity function/criterion)，分類和回歸任務(wù)一般采用不同的不純度函數(shù)，在分類和回歸常采用的不純度函數(shù)有以下四種：圖 1. 四種不同的impurity function- 備注1：
Gini不純度只適用于分類任務(wù)，
為當(dāng)前節(jié)點的訓(xùn)練樣本集合， 為當(dāng)前節(jié)點訓(xùn)練樣本中目標變量取值k( )出現(xiàn)的概率。- 備注2：
信息熵只適用于分類任務(wù)， 其他同備注1。
- 備注3
MSE只適用于回歸任務(wù)，
為當(dāng)前節(jié)點樣本目標變量的平均值。- 備注4
MAE只適用于回歸任務(wù)，
為當(dāng)前節(jié)點樣本目標變量的平均值。- 備注5
在sklearn內(nèi)部，DecisionTreeClassifier, RandomForestClassifier等基于決策樹的分類模型默認使用'gini'作為impurity function，也可通過criterion參數(shù)指定為'entropy' ；而DecisionTreeRegressor, RandomForestRegressor等基于決策樹的回歸模型默認使用'mse'作為impurity function，也可通過criterion參數(shù)指定為'mae'。
決策樹中某一節(jié)點的訓(xùn)練過程在數(shù)學(xué)等價于下面優(yōu)化問題：
(2-2)即尋找
最小的切分變量和切分點。在本文實現(xiàn)的回歸隨機森林中，
選用MSE，即針對某一切分點 (2-3)2.1.2 訓(xùn)練過程
RFR訓(xùn)練過程示意圖如圖2所示。
圖 2. RFR訓(xùn)練示意圖說明：
- bootstrap[1]是對輸入訓(xùn)練樣本集合 進行N次放回重復(fù)抽樣得到樣本集合 (N一般等于 的大小[2]，具體可以查看sklearn的實現(xiàn)代碼)，這暗示了 中的某些樣本會在 中重復(fù)出現(xiàn)多次。在對 進行bootstrap采樣時， 中每個樣本會以63.3%的概率被抽中。為什么是63.3%？考慮N次放回重復(fù)抽樣，每次抽樣每個樣本被選中的概率為 ，進行N次抽樣，被選中概率為：
具體可以參考[3-4]。
- subsample會根據(jù)輸入?yún)?shù)sample_sz的大小從 中采樣sample_sz個樣本組成subsample樣本集合 (在sklearn中，subsample大小就是輸入訓(xùn)練樣本集合 的大小N，不可以通過輸入?yún)?shù)指定subsample的大小，詳見這里)。
- 為什么要進行bootstrap? 集成學(xué)習(xí)（包括隨機森林）的目的是為了使用多個不同的子模型來增加最終模型預(yù)測結(jié)果的魯棒性和穩(wěn)定性(即減小方差)，如果多個子模型都采用同樣的數(shù)據(jù)集訓(xùn)練，那么訓(xùn)練得出的子模型都是相同的，集成學(xué)習(xí)將變得沒有意義，為了能從原輸入訓(xùn)練數(shù)據(jù)集得到多個不同的數(shù)據(jù)集，需要使用bootstrap來從原數(shù)據(jù)集采樣出不同的子數(shù)據(jù)集來訓(xùn)練子模型。
圖2中n個回歸決策樹的訓(xùn)練過程如圖3所示。
圖 3. 回歸決策樹訓(xùn)練過程2.2 模型預(yù)測
RFR的預(yù)測結(jié)果是由內(nèi)部所有二叉決策樹的預(yù)測結(jié)果取平均值得到的。二叉決策樹的預(yù)測過程主要分為以下步驟：
2.3 計算feature importance
特征的重要性表示特征對預(yù)測結(jié)果影響程度，某一特征重要性越大，表明該特征對預(yù)測結(jié)果的影響越大，重要性越小，表明該特征對預(yù)測結(jié)果越小。RFR中某一特征的重要性，是該特征在內(nèi)部所有決策樹重要性的平均值，而在決策樹中，計算某一個特征的重要性可以采用以下三種方法：
方法1
步驟如下：
- 使用訓(xùn)練數(shù)據(jù)訓(xùn)練模型；
- 計算訓(xùn)練數(shù)據(jù)在模型上依據(jù)某一metric的score，記為 （在回歸中，可以選用r2）；
- 遍歷訓(xùn)練數(shù)據(jù)集中的每一個feature，每次在原訓(xùn)練數(shù)據(jù)集的基礎(chǔ)上將對應(yīng)的feature進行shuffle操作，然后使用模型得到shuffle后數(shù)據(jù)集的score，記為 ，最后通過 計算出第 個feature的重要性。
方法2
如果一個feature很重要，那么從數(shù)據(jù)集中去除該feature后，模型的prediction error會增大；相反，如果一個feature不重要，那么從數(shù)據(jù)集中去除后，模型的prediction error不會變化太大。該方法的步驟如下：
- 使用原數(shù)據(jù)集訓(xùn)練模型，并評估模型在訓(xùn)練數(shù)據(jù)集上的prediction error，記為 ;
- 遍歷訓(xùn)練數(shù)據(jù)集中的每一個feature，每次從原訓(xùn)練數(shù)據(jù)集的基礎(chǔ)上去除該feature，然后使用得到數(shù)據(jù)集訓(xùn)練模型，最終計算出prediction error，記為 ，最后通過 計算出第 個feature的重要性。
該方法和方法1有點類似，但該方法在計算每個feature的重要性的時候都需要重新訓(xùn)練模型，會增加額外的計算量，在實際應(yīng)用一般不采用該方法。
方法3
該方法是sklearn內(nèi)部樹模型計算feature重要性時采用的方法。即某一節(jié)點
的重要性為 (2-5)其中，
分別為節(jié)點 以及其左右子節(jié)點中訓(xùn)練樣本個數(shù)與總訓(xùn)練樣本數(shù)目的比例， 分為為節(jié)點 以及其左右子節(jié)點的不純度。知道每一個節(jié)點的重要性之后，即可通過公式(2-6)得出某一feature的重要性。 (2-6)為了使所有feature的重要性加起來等于1，需要每一feature的重要性進行normalization，即公式(2-7)
(2-7)方法3在sklearn中的實現(xiàn)，請查看_tree.pyx。
在本文實現(xiàn)的RFR中，同時實現(xiàn)了方法1和方法3。
三、回歸隨機森林實現(xiàn)
3.1 代碼
代碼有點長，不想看的可以直接跳過。
from multiprocessing import Pool, cpu_countclass RFRegressor(object):"""Random Forest Regressor"""def __init__(self, n_trees, sample_sz, min_leaf_sz=5, n_jobs=None, max_depth=None):self._n_trees = n_treesself._sample_sz = sample_szself._min_leaf_sz = min_leaf_szself._n_jobs = n_jobsself._max_depth = max_depthself._trees = [self._create_tree() for i in range(self._n_trees)]def _get_sample_data(self, bootstrap=True):"""Generate training data for each underlying decision treeParameters----------bootstrap: boolean value, True/FalseThe default value is True, it would bootstrap sample from input training data. If False, the exclusive sampling willbe performed.Returns-------idxs: array-like objectReturn the indices of sampled data from input training data"""if bootstrap:idxs = np.random.choice(len(self._X), self._sample_sz)else:idxs = np.random.permutation(len(self._X))[:self._sample_sz]return idxsdef _create_tree(self):"""Build decision treeeReturns-------dtree : DTreeRegressor object"""return DTreeRegressor(self._min_leaf_sz, self._max_depth)def _single_tree_fit(self, tree):"""Fit the single underlying decision treeParameters----------tree : DTreeRegressor objectReturns-------tree : DTreeRegressor object"""sample_idxs = self._get_sample_data()return tree.fit(self._X.iloc[sample_idxs], self._y[sample_idxs])def fit(self, x, y):"""Train a forest regressor of trees from the training set(x, y)Parameters----------x : DataFrame,The training input samples.y : Series or array-like objectThe target values."""self._X = xself._y = yif self._n_jobs:self._trees = self._parallel(self._trees, self._single_tree_fit, self._n_jobs)else:for tree in self._trees:self._single_tree_fit(tree)def predict(self, x):"""Predict target values using trained modelParameters---------x : DataFrame or array-like objectinput samplesReturns-------ypreds : array-like objectpredicted target values"""all_tree_preds = np.stack([tree.predict(x) for tree in self._trees]) return np.mean(all_tree_preds, axis=0)def _parallel(self, trees, fn, n_jobs=1):"""Parallel jobs executionParameters----------trees : list-like objecta list-like object contains all underlying treesfn : function-like objectmap functionn_jobs : integerThe number of jobs.Returns-------result : list-like objecta list-like result object for each call of map function `fn`"""try:workers = cpu_count()except NotImplementedError:workers = 1if n_jobs:workders = n_jobspool = Pool(processes=workers)result = pool.map(fn, trees)return list(result)@propertydef feature_importances_(self):"""Calculate feature importanceReturns-------self._feature_importances : array-like objectthe importance score of each feature """if not hasattr(self, '_feature_importances'):norm_imp = np.zeros(len(self._X.columns))for tree in self._trees:t_imp = tree.calc_feature_importance()norm_imp = norm_imp + t_imp / np.sum(t_imp)self._feature_importances = norm_imp / self._n_treesreturn self._feature_importances@propertydef feature_importances_extra(self):"""Another method to calculate feature importance"""norm_imp = np.zeros(len(self._X.columns))for tree in self._trees:t_imp = tree.calc_feature_importance_extra()norm_imp = norm_imp + t_imp / np.sum(t_imp)norm_imp = norm_imp / self._n_treesimp = pd.DataFrame({'col':self._X.columns, 'imp':norm_imp}).sort_values('imp', ascending=False)return impclass DTreeRegressor(object):def __init__(self, min_leaf_sz, max_depth=None):self._min_leaf_sz = min_leaf_szself._split_point = 0self._split_col_idx = 0self._score = float('inf')self._sample_sz = 0self._left_child_tree = Noneself._right_child_tree = Noneself._feature_importances = []self._node_importance= 0if max_depth is not None:max_depth = max_depth - 1self._max_depth = max_depthdef fit(self, x, y):self._X = xself._y = yself._col_names = self._X.columnsself._feature_importances = np.zeros(len(self._col_names))self._sample_sz = len(self._X)self._val = np.mean(self._y)if self._max_depth is not None and self._max_depth < 2:return selfself._find_best_split()return selfdef _calc_mse_inpurity(self, y_squared_sum, y_sum, n_y):"""Calculate Mean Squared Error impurityThis is just the recursive version for calculating varianceParameters----------y_squared_sum: float or int , the sum of y squared y_sum: float or int , the sum of y valuen_y: int, the number of samplesReturns-------"""dev = (y_squared_sum / n_y) - (y_sum / n_y) ** 2return devdef _find_best_split(self):for col_idx in range(len(self._col_names)):self._find_col_best_split_point(col_idx)self._feature_importances[self._split_col_idx] = self._node_importanceif self.is_leaf:return left_child_sample_idxs = np.nonzero(self.split_col <= self.split_point)[0]right_child_sample_idxs = np.nonzero(self.split_col > self.split_point)[0]self._left_child_tree = (DTreeRegressor(self._min_leaf_sz, self._max_depth).fit(self._X.iloc[left_child_sample_idxs], self._y[left_child_sample_idxs]))self._right_child_tree = (DTreeRegressor(self._min_leaf_sz, self._max_depth).fit(self._X.iloc[right_child_sample_idxs], self._y[right_child_sample_idxs]))def _find_col_best_split_point(self, col_idx):x_col = self._X.values[:, col_idx]sorted_idxs = np.argsort(x_col)sorted_x_col = x_col[sorted_idxs]sorted_y = self._y[sorted_idxs]lchild_n_samples = 0lchild_y_sum = 0.0lchild_y_square_sum = 0.0rchild_n_samples = self._sample_szrchild_y_sum = sorted_y.sum()rchild_y_square_sum = (sorted_y ** 2).sum()node_y_sum = rchild_y_sumnode_y_square_sum = rchild_y_square_sum for i in range(0, self._sample_sz - self._min_leaf_sz):xi, yi = sorted_x_col[i], sorted_y[i]rchild_n_samples -= 1rchild_y_sum -= yirchild_y_square_sum -= (yi ** 2)lchild_n_samples += 1lchild_y_sum += yilchild_y_square_sum += (yi ** 2)if i < self._min_leaf_sz or xi == sorted_x_col[i+1]:continuelchild_impurity = self._calc_mse_inpurity(lchild_y_square_sum,lchild_y_sum, lchild_n_samples)rchild_impurity = self._calc_mse_inpurity(rchild_y_square_sum,rchild_y_sum, rchild_n_samples)split_score = (lchild_n_samples * lchild_impurity + rchild_n_samples * rchild_impurity) / self._sample_szif split_score < self._score:self._score = split_scoreself._split_point = xiself._split_col_idx = col_idxself._node_importance = (self._sample_sz * (self._calc_mse_inpurity(node_y_square_sum, node_y_sum, self._sample_sz) - split_score))def predict(self, x):if type(x) == pd.DataFrame:x = x.valuesreturn np.array([self._predict_row(row) for row in x])def _predict_row(self, row):if self.is_leaf:return self._valt = (self._left_child_tree if row[self._split_col_idx] <= self.split_point else self._right_child_tree)return t._predict_row(row)def __repr__(self):pr = f'sample: {self._sample_sz}, value: {self._val}rn'if not self.is_leaf:pr += f'split column: {self.split_name}, split point: {self.split_point}, score: {self._score} 'return prdef calc_feature_importance(self):if self.is_leaf:return self._feature_importancesreturn (self._feature_importances + self._left_child_tree.calc_feature_importance()+ self._right_child_tree.calc_feature_importance())def calc_feature_importance_extra(self):imp = []o_preds = self.predict(self._X.values)o_r2 = metrics.r2_score(self._y, o_preds)for col in self._col_names:tmp_x = self._X.copy()shuffle_col = tmp_x[col].valuesnp.random.shuffle(shuffle_col)tmp_x.loc[:,col] = shuffle_coltmp_preds = self.predict(tmp_x.values)tmp_r2 = metrics.r2_score(self._y, tmp_preds)imp.append((o_r2 - tmp_r2))imp = imp / np.sum(imp)return imp@propertydef split_name(self):return self._col_names[self._split_col_idx]@propertydef split_col(self):return self._X.iloc[:, self._split_col_idx]@propertydef is_leaf(self):return self._score == float('inf')@propertydef split_point(self):return self._split_point3.2 測試
訓(xùn)練與預(yù)測：
使用sklearn內(nèi)部的RandomForestRegressor的結(jié)果：
需要注意的是，上面兩次訓(xùn)練的樣本不同，所以和實際情況有點出入，但大體趨勢是對。
feature importance:
RFRegressor采用方法3計算的feature重要性：
RandomForestRegressor計算的feature重要性：
RFRegressor采用方法1計算的feature重要性：
參考與薦讀：
[1] Bootstrap Sample: Definition, Example
[2] Number of Samples per-Tree in a Random Forest
[3] Why on average does each bootstrap sample contain roughly two thirds of observations?
[4] Random Forests out-of-bag sample size
[5] Random Forests Leo Breiman and Adele Cutler
[6] Why is Random Forest with a single tree much better than a Decision Tree classifier?
[7] Winning the bias-variance tradeoff
[8] subsample-bootstrapping
[9] How can SciKit-Learn Random Forest sub sample size may be equal to original training data size?
[10] the difference between bagging and boosting
[11] The Mathematics of Decision Trees, Random Forest and Feature Importance in Scikit-learn and Spark
[12] Sklearn Decision Trees documentation
總結(jié)
以上是生活随笔為你收集整理的r语言随机森林回归预测_从零实现回归随机森林的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： Linux串口设置及编程（转）
- 下一篇： 不关注公众号可以获取openid吗_微信