機器自己摸索決策邏輯

在【Kaggle】房價預測模型在房產(chǎn)投資場景的應(yīng)用一文中，我提到隨機森林（Random Forest）算法模型具有良好的數(shù)據(jù)解釋性，本文就用從零開始寫該算法的方式，希望能徹底講清楚隨機森林的工作原理。

模型的另一種表達方式是函數(shù)：。所謂的數(shù)據(jù)解釋性，指的就是人能理解是如何通過x生成y的。換句話說，人能理解的內(nèi)容。比如，統(tǒng)計學最愛的邏輯回歸直線，它只有這2個變量，就是這條直線的斜率，是直線的截距，通過調(diào)整直線斜率和截距，就可以把兩類事物區(qū)分開來，或是預測y的區(qū)域。

隨機森林的基礎(chǔ)是決策樹，從結(jié)構(gòu)上看，它是決策樹的集合。決策樹是一種基于因果關(guān)系，通過邏輯論證推理來學習知識的模型。之所以說人能理解隨機森林，就是因為決策樹的方法論也是人理解和探索世界的方法論，只是現(xiàn)在人把這些工作交給機器來做了而已。

隨機和森林

隨機森林模型的工作原理并不復雜，它的核心思想都寫在它的腦門上了，就是“隨機”和“森林”。
，為數(shù)據(jù)和決策增加更多隨機概率。

• “森林”指的是由N棵決策樹集成（Ensembling）在一起，這樣可以降低單棵決策樹的預測方差（Variance）。
• “隨機”指的是為模型增加更多隨機性，主要體現(xiàn)在以下兩方面：
  • 數(shù)據(jù)采樣的隨機性。隨機森林采用自舉法（bootstrapping）的重采樣方法，讓每顆決策樹的采樣數(shù)據(jù)盡可能不同。
  • 決策的隨機性。決策樹的工作就是做決策，所謂的為決策增加隨機性，指的是增加“差決策”被選擇的概率。在做選擇時，有好不選而選差的，這聽起來很荒謬，但這就好像是小說《天龍八部》里虛竹破珍瓏棋局。虛竹恰恰是自殺式地下了一手任何會下棋的人都不會考慮的錯棋之后，才迎來真正的好棋，而那些一直在理所當然地下“好棋”的高手都掉進了珍瓏棋局的陷阱里，棋路越下越窄，最后被逼進絕路。

Let's get to work

對隨機森林的工作原理有了初步了解后，讓我們開始寫代碼吧。數(shù)據(jù)集還是使用【Kaggle】房價預測模型在房產(chǎn)投資場景的應(yīng)用的數(shù)據(jù)，為了驗證模型的準確性，用sklearn的隨機森林模型（RandomForestRegressor）作為對標模型。

首先，為了簡化模型，我只用"OverallQual"和"GrLivArea"這兩列數(shù)據(jù)，也就是說，模型現(xiàn)在只需要處理2個特征。

x = train_df[['OverallQual', 'GrLivArea']].copy()
trn_x, val_x, trn_y, val_y = train_test_split(x, y, test_size=0.3)
n = 1000
trn_x, trn_y = trn_x[:n], trn_y[:n]

隨機森林模型

class RandomForest():
  def __init__(self, x, y, nr_trees, sample_sz=1000):
    self.x, self.y, self.nr_trees, self.sample_sz = x, y, nr_trees, sample_sz
    self.trees = [self.create_tree() for _ in range(nr_trees)]
  
  def create_tree(self):
    idxs = np.random.permutation(len(self.x))[:self.sample_sz]
    return DecisionTree(self.x.iloc[idxs], self.y[idxs], min_leaf=5)

可以看到，隨機森林模型參數(shù)有兩個：nr_trees（決策樹數(shù)）和sample_sz（樣本采樣數(shù)）。前文提到，決策樹會采用重采樣（bootstrapping）來保證采樣的隨機性，然而我們的模型并沒有采用重采樣，而是采用隨機不重復采樣："idxs = np.random.permutation(len(self.x))[:self.sample_sz]"。

