統計學習1：朴素貝葉斯模型(Numpy實現)

模型

生成模型介紹

我們定義樣本空間為\(\mathcal{X} \subseteq \mathbb{R}^n\)，輸出空間為\(\mathcal{Y} = \{c_1, c_2, ..., c_K\}\)。\(\textbf{X}\)為輸入空間上的隨機向量，其取值為\(\textbf{x}\)，滿足\(\textbf{x} \in \mathcal{X}\)；\(Y\)為輸出空間上的隨機變量，設其取值為\(y\)，滿足\(y \in \mathcal{Y}\)。我們將容量為\(m\)的訓練樣本表示為:

\[\begin{aligned} D = \{\{\textbf{x}^{(1)}, y^{(1)}\}, \{\textbf{x}^{(2)}, y^{(2)}\},..., \{\textbf{x}^{(m)}, y^{(m)}\}\} \end{aligned}\tag{1} \]

\[ \]

我們遵循機器學習的一個基本假設，即訓練樣本是從一個未知的總體分布\(P(\textbf{X} = \textbf{x}, Y=y)\)中采樣產生，且訓練樣本獨立同分布。

我們采取概率模型的視角，即將分類模型表示為條件概率分布\(P(Y=y|\textbf{X}=\textbf{x})\)。而依據分布\(P(Y=y|\textbf{X}=\textbf{x})\)的求解可將模型分為判別模型和生成模型。判別模型直接對條件概率分布\(P(Y=y|\textbf{X}=\textbf{x})\)進行參數估計(估計方法可采用極大似然估計或貝葉斯估計)；而生成模型則利用條件概率公式\(P(Y=y|\textbf{X}=\textbf{x}) = \frac{P(\textbf{X}=\textbf{x}, Y=y)}{P(\textbf{X}=\textbf{x})}\)來計算分布。分子\(P(\textbf{X}=\textbf{x}, Y=y)\)是一個聯合概率分布，能夠還原出聯合概率分布\(P(\textbf{X}=\textbf{x}, Y=y)\)是生成模型的一大特性。

朴素貝葉斯模型推導

我們對分子繼續運用條件概率公式，進一步得到

\[\begin{aligned} P(Y=y|\textbf{X}=\textbf{x}) = \frac{P(\textbf{X}=\textbf{x}|Y=y)P(Y=y)}{P(\textbf{X}=\textbf{x})} \end{aligned} \tag{2} \]

這個公式即大名鼎鼎的貝葉斯公式。 這里我們采用貝葉斯學派的視角，將\(P(Y=y)\)稱為先驗概率分布，表示在數據觀測之前對\(Y\)的信念；\(P(Y = y|\textbf{X}=\textbf{x})\)稱為后驗概率分布，表示經過觀測數據\(\textbf{X}\)(也稱“證據”)校正后對\(Y\)的信念。注意不要和和貝葉斯估計中參數\(\theta\)的先驗和后驗分布搞混了，貝葉斯估計也應用了貝葉斯公式，但先驗概率分布和后驗概率分布的實際含義與這里完全不同。

我們再將分母運用全概率公式展開，我們得到

\[\begin{aligned} P(Y=y|\textbf{X}=\textbf{x})= \frac{P(\textbf{X}=\textbf{x}|Y=y)P(Y=y)}{\sum_{y\in \mathcal{Y}}P(\textbf{X}=\textbf{x}|Y=y)P(Y=y)} \end{aligned}\tag{3} \]

這意味着我們只需要學習概率分布\(P(Y=y)\)和\(P(\textbf{X}=\textbf{x}|Y=y)\)，而無需關心\(P(\textbf{X}=\textbf{x})\)。

將隨機向量\(\textbf{x}\)沿着其特征維度展開，我們繼續得到

\[\begin{aligned} P(\textbf{X} = \textbf{x} | Y = y) = P(X_1 = x_1, ..., X_n = x_n | Y=y), \quad n \text{為特征維度} \end{aligned}\tag{4} \]

