傳送門

《統計學習方法》讀書筆記——機器學習常用評價指標
《統計學習方法》讀書筆記——感知機（原理+代碼實現）
《統計學習方法》讀書筆記——K近鄰法（原理+代碼實現）
《統計學習方法》讀書筆記——樸素貝葉斯法（公式推導+代碼實現）

寫在前面

樸素貝葉斯法與貝葉斯估計是不同的概念。

?

損失函數與風險函數

損失函數用于度量一次預測的好壞；

風險函數用于度量平均意義下模型的好壞。

全概率公式與逆概率公式

設$A_1,A_2,...,A_n$A1?,A2?,...,An?為一組完備事件組，則對任一事件$B$B，有如下全概率公式：
$P(B) =\sum_{i=1}^nP(A_i)P(B|A_i)$P(B)=i=1∑n?P(Ai?)P(B∣Ai?)
若$P(B)>0$P(B)>0，則有如下貝葉斯公式，或稱逆概率公式：
$P(A_i|B)=\frac{P(A_iB)}{P(B)}=\frac{P(A_i)P(B|A_i)}{\sum\limits_{j=1}^nP(A_j)P(B|A_j)}$P(Ai?∣B)=P(B)P(Ai?B)?=j=1∑n?P(Aj?)P(B∣Aj?)P(Ai?)P(B∣Ai?)?

先驗概率與后驗概率

常稱$P(A_i)$P(Ai?)為事件$A_i$Ai?發生的先驗概率，而稱$P(A_i|B)$P(Ai?∣B)為事件$A_i$Ai?發生的后驗概率。(概率論教程P25)

樸素貝葉斯法

“如果你不知道怎樣踢球，就往球門方向踢 ” ——施拉普納

要搞明白樸素貝葉斯法的原理，要先知道“球門方向”在哪里。
設輸入空間$\mathcal{X}\subseteq{\bm{R}^n}$X?Rn為$n$n維向量的集合，輸入空間為類標記集合$\mathcal{Y} = \{c_1,c_2,...,c_K\}$Y={c1?,c2?,...,cK?}。輸入為特征向量$x\in\mathcal{X}$x∈X，輸入為類別$y\in\mathcal{Y}$y∈Y。$X$X是定義在輸入空間上的隨機變量，$Y$Y是定義在輸出空間$\mathcal{Y}$Y上的隨機變量。對于一個輸入的數據$x$x，要預測$x$x所屬的類別，相當于求以下概率的值：
$P(Y=c_k|X=x)$P(Y=ck?∣X=x)
根據貝葉斯公式（逆概率公式），可得：
$P(Y=c_k|X=x) = \frac{P(X=x|Y=c_k)P(Y=c_k)}{\sum\limits_k P(X=x|Y=c_k)P(Y=c_k)}$P(Y=ck?∣X=x)=k∑?P(X=x∣Y=ck?)P(Y=ck?)P(X=x∣Y=ck?)P(Y=ck?)?

插入部分
上述的貝葉斯公式中，$P(Y=c_k)$P(Y=ck?)容易求得，統計數據集中各個類別樣本的數量就可以計算出來。
因為$x$x是向量，即$x = (x^{(1)},x^{(2)},...,x^{(n)})$x=(x(1),x(2),...,x(n))，因此將$P(X=x|Y=c_k)$P(X=x∣Y=ck?)部分展開來就是$P(X^{(1)}=x^{(1)},X^{(2)}=x^{(2)},...,X^{(n)}=x^{(n)} |Y=c_k)$P(X(1)=x(1),X(2)=x(2),...,X(n)=x(n)∣Y=ck?)。
這樣看來，假設對mnist數據集來說，每一張圖片有$K=10$K=10種可能的預測結果（即$c_k,k=1,2,...,10$ck?,k=1,2,...,10），每一張圖片由28*28=784個像素點組成(即$x^{(j)}$x(j),$j=1,2,...,784$j=1,2,...,784)，每一個像素點的取值范圍為$0<c^{(j)}<255$0<c(j)<255，即$x^{(j)}$x(j)的可值為$S_j=255$Sj?=255個。那么要計算的數據量為$K\prod\limits_{j=1}^n=10\times255^{784}$Kj=1∏n?=10×255784，無法計算。
樸素貝葉斯法對條件概率分布作了條件獨立性的假設，樸素貝葉斯法也由此得名，即
$\begin{aligned} P(X=x|Y=c_k) &= P(X^{(1)}=x^{(1)},X^{(2)}=x^{(2)},...,X^{(n)}=x^{(n)} |Y=c_k)\\ &=P(X^{(1)}=x^{(1)}|Y=c_k)P(X^{(2)}=x^{(2)}|Y=c_k)...P(X^{(n)}=x^{(n)} |Y=c_k)\\ &=\prod\limits_{j=1}^nP(X^{(j)}=x^{(j)}|Y=c_k) \end{aligned}$P(X=x∣Y=ck?)?=P(X(1)=x(1),X(2)=x(2),...,X(n)=x(n)∣Y=ck?)=P(X(1)=x(1)∣Y=ck?)P(X(2)=x(2)∣Y=ck?)...P(X(n)=x(n)∣Y=ck?)=j=1∏n?P(X(j)=x(j)∣Y=ck?)?
這一假設使樸素貝葉斯法變得簡單，但會犧牲一定的分類準確率。
（因為向量的特征之間大概率是不獨立地，如果我們獨立了，會無法避免地拋棄一些前后連貫的信息，比方說我說“三人成？”，后面大概率就是個”虎“，這個虎明顯依賴于前面的三個字。）