之所以不使用重采樣，是因為隨機森林是20多年前發(fā)明的，重采樣只是當時小數(shù)據(jù)時代的不得已之舉，現(xiàn)在是大數(shù)據(jù)時代，已經(jīng)不存在數(shù)據(jù)不足的問題，隨機不重復采樣才是更好之選。

這里將采樣數(shù)設(shè)置為1000（sample_sz=1000），是為了保證我們的模型和對標模型沒有數(shù)據(jù)集的差異（訓練集大小就是1000）。

決策樹模型

決策樹是一個根據(jù)事物屬性對事物分類的樹形結(jié)構(gòu)。比如病人小李去醫(yī)院看病，醫(yī)生會先問“哪里不舒服啦”，小李說“體乏、流鼻涕、咳嗽”，這時醫(yī)生會讓小李先去量體溫并做血常規(guī)檢測，因為他判斷小李很可能得了感冒，但也不排除是流感甚至肺炎的可能。如何判斷病人具體得到那種病呢，醫(yī)生這時候就會用到一個決策樹：體溫低于38.2°，血常規(guī)正常，那就是一般感冒；體溫高于38.2°，病情超過1周了，那可能是流感，加上血常規(guī)某項指標低于標準值，那還需要讓小李再去做另一項肺部檢查才能排除他是否得了肺炎。

醫(yī)生之所以能根據(jù)體溫和檢查報告來判斷病人得了什么病，是因為他學了幾年的醫(yī)學知識，從已知的病理出發(fā)：感冒、流感和肺炎都會產(chǎn)生小李描述的癥狀，因此先假設(shè)小李得了這三種病，接著收集數(shù)據(jù)來驗證假設(shè)。

決策樹則是采用逆向?qū)W習的方法，從病例的已知結(jié)果出發(fā)逆向找出因果鏈條。雖然兩者的學習方法不同，但結(jié)果是相同的，都是找出因果關(guān)系。

class DecisionTree():
  def __init__(self, x, y, idxs=None, min_leaf=1):
    if idxs is None: idxs = np.arange(len(x))
    self.x, self.y, self.min_leaf, self.idxs = x.reset_index(drop=True), y, min_leaf, idxs
    self.n, self.c = len(idxs), x.shape[1]
    self.val = self.y[idxs].mean()
    self.score = float('inf')
    self.find_split()
  
  def find_split(self):
    for c in range(self.c):
      self.find_better_split(c)
  
  def find_better_split(self, var_idx):
    x, y = self.x.iloc[self.idxs, var_idx].values, self.y[self.idxs]
    for i in range(self.n - self.min_leaf):
      l_idxs = (x <= x[i])
      r_idxs = x > x[i]
      if r_idxs.sum() == 0: continue
      l_std = y[l_idxs].std()
      r_std = y[r_idxs].std()
      score = l_std * l_idxs.sum() + r_std * r_idxs.sum()
      if score < self.score:
        self.score, self.var_idx, self.split_val = score, var_idx, x[i]
        
  @property
  def split_var(self): return self.x.columns[self.var_idx]
  
  @property
  def is_leaf(self): return self.score == float('inf')
  
  def __repr__(self):
    txt = f'samples: {self.n}, value: {self.val}, score: {self.score}, '
    if not self.is_leaf:
      txt += f'split variable: {self.split_var}, split value: {self.split_val}'
    return txt

find_better_split()是決策樹進行決策的函數(shù)，var_idx是column id，首先需要拿出所有樣本的同一列作為x，以小李看病例子為例，x表示1000個病人的體溫（或血常規(guī)檢查報告），y表示這1000個病人對應(yīng)的病癥。