這里我們為了簡單起見，假設樣本屬性是離散的，第\(j\)個屬性\(x_j\)的屬性集為\(A_j=\{a_{j1}, a_{j2},...,a_{jl},..., a_{j{N_j}}\}\)，滿足\(x_j \in A_j\)。可以看出，條件概率分布\(P(\textbf{X}=\textbf{x}|Y=y)\)的參數總量是指數級的(\(x_j\)的屬性集\(A_j\)大小為\(N_j\)，\(j=1, 2, ..., n\)，\(Y\)可取值有\(K\)個，那么參數個數為\(K \prod_{j=1}^{n}N_j\))，不能對其直接進行參數估計。

因此，我們決定對原本擁有指數級參數數量的分布進行拆分。這里，朴素貝葉斯法做出了條件獨立性假設：樣本特征在類確定的條件下條件獨立(這也是“朴素”(Naive)一詞的得名)。這樣我們就能將原本擁有龐大參數的概率分布進行拆分：

\[\begin{aligned} P(\textbf{X} = \textbf{x} | Y=y) = P(X_1 = x_1, ..., X_n = x_n | Y=y) = \prod_{j=1}^nP(X_j = x_j | Y = y) \end{aligned}\tag{5} \]

這樣，我們就可以對\(P(\textbf{X} | Y=y)\)分布進行高效的參數估計。之后，我們對於輸入樣本\(\textbf{x}\)，計算概率分布\(P(Y=y|\textbf{X}=\textbf{x})\)：

\[\begin{aligned} P(Y=y|\textbf{X}=\textbf{x})= \frac{P(\textbf{X}=\textbf{x}|Y=y)P(Y=y)}{\sum_{ y \in \mathcal{Y}}P(\textbf{X}=\textbf{x}|Y=y)P(Y=y)} = \frac{\prod_{j=1}^nP(X_j = x_j | Y=y)P(Y=y)}{\sum_{y \in \mathcal{Y}}[ \prod_{j=1}^nP(X_j = x_j | Y=y)P(Y=y) ]} \end{aligned}\tag{6} \]

我們采取后驗概率最大化原則(即最終的輸出分類取使條件概率最大的那個)，設\(f(\textbf{x})\)為分類決策函數，即

\[\begin{aligned} y = f(\textbf{x}) = \underset{y}{\arg\max} P(Y = y|\textbf{X}=\textbf{x})=\underset{y}{\arg\max}\frac{\prod_{j=1}^nP(X_j = x_j | Y=y)P(Y=y)}{\sum_{y \in \mathcal{Y}}[ \prod_{j=1}^nP(X_j = x_j | Y=y)P(Y=y) ]} \end{aligned}\tag{7} \]

我們發現，不管\(y\)取何值，式\((7)\)中分母總是恆定的，因此我們可以將式\((7)\)化簡為

\[\begin{aligned} y = f(\textbf{x}) = \underset{y}{\arg\max} P(Y=y|\textbf{X}=\textbf{x}) = \underset{y}{\arg\max}\prod_{j=1}^nP(X_j = x_j | Y=y)P(Y=y) \end{aligned}\tag{8} \]

這就是朴素貝葉斯模型分類決策函數的最終表達式。

參數估計

極大似然估計

如式\((8)\)中所述，我們需要對先驗概率分布\(P(Y=y)\)和條件概率分布\(P(X_j = x_j|Y=y)\)進行參數估計。根據極大似然估計(具體的推導過程可以參見李航《統計學習方法》中的習題解答)，我們可以運用訓練集\(D\)將先驗概率分布\(P(Y=y)\)估計為

\[\begin{aligned} P(Y=y) = \frac{\sum_{i=1}^m(y^{(i)} = y)}{m}, \quad m \text{為}D \text{中樣本個數} \end{aligned} \tag{9} \]

同樣，條件概率分布\(P(X_j = x_j|Y=y)\)的估計為

\[\begin{aligned} P(X_j = x_j | Y=y) = \frac{P(X_j = x_j, Y = y)}{P(Y = y)} = \frac{\sum_{i=1}^{m}I(x_j^{(i)} =x_j, y^{(i)}=y)}{\sum_{i=1}^{m}I(y^{(i)}=y)} , \quad m \text{為}D \text{中樣本個數} \end{aligned} \tag{10} \]

貝葉斯估計(平滑修正)