那么上面的貝葉斯公式就可以轉換為以下形式：
$\begin{aligned} P(Y=c_k|X=x) &= \frac{P(X=x|Y=c_k)P(Y=c_k)}{\sum\limits_k P(X=x|Y=c_k)P(Y=c_k)} \\ &= \frac{P(Y=c_k)\prod\limits_{j=1}^nP(X^{(j)}=x^{(j)}|Y=c_k))}{\sum\limits_k P(Y=c_k)\prod\limits_{j=1}^nP(X^{(j)}=x^{(j)}|Y=c_k)} \end{aligned}$P(Y=ck?∣X=x)?=k∑?P(X=x∣Y=ck?)P(Y=ck?)P(X=x∣Y=ck?)P(Y=ck?)?=k∑?P(Y=ck?)j=1∏n?P(X(j)=x(j)∣Y=ck?)P(Y=ck?)j=1∏n?P(X(j)=x(j)∣Y=ck?))??
這是樸素貝葉斯法分類的基本公式。于是，樸素貝葉斯分類器可表示為
$y= \argmax_{c_k}\frac{P(Y=c_k)\prod\limits_{j=1}^nP(X^{(j)}=x^{(j)}|Y=c_k))}{\sum\limits_k P(Y=c_k)\prod\limits_{j=1}^nP(X^{(j)}=x^{(j)}|Y=c_k)}$y=ck?argmax?k∑?P(Y=ck?)j=1∏n?P(X(j)=x(j)∣Y=ck?)P(Y=ck?)j=1∏n?P(X(j)=x(j)∣Y=ck?))?
其中分母對所有$c_k$ck?都是相同的，所以
$y= \argmax_{c_k}P(Y=c_k)\prod\limits_{j=1}^nP(X^{(j)}=x^{(j)}|Y=c_k))$y=ck?argmax?P(Y=ck?)j=1∏n?P(X(j)=x(j)∣Y=ck?))
該公式即為樸素貝葉斯分類器。
至于式子里面的兩項具體怎么求，我們首先看第一項。
$P\left(Y=c_{k}\right)=\frac{\sum_{i=1}^{N} I\left(y_{i}=c_{k}\right)}{N}, \quad k=1,2, \cdots, K$P(Y=ck?)=N∑i=1N?I(yi?=ck?)?,k=1,2,?,K
N為訓練樣本的數目，假設我們手里現在有100個樣本，那N就是100。
分子中$I(...)$I(...)是指示函數，括號內條件為真時指示函數為1，反之為0。分子的意思求在這100個樣本里有多少是屬于$c_k$ck?分類的。
再看第二項
$P(X^{(j)}=a_{j l} |Y=c_{k})=\frac{\sum_{i=1}^{N} I(x_{i}^{(j)}=a_{j l},y_{i}=c_{k})}{\sum_{i=1}^{N} I(y_{i}=c_{k})} \\$P(X(j)=ajl?∣Y=ck?)=∑i=1N?I(yi?=ck?)∑i=1N?I(xi(j)?=ajl?,yi?=ck?)?
其中$j=1,2, ..., n ; \quad l=1,2, ..., S_{j};\quad k=1,2,..., K$j=1,2,...,n;l=1,2,...,Sj?;k=1,2,...,K
構造原理和上面第一項的一樣，也是通過指示函數來計數得出概率。那么兩項都能求了，給出樸素貝葉斯算法。

但是，步驟（2）中，那么多概率連乘，如果其中有一個概率為0怎么辦？那整個式子直接就是0了，這樣不對。所以我們連乘中的每一項都得想辦法讓它保證不是0。