決策樹是逆向?qū)ふ乙蚬湕l的，然而樹中的鏈條有很多，假設(shè)鏈條上有3個節(jié)點，每個節(jié)點有True和False兩個分叉，能形成條鏈，機器是如何判斷哪條鏈才是正確的因果關(guān)系鏈呢？

其實機器并沒有更好的辦法，它只能一個個試。比如它先以樣本1的體溫為參照點，把所有體溫小于等于樣本1的樣本歸為一類，把所有體溫大于樣本2的樣本歸為另一類。因為已知樣本的真實分類，因此可以統(tǒng)計這種分類方法的準確率水平，用變量score表示。用同樣的方法，再拿樣本2、樣本3......樣本1000作為參照點，最終選出score值最小的那個樣本，它的體溫就是最準確的分類特征，例如體溫小于等于就是普通感冒，大于這個溫度就是流感或肺炎。

score稱為“混雜度”，它表示樣本里分類的純度，score越小表示分類越準。因為你不可能只經(jīng)過一次分類就能把樣本都分完整了，因此，每次分類，或者說決策的過程，就是減少score值的過程。可以把計算score的函數(shù)理解為損失函數(shù)。

隨機森林提供了很多計算score的函數(shù)，如mse（默認值）、cross entropy、gini。MSE的平方根就是RMSE（root mean square error），它指的是樣本觀察值與均值的平均距離，值越小表示模型效果越好。RMSE的另一種近似的表達是標準差，它也是觀察值與均值的距離，值越小表示模型越好。相比RMSE，標準差的計算量要小，因此這里用標準差來計算score。

min_leaf可以用來調(diào)控參與決策的樣本數(shù)，該值越大，表示參與決策的樣本數(shù)越少，決策數(shù)也會跟著減少。在這里它的默認是5，但你不要以為只是少了5個樣本而已。因為決策樹需要做大量的決策（分叉），如果每次決策都少5個樣本，那結(jié)果會少許多決策。min_leaf=2，每顆決策樹的決策數(shù)是，min_leaf每增加一倍，決策數(shù)差不多就會減少一半。

既然決策樹的工作就是不斷決策，那么這個決策循環(huán)什么時候結(jié)束呢？當無法再對現(xiàn)在的樣本分類的時候，就是決策樹的終點，它們也稱為葉節(jié)點，可以用is_leaf來查詢是否是葉節(jié)點。葉子節(jié)點的score == inf，即無限大，表示它們不參與決策。

單次決策

rf = RandomForest(trn_x, trn_y, 1)
print(rf.trees[0])

samples: 1000, value: 12.030368441695023, score: 288.68995409377146,
split variable: OverallQual, split value: 6

Figure 1

前面用了很長的篇幅來解釋決策樹的工作原理，接了來就是驗證決策樹模型的表現(xiàn)了。Figure1是對標模型的決策樹圖，可以看到，兩個模型都是用OverallQual這個特征來分類，區(qū)別在于"<= 6.5"和"<=6"。很奇怪，明明OverallQual都是整數(shù)，6.5是從何而來的？我猜測"6.5 = 6 + 誤差修正"，它這樣做的目的是希望讓分類更精確。

%timeit RandomForest(trn_x, trn_y, 1)
10 loops, best of 3: 133 ms per loop

%timeit m.fit(trn_x, trn_y)
1000 loops, best of 3: 1.32 ms per loop

%timeit是用來測試程序運行所花費的時間的，相比%time，它取多次運行的均值，因此更精確?？梢钥吹?，現(xiàn)在這個決策樹運行一次需要133ms，比對標模型慢100倍。然后我用%prun來查看每個模塊具體的時間花銷，不出所料，問題在find_better_split()。

x <= x[i]以及x > x[i]的底層代碼也是用for循環(huán)實現(xiàn)的，那么find_better_split()的時間復雜度就是。顯然，如果x是已經(jīng)升序排序好的，那x <= x[i]就可以少做很多比較，而快速排序的時間復雜度是，遠小于。