觀察式\((10)\)可知，如果訓練集中屬性值\(x_j\)和類\(y\)沒有同時出現過，即\(P(X_j=x_j, Y=y)=0\)，那么\(P(X_j = x_j | Y=y)=0\)會直接導致連乘式。這就意味着不管其他屬性如何，哪怕其他屬性明顯符合要求，樣本\(\prod_{j=1}^nP(X_j = x_j | Y=y)=0\) ，\(\textbf{x}\)屬於類\(y\)的概率都會被判為0，這明顯不太合理。

因此，為了避免其他屬性攜帶的信息被訓練集中未出現的屬性值“抹去”，我們采用貝葉斯估計，等價於在估計概率值時通常進行“平滑”(smoothing)(具體的推導過程可以參見李航《統計學習方法》中的習題解答)。即令式\((10)\)修正為

\[\begin{aligned} P_{\lambda}(X_j = x_j|Y=y) = \frac{\sum_{i=1}^{m}I(x_j^{(i)} =x_j, y^{(i)}=y) + \lambda}{\sum_{i=1}^{m}I(y^{(i)}=y)+N_j \lambda}, \quad \lambda > 0, \quad N_j\text{為屬性}x_j\text{可能的取值數} \end{aligned} \tag{11} \]

我們常取\(\lambda=1\)，這時稱為拉普拉斯平滑(Laplacian smoothing)。

類似地，式\((9)\)中先驗概率被修正為：

\[\begin{aligned} P_\lambda(Y=y) = \frac{\sum_{i=1}^m(y^{(i)} = y)+\lambda}{m+K\lambda}, \quad \lambda > 0, \quad K\text{為標簽}y\text{可能的取值數} \end{aligned} \tag{12} \]

可以看出，拉普拉斯平滑解決了訓練集樣本數不足導致的概率值為 0 的問題。拉普拉斯修正實際上假設了屬性值與類別均勻分布，這是在參數估計的過程中額外引入的關於數據的先驗 (prior)。當樣本容量趨近於無窮時，我們發現修正過程所引入的先驗的影響也趨近於 0，使得計算的概率值趨近於實際的概率值。

算法

在實際的應用中，朴素貝葉斯模型有兩種訓練方式。

若使用的場景對模型的預測速度要求較高，在給定訓練集\(D\)的情況下，我們將概率分布\(P_\lambda(Y=y)\)和概率分布\(P_{\lambda}(X_j = x_j|Y=y)\)所有可能的取值(\(y\in \mathcal{Y}\)，\(x_j \in A_j\)，\(A_j\)為第\(j\)個樣本屬性的取值集合)都計算出來存好，然后在測試樣本\(\textbf{x}^{*}\)來了之后，通過“查表”的方式將對應的概率值檢索出來，然后再對其類別進行判別。這樣，我們計算概率分布\(P_\lambda(Y=y)\)和概率分布\(P_{\lambda}(X_j = x_j|Y=y)\)所有可能取值的過程即對朴素貝葉斯模型進行顯式訓練的過程。
朴素貝葉斯的訓練和測試算法如下：

如果我們不斷有新的訓練數據產生，可以采用“懶惰學習”(lazy learning)的方法，先不進行任何訓練，測試樣本來了之后再依照測試樣本的屬性\(x_j^{*}\)和當前數據集的狀況來計算單點概率，這樣可以避免對所有可能的屬性都計算單點概率。若訓練數據不斷增加，則可在現有計算結果的基礎上，僅僅對新增樣本的屬性值所涉及的單點概率進行計數修正，這樣可以實現“增量學習”。

代碼實現

使用Python語言實現朴素貝葉斯算法如下（這里我們采用Iris數據集進行測試）。