代碼實現

# coding=utf-8
import numpy as np
import time
def loadData(fileName):
    print('start to read data：' + fileName)
    dataArr = []
    labelArr = []
    fr = open(fileName, 'r')
    # 將文件按行讀取
    for line in fr.readlines():
        # 對每一行數據按切割符','進行切割，返回字段列表
        curLine = line.strip().split(',')
        # 將mnist圖像二值化
        dataArr.append([int(int(num) > 128) for num in curLine[1:]])
        labelArr.append(int(curLine[0]))
    # 返回data和label
    return dataArr, labelArr


def NaiveBayes(Py, Px_y, x):
    featrueNum = 784
    classNum = 10
    P = [0] * classNum
    for i in range(classNum):
        sum = 0
        for j in range(featrueNum):
            sum += Px_y[i][j][x[j]]
        P[i] = sum + Py[i]
    return P.index(max(P))


def model_test(Py, Px_y, testDataArr, testLabelArr):
    errorCnt = 0
    for i in range(len(testDataArr)):
        predict = NaiveBayes(Py, Px_y, testDataArr[i])
        if predict != testLabelArr[i]:
            errorCnt += 1
    return 1 - (errorCnt / len(testDataArr))


def getAllProbability(trainDataArr, trainLabelArr):
    featureNum = 784
    classNum = 10
    lambbaVal = 1

    #初始化先驗概率分布存放數組，后續計算得到的P(Y = 0)放在Py[0]中，以此類推
    #數據長度為10行1列
    Py = np.zeros((classNum, 1))
    #對每個類別進行一次循環，分別計算它們的先驗概率分布
    #計算公式為書中"4.2節 樸素貝葉斯法的參數估計 公式4.8"
    for i in range(classNum):
        # 計算先驗概率
        Py[i] = ((np.sum(np.mat(trainLabelArr) == i)) + lambbaVal) / (len(trainLabelArr) + classNum * lambbaVal)
    # 防止數據過小向下溢出使用log()
    # 在似然函數中通常會使用log的方式進行處理
    Py = np.log(Py)

    # 計算條件概率 Px_y=P（X=x|Y = y）
    # 計算分子部分
    # 2表示一個像素點可取值的個數（因為對數據做了二值化處理）
    Px_y = np.zeros((classNum, featureNum, 2))
    for i in range(len(trainLabelArr)):
        #獲取當前循環所使用的標記
        label = trainLabelArr[i]
        #獲取當前要處理的樣本
        x = trainDataArr[i]
        #對該樣本的每一維特診進行遍歷
        for j in range(featureNum):
            #在矩陣中對應位置加1
            #這里還沒有計算條件概率，先把所有數累加，全加完以后，在后續步驟中再求對應的條件概率
            Px_y[label][j][x[j]] += 1


    # 計算分母部分
    for label in range(classNum):
        for j in range(featureNum):
            #獲取y=label，第j個特診為0的個數
            Px_y0 = Px_y[label][j][0]
            #獲取y=label，第j個特診為1的個數
            Px_y1 = Px_y[label][j][1]
            #對式4.10的分子和分母進行相除，再除之前依據貝葉斯估計，分母需要加上2（為每個特征可取值個數）
            #分別計算對于y= label，x第j個特征為0和1的條件概率分布
            Px_y[label][j][0] = np.log((Px_y0 + 1) / (Px_y0 + Px_y1 + 2))
            Px_y[label][j][1] = np.log((Px_y1 + 1) / (Px_y0 + Px_y1 + 2))

    #返回先驗概率分布和條件概率分布
    return Py, Px_y


if __name__ == "__main__":
    start = time.time()
    # 獲取訓練集、測試集
    trainData, trainLabel = loadData('./mnist/mnist_train.csv')
    testData, testLabel = loadData('./mnist/mnist_test.csv')

    #開始訓練，學習先驗概率分布和條件概率分布
    print('start to train')
    Py, Px_y = getAllProbability(trainData, trainLabel)

    #使用習得的先驗概率分布和條件概率分布對測試集進行測試
    print('start to test')
    accuracy = model_test(Py, Px_y, testData, testLabel)

    print('the accuracy is:', accuracy)
    print('time span:', time.time() -start)

輸出結果

start to read data：./mnist/mnist_train.csv
start to read data：./mnist/mnist_test.csv
start to train
start to test
the accuracy is: 0.8432999999999999
time span: 129.84226727485657

參考

原理：《統計學習方法》
代碼： https://github.com/Dod-o/Statistical-Learning-Method_Code