除此之外，np.std()也比我想象中更花時間（14ms），因此，可以換一種計算方法來提升速度。Figure 2是標準的標準差計算公式，F(xiàn)igure 3是我們要采用的新公式。新公式不需要計算每個觀察值和均值的差，計算量更少，只需要求和的均值（E[X] == 均值）即可。

Figure 2

Figure 3

def calc_std(n, s, s2): return np.sqrt(s2/n - (s/n)**2)

def find_better_split(self, var_idx):
  x, y = self.x.iloc[self.idxs, var_idx].values, self.y[self.idxs]
  idxs = np.argsort(x)
  x_sorted, y_sorted = x[idxs], y[idxs]
  r_cnt, r_s, r_s2 = len(x_sorted), y_sorted.sum(), (y_sorted**2).sum()
  l_cnt, l_s, l_s2 = 0., 0., 0.

  for i in range(self.n - self.min_leaf):
    l_cnt += 1; r_cnt -= 1
    l_s += y_sorted[i]; r_s -= y_sorted[i]
    l_s2 += y_sorted[i]**2; r_s2 -= y_sorted[i]**2
    if x_sorted[i] == x_sorted[i + 1]: continue

    l_std = calc_std(l_cnt, l_s, l_s2)
    r_std = calc_std(r_cnt, r_s, r_s2)
    score = l_cnt * l_std + r_cnt * r_std
    if score < self.score:
      self.score, self.var_idx, self.split_val = score, var_idx, x_sorted[i]

可以看到，排序后的x，不會重復用同一類的樣本做參考點重復決策，少做了很多無用功，因此運行一次程序只需要8ms。

rf = RandomForest(trn_x, trn_y, 1)
print(rf.trees[0])
%timeit RandomForest(trn_x, trn_y, 1)

samples: 1000, value: 12.030368441695023, score: 288.6899540938478, split variable: OverallQual, split value: 6
100 loops, best of 3: 8.63 ms per loop

多次決策

def find_split(self):
  for c in range(self.c):
    self.find_better_split(c)
  if not self.is_leaf:
#     pdb.set_trace()
    x = self.x.iloc[:, self.var_idx].values
    l_idxs = np.nonzero(x > self.split_val)[0]
    r_idxs = np.nonzero(x <= self.split_val)[0]
    self.lb = DecisionTree(self.x.iloc[l_idxs], self.y[l_idxs], min_leaf=5)
    self.rb = DecisionTree(self.x.iloc[r_idxs], self.y[r_idxs], min_leaf=5)

讀到這里，如果你有編程經(jīng)驗，那會很容易從決策樹聯(lián)想到二叉樹，經(jīng)典的二叉樹是通過遞歸實現(xiàn)的。這里也是借鑒了遞歸的思想，在決策樹中生成決策樹，直到到達葉子節(jié)點為止。

t = RandomForest(trn_x, trn_y, 1).trees[0]
t.rb, t.lb, t.lb.lb, t.lb.rb

(samples: 617, value: 11.809290068879918, score: 151.00603314209667,
split variable: GrLivArea, split value: 1376,
 samples: 383, value: 12.386518196334503, score: 90.79697982952005,
split variable: OverallQual, split value: 7,
 samples: 162, value: 12.584355513747466, score: 36.72007562937957,
split variable: GrLivArea, split value: 1915,
 samples: 221, value: 12.241497176330432, score: 40.144323458391156,
split variable: GrLivArea, split value: 1935)

從兩個模型的測試結(jié)果來看，兩個模型的決策樹分類結(jié)果都是相同的，只是在分類值上略有不同。要驗證模型的實際預測效果還是需要通過 score。

R^2 Score

def predict(self, x):
  return [self.predict_row(xi[1]) for xi in x.iterrows()]