from numpy.lib.index_tricks import c_
from sklearn.datasets import load_iris
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score
import numpy as np
from copy import deepcopy
class NaiveBayes():
    def __init__(self, A, C, lambda_v=1):
        #常常取lambda_v=1，此時稱為拉普拉斯平滑
        self.K = len(C)
        self.A = deepcopy(A)
        self.C = deepcopy(C)
        self.y_cnt = {} #記錄sum I(y_i = y)
        self.y_prob = {}  # 記錄平滑后的P_lambda(Y = y)
        self.lambda_v = lambda_v # 平滑參數
        self.x_condition_y_prob = {} 
        #記錄P_lambda(X_j = x_j | Y = y),這是個字典-列表-字典嵌套
        # 初始化y_cnt和y_prob
        for c_k in self.C:
            self.y_cnt[c_k], self.y_prob[c_k] = 0, 0
        # 初始化attr_cnt和x_condition_y_prob
        for c_k in self.C:
            self.x_condition_y_prob[c_k] = [dict() for j in range(len(self.A))]
            for j in range(len(self.A)):
                for a_j_l in self.A[j]:
                    self.x_condition_y_prob[c_k][j][a_j_l] = 0

    def fit(self, X_train, y_train):
        try:
            assert(X_train.shape[0] == y_train.shape[0])
        except:
            AssertionError("input dimension 0 should be the same!")

        # 記錄 sum I(y_i = y)
        for y in y_train:
            try:
                self.y_cnt[y] += 1 
                #如果訓練集中沒出現過，此處y_cnt[y]將一直為0
            except:
                raise ValueError("invalid label!")   
        # 計算P_lambda(Y = c_k)
        # 並對所有概率進行平滑處理
        for c_k in self.C:
            self.y_prob[c_k] = self.y_cnt[c_k] 
            self.y_prob[c_k] += self.lambda_v
            self.y_prob[c_k] /= (X_train.shape[0] + self.K * self.lambda_v)
    
        # 記錄sum I(x_j(i) = x_j, y(i) = y)
        # 遍歷每一個訓練樣本
        for x, y in zip(X_train, y_train):
            # 遍歷訓練樣本中的每一個屬性
            for j, attr in enumerate(x):
                try:
                    self.x_condition_y_prob[y][j][attr] += 1 
                    #如果訓練集中沒出現過，此處self.x_condition_y_prob[y][attr] 將一直為0
                except:
                    raise ValueError("invalid attribution %.1f!" % attr)
        # 計算P_lambda(X_j = x_j | Y = c_k)
        # 並對所有樣本進行平滑處理
        for c_k in self.C:
            c_cnt = self.y_cnt[c_k] 
            # 遍歷所有屬性的屬性集合
            for j in range(len(self.A)):
                # 遍歷每一個屬性的取值集合
                for a_j_l in self.A[j]:
                    self.x_condition_y_prob[c_k][j][a_j_l] += self.lambda_v
                    self.x_condition_y_prob[c_k][j][a_j_l]/=(c_cnt + len(self.A[j])*self.lambda_v)             
            
    def pred(self, X_test):
        # 遍歷測試集中的每一個樣本
        y_pred = []
        for x in X_test:
            max_ck = self.C[0]
            max_prob = -1

            # 考察每一種可能的標簽
            for c_k in self.C:
                prob = 1.0
                # 考察該測試樣本的每一個屬性
                for j, attr in enumerate(x):
                    prob *= self.x_condition_y_prob[c_k][j][attr]

                prob *= self.y_prob[c_k]
                if prob == 0:
                    raise ValueError("prob underflow!")
                if prob > max_prob:
                    max_prob = prob
                    max_ck = c_k
            y_pred.append(max_ck)
        return np.array(y_pred)
                
# 獲得特征和類標記的結合
def get_A_and_C(X, y):
    # A為第j各屬性的可能取值集合
    A, C = [set() for j in range(X.shape[1])], set()
    # 遍歷每個樣本的特征向量
    for x in X:
        for j, attr in enumerate(x):
            A[j].add(attr)
    # 遍歷每個樣本的標簽
    for c_k in y:
        C.add(c_k)
    A = [ list(A_j) for A_j in A]
    return list(A), list(C)

if __name__ == "__main__":
    X, y = load_iris(return_X_y=True)
    A, C = get_A_and_C(X, y)
    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.5, random_state=0)
    clf = NaiveBayes(A, C)
    clf.fit(X_train, y_train)
    y_pred = clf.pred(X_test)
    acc_score = accuracy_score(y_test, y_pred)
    print("The accuracy is: %.1f" % acc_score)

最終的測試結果如下：

The accuracy is: 0.9

可以看出我們實現的算法在Iris數據集上取得了90%的精度，說明我們算法的效果不錯。

參考文獻

• [1] 李航. 統計學習方法(第2版)[M]. 清華大學出版社, 2019.
• [2] 周志華. 機器學習[M]. 清華大學出版社, 2016.
• [3] Calder K. Statistical inference[J]. New York: Holt, 1953.