def predict_row(self, xi):
  if self.is_leaf: return self.val
  t = self.rb if xi[self.var_idx] <= self.split_val else self.lb
  return t.predict_row(xi)

DecisionTree.predict = predict
DecisionTree.predict_row = predict_row
metrics.r2_score(t.predict(val_x[:5]), val_y[:5]), metrics.r2_score(m.predict(val_x[:5]), val_y[:5])

(0.9243634901216589, 0.8916108333330888)

是衡量模型解釋數(shù)據(jù)變化的能力，也就是x變了y會怎么變，它的取值范圍是[]，值越大表示模型預測能力越強。我們的模型比對標模型表現(xiàn)得還要好！

Final model

萬事俱備，現(xiàn)在是時候用完整的數(shù)據(jù)來檢驗完整的模型，在此之前，需要為隨機森林添加predict()。

def predict(self, x):
  return np.mean([t.predict(x) for t in self.trees], 0)

RandomForest.predict = predict

另外，還記得前文提到的兩個隨機性嗎？采樣隨機性已經(jīng)實現(xiàn)了，就剩決策隨機性了。所謂的放棄“好決策”擁抱“壞決策”，實際上就是隨機地丟棄部分待決策的特征。max_features=0.5的意思就是隨機丟棄50%待決策特征。

class DecisionTree():
......
  def find_split(self):
#     pdb.set_trace()
    idxs = np.random.permutation(self.c)[:int(self.c * self.max_features)]
    for c in idxs:
      self.find_better_split(c)
......

最后，我們用同一個驗證集（val_x）分別驗證兩個模型。可以看到，我們這個不到80行的模型表現(xiàn)得比sklearn的模型還要好！但前者的運行時間卻是后者的1400多倍！

%time rf = RandomForest(trn_x, trn_y, 20, min_leaf=2, max_features=0.5, sample_sz=800)
metrics.r2_score(rf.predict(val_x), val_y)

CPU times: user 3min 21s, sys: 388 ms, total: 3min 22s
Wall time: 3min 22s
0.8489631454261701

m = RandomForestRegressor(n_estimators=20, n_jobs=1, min_samples_leaf=2, max_features=0.5)
%time m.fit(trn_x, trn_y)
metrics.r2_score(m.predict(val_x), val_y)

CPU times: user 138 ms, sys: 1.99 ms, total: 140 ms
Wall time: 140 ms
0.8454547075884457

Speed Up & Cython

我們自己編寫的模型之所以要比sklearn模型要慢得多，主要的原因在于sklearn是用Cython寫的。Python和Cython雖然語法上相差不大，但它們生成的程序有本質(zhì)區(qū)別。Python程序的運行依托于Python解釋器，而Cython像C語言一樣會編譯成機器碼，因此Cython程序的運行速度遠超Python程序。

通過運行一段簡單代碼，你就可以看到它們之間在執(zhí)行速度上質(zhì)的區(qū)別：

def foo1(n):
  i = 0
  p = 2
  while i < n:
    i += 1
    p = p * i
%timeit foo1(100)

100000 loops, best of 3: 12.4 μs per loop

%load_ext Cython
%%cython
def foo2(int n):
  cdef int i = 0
  cdef int p = 2
  while i < n:
    i += 1
    p = p * i
%timeit foo2(100)

10000000 loops, best of 3: 57.9 ns per loop

可以看到，對于這個簡單到爆的函數(shù)，python運行一次需要12皮秒，而編譯成Cython程序運行一次只需不到60納秒，兩者相差200倍！

因此，如果你有興趣，歡迎你用Cython重寫這段不到80行的程序，我相信它不會比sklearn模型慢太多。

END

本文通過從零開始編寫隨機森林的方式，解析隨機森林的核心原理：決策樹集成和采樣/決策隨機化，并深入分析了決策的工作原理和幾個核心參數(shù)：min_leaf、max_features，最后還簡單介紹了通過Cython可以從根本上改善整個模型的運行速度。