程式碼+剖析|樸素貝葉斯原理剖析及實現
前言
《統計學習方法》一書在前幾天正式看完，由於這本書在一定程度上對於初學者是有一些難度的，趁著熱乎勁把自己走過的彎路都寫出來，後人也能走得更順暢一點。
以下是我的github地址，其中有《統計學習方法》書中所有涉及到的演算法的實現，也是在寫部落格的同時編寫的。在編寫宗旨上是希望所有人看完書以後，按照程式的思路再配合書中的公式，能知道怎麼講所學的都應用上。（有點自戀地講）程式中的註釋可能比大部分的部落格內容都多。希望大家能夠多多捧場，如果可以的話，為我的github打顆星，也能讓更多的人看到。
github：
GitHub|手寫實現李航《統計學習方法》書中全部演算法
正文
樸素貝葉斯的直觀理解
在網上曾經有一個有意思的機率討論，題目是這樣的（我相信大家會願意看一看這種有趣的問題）：
有三張彩票，一張有獎。你買了一張，老闆在自己剩餘的兩張中刮開了一張，沒中。這時候他要用剩下的一張和你換，你換不換？換和不換你中獎的機率一樣嗎？（你可以思考一下，然後看我下面的回答）
—————————————————————————-
從直覺上來講，彩票中獎的機率是1/3，你最先抽了一張，不管咋操作，中獎的機率應該都是1/3。這時候老闆排除掉了一張沒中獎的，剩下兩張必有一張中獎，所以機率是1/2。換和不換應該都一樣。你是這麼答的嗎？
這時候需要引申出貝葉斯了，貝葉斯在機率的計算中引入了先驗機率，即在知道了一些事情的情況下，機率會發生變化，按照貝葉斯的解答應該是怎樣的呢？
首先我們要引入機率的一些寫法（事實上如果你不太明白這些公式的話，建議適當補充一些機率的課程，在機器學習中非常有用）
補充知識：P（A|B）表示在B發生的情況下A發生的機率。P(A|B) = [P(B|A)*P(A)] / P(B)。
—————————————————————————-
我們假設A表示你手裡目前這張是中獎的，B表示老闆刮出來的那張是沒中獎的（雖然題目中說明了老闆刮出來的沒中，但還是需要將情況假設出來）
那麼由上面可得：
P（A|B） = ［P（B|A）*P（A）］ / P（B）
在式中P（B|A）表示在A發生的情況下B發生的機率，代入題目中也就是如果你手裡這張是中獎的，那老闆刮出來那張沒中獎的機率。很顯然這是必然的，因為只有一張是中獎的，所以P（B|A）=1，也就是必然的事件。P（A）表示你手裡目前這張是中獎的機率，在計算P（A）時因為沒有給出前提條件，所以P（A）=1/3，也就是三張彩票，你手裡中獎的機率是1/3。P（B）表示老闆刮出來那張是沒中獎的機率，題目中已經給出前提條件老闆刮出來的是沒中獎的，所以機率P（A）=1。
那麼將數代進式子裡面：P（A|B） = ［P（B|A）*P（A）］ / P（B） = ［1 * 1/3］ / 1 = 1/3，也就是說在老闆刮出來一張沒中獎的彩票前提下，你手裡這張中獎的機率是1/3，老闆手裡剩餘那張是2/3。為了中獎機率最大化，應該和老闆手裡的彩票交換。
是不是感覺不對勁，三張彩票我自己抽了一張，我中不中獎怎麼還和老闆刮不刮自己的彩票有關係了呢？P（A|B）這個演算法是錯的吧？當時論壇也爭論了很久，但十年前的人們還不太會使用貝葉斯來武裝自己，有一部分用直覺，也就是腦袋演算法，堅持就是1/2。另一部分去繞來繞去分析，企圖找到直覺沒考慮到的地方來證明是1/3和2/3。後來有位程式設計師寫了演算法來模擬，結果確實是1/3和2/3。有興趣的同學可以透過這個連結自行跑一下程式：一道有爭議的機率題
至於你們問我這個現象該怎麼解釋，可能需要很長的篇幅，但在本文中這不是重點展開的地方，這個例子和今天的樸素貝葉斯分類器有什麼關係？很抱歉，沒有關係，我只是想皮一下。（開個玩笑，我覺得這對於讀者來說，理解前提條件對於機率來說是一種什麼概念就可以了。）
好了，我們開始正式地講解樸素貝葉斯分類器是個啥玩意了。
機器學習演算法中我們通常在這樣的情況下去使用樸素貝葉斯分類器。
假設我們有一個手寫的資料集樣本，裡面有100條記錄，記錄0-10是10個人各自寫的0，記錄11-20是10個人各自寫的1，……。。以此類推一共100條記錄。那麼這時候外頭有個人進來寫了個數字X，怎麼識別出來它寫的幾呢？沒學習過機器學習的人可能也能提出這樣一種方法：我們只要把寫的那個數字和0-9進行匹配，那個匹配度最高就是哪個數啦。沒錯，樸素貝葉斯用的就是這樣樸素的思想（開玩笑，這裡的樸素可不是這個意思）。
樸素貝葉斯工作原理：假設Y表示是數字幾，寫的那個數叫X，那麼我們可以透過某種方法求P（Y = 0 | X），P（Y = 1 | X），……。，P（Y = 9| X）。其中P（Y = 1 | X）表示在給定手寫數字X的情況下，它是1的機率。這樣得出10個數字各自的機率，哪個高就最有可能是哪個。那麼咋求這個P（Y = 1 | X），彆著急，這就是後文重點討論的內容，在直觀理解中無須涉及。
總之，樸素貝葉斯分類器就是一個對所有可能性求機率的模型，最後輸出結果中哪種可能性高就輸出哪種。核心公式是P（Y | X），至於P（Y | X） = ？，這部分怎麼求？我們現在開始用數學的知識開始討論。（我儘可能簡化一些公式，但本章節相對於前幾章來說，公式有些稍多。讀者需要習慣，學會樸素貝葉斯分類器，我們一隻腳也就正式踏入了機器學習的大門）
樸素貝葉斯分類器的數學角度（配合《統計學習方法》食用更佳）
我們剛剛提到樸素貝葉斯最終需要計算得到的公式是P（Y=？|X），然後比較不同機率的大小，選擇機率最大的一類輸出。不妨將不同的類設為Ck，在手寫識別中C0就表示數字0，C1表示數字1，以此類推。那麼上式也就可以寫成了P（Y=Ck|X=x），即在給定樣本X為x的情況下，Y為Ck的機率。依據上文的貝葉斯公式，我們可以將該式轉換為下圖中的第一個等號右邊這樣（依據貝葉斯公式定義）。
在第一個等號右邊的式子中，分子的左邊P（X=x|Y=y）表示在給定數字幾的情況下，出現該樣本的機率，比如說隨便給了一個數字不知道是幾，然後被通知這個數是1，問你紙上會出現這個圖案的機率（也就是說告訴你這個數是1，那麼紙上出現的那團圖案長得像1的可能性就大，長得像8的可能性就小），這個在給定樣本進行訓練的情況下是可以統計出來的（具體怎麼統計會在後文講解）。P（Y=Ck）表示每個類別出現的機率，比如說10個數裡面有5個1，那P（Y=C1）=5/10，也可以透過給定訓練樣本進行統計得到。
式子唯一還不太明白的就是分母的P（X=x），給定一個樣本的機率？既然不清楚，那就把這個式子等價轉換一下，將分母轉換成與第二個等號右邊一樣。P（X=x|Y=Ck）P（Y=Ck）=P（X=x， Y=Ck），表示對於每一類Y，都求一次X=x的機率，那麼總的求和以後就變成了P（X=x）。（這部分不太明白的建議再補充一些機率的課程）。那麼最初的式子就變成了求第二個等號右邊的式子。
再往下轉換，就變成了下圖這樣。
和上圖的式子進行比對，其實就是把P（X=x|Y=Ck）這一項變成了連乘，至於為什麼能連乘，下圖有詳細說明：
為什麼可以把裡面的直接拆開來連乘？機率老師不是說過只有相互獨立才能直接拆嗎？是的，樸素貝葉斯分類器對條件機率分佈做出了條件獨立性的假設。為啥？因為這樣能算，就這麼簡單，如果條件都不獨立，後面咋整？（讀者：那你這不嚴謹啊）。emmm…。事實上是這樣，向量的特徵之間大機率是不獨立地，如果我們假設獨立了，會無法避免地拋棄一些前後連貫的資訊（比方說我說“好事成_”，後面大機率就是個”雙“，這個雙明顯依賴於前面的三個字）。在建立模型時如果這些都考慮進去，會讓模型變得很複雜。後來前人說那我們試試不管它們，強行獨立。誒發現效果還不錯誒，那就這麼用吧。這就是計算機科學家和數學家的分歧所在。
上圖中P（X=x|Y=Ck）轉換成能求的式子了以後，那麼就是比較Y為不同Ck的情況下哪個機率最大，那就表示屬於哪個類的可能性最大。所以前頭式子前頭加上一個argmax，表示求讓後式值最大的Ck。
然後由於下圖中圈出來這一項是在Y為不同Ck情況下的連乘，所以不管k為多少，所有Ck連乘結果肯定是一致的，在比較誰的值最大時，式子裡面的常數無法對結果的大小造成影響，可以去掉。
就變成了下面這樣：
這就是最終簡化版的樸素貝葉斯分類器。至於式子裡面的兩項具體怎麼求，我們首先看第一項。
N為訓練樣本的數目，假設我們手裡現在有100個樣本，那N就是100。分子中I是指示函式，括號內條件為真時指示函式為1，反之為0。分子啥意思呢？就是說對於Ck，在這100個樣本離有多少是Ck分類的，就比如Ck為C1時表示數字類別1，這時候就是看這100個樣本里有多少個數字1，處於總的100，就是類別1的機率，也就是P（Y=C1）。簡單不？相當簡單。我們再看第二項。
這個我覺得不用細說了，構造原理和上面第一項的一樣，也是透過指示函式來計數得出機率。那麼兩項都能求了，給出樸素貝葉斯演算法。
樸素貝葉斯演算法
可以看到在步驟中首先根據給定的訓練樣本求先驗機率和條件機率，也就是上文求的那兩個式子，然後訓練計數後給定一個樣本，計算在不同類別下該樣本出現的機率，求得最大值即可。
關於樸素貝葉斯和貝葉斯：
還記得原先求不同類別下的最大機率這一式子嗎？裡面有很多的連乘記得嗎？
這裡提出了一個問題，那麼多機率連乘，如果其中有一個機率為0怎麼辦？那整個式子直接就是0了，這樣不對。所以我們連乘中的每一項都得想辦法讓它保證不是0，哪怕它原先是0，（如果原先是0，表示在所有連乘項中它機率最小，那麼轉換完以後只要仍然保證它的值最小，對於結果的大小來說沒有影響。）。這裡就使用到了貝葉斯估計。
它透過分母和分子各加上一個數來保證始終不為0。為什麼要分母加Sj或者K倍λ，以及分子要加個1倍λ，讀者可以試驗一下，這樣處理以後，所有機率的總和仍然是1，所以為什麼這麼加，只不過這麼加完以後，機率仍然是機率，總和仍然為1，但所有項都大於0。
舉個例子：
關於《統計學習方法》的一些補充：
使用貝葉斯估計雖然保證了所有連乘項的機率都大於0，不會再出現某一項為0結果為0的情況。但若一個樣本資料時高維的，比如說100維（100其實並不高），連乘項都是0-1之間的，那100個0-1之間的數相乘，最後的數一定是非常非常小了，可能無限接近於0。對於程式而言過於接近0的數可能會造成下溢位，也就是精度不夠表達了。所以我們會給整個連乘項取對數，這樣哪怕所有連乘最後結果無限接近0，那取完log以後數也會變得很大（雖然是負的很大），計算機就可以表示了。同樣，多項連乘取對數，對數的連乘可以表示成對數的相加，在計算上也簡便了。所以在實際運用中，不光需要使用貝葉斯估計（保證機率不為0），同時也要取對數（保證連乘結果不下溢位）。
有人可能關心取完log以後結果會不會發生變化，答案是不會的。log在定義域內是遞增函式，log（x）中的x也是遞增函式（x在x軸的最大處，x值也就越大，有些拗口，但最終導致x越大，取完對數的log（x）同樣越大，可以仔細想想）。在單調性相同的情況下，連乘得到的結果大，log取完也同樣大，並不影響不同的連乘結果的大小的比較。好啦，上程式碼。
貼程式碼：（建議去本文最上方的github連結下載，有書中所有演算法的實現以及詳細註釋）
# coding=utf-8
# Author：Dodo
# Date：2018-11-17
# Email：lvtengchao@pku。edu。cn
# Blog：www。pkudodo。com
‘’‘
資料集：Mnist
訓練集數量：60000
測試集數量：10000
————————————————
執行結果：
正確率：84。3%
執行時長：103s
’‘’
import
numpy
as
np
import
time
def
loadData
（
fileName
）：
‘’‘
載入檔案
：param fileName：要載入的檔案路徑
：return： 資料集和標籤集
’‘’
#存放資料及標記
dataArr
=
［］；
labelArr
=
［］
#讀取檔案
fr
=
open
（
fileName
）
#遍歷檔案中的每一行
for
line
in
fr
。
readlines
（）：
#獲取當前行，並按“，”切割成欄位放入列表中
#strip：去掉每行字串首尾指定的字元（預設空格或換行符）
#split：按照指定的字元將字串切割成每個欄位，返回列表形式
curLine
=
line
。
strip
（）
。
split
（
‘，’
）
#將每行中除標記外的資料放入資料集中（curLine［0］為標記資訊）
#在放入的同時將原先字串形式的資料轉換為整型
#此外將資料進行了二值化處理，大於128的轉換成1，小於的轉換成0，方便後續計算
dataArr
。
append
（［
int
（
int
（
num
）
>
128
）
for
num
in
curLine
［
1
：］］）
#將標記資訊放入標記集中
#放入的同時將標記轉換為整型
labelArr
。
append
（
int
（
curLine
［
0
］））
#返回資料集和標記
return
dataArr
，
labelArr
def
NaiveBayes
（
Py
，
Px_y
，
x
）：
‘’‘
透過樸素貝葉斯進行機率估計
：param Py： 先驗機率分佈
：param Px_y： 條件機率分佈
：param x： 要估計的樣本x
：return： 返回所有label的估計機率
’‘’
#設定特徵數目
featrueNum
=
784
#設定類別數目
classNum
=
10
#建立存放所有標記的估計機率陣列
P
=
［
0
］
*
classNum
#對於每一個類別，單獨估計其機率
for
i
in
range
（
classNum
）：
#初始化sum為0，sum為求和項。
#在訓練過程中對機率進行了log處理，所以這裡原先應當是連乘所有機率，最後比較哪個機率最大
#但是當使用log處理時，連乘變成了累加，所以使用sum
sum
=
0
#獲取每一個條件機率值，進行累加
for
j
in
range
（
featrueNum
）：
sum
+=
Px_y
［
i
］［
j
］［
x
［
j
］］
#最後再和先驗機率相加（也就是式4。7中的先驗機率乘以後頭那些東西，乘法因為log全變成了加法）
P
［
i
］
=
sum
+
Py
［
i
］
#max（P）：找到機率最大值
#P。index（max（P））：找到該機率最大值對應的所有（索引值和標籤值相等）
return
P
。
index
（
max
（
P
））
def
test
（
Py
，
Px_y
，
testDataArr
，
testLabelArr
）：
‘’‘
對測試集進行測試
：param Py： 先驗機率分佈
：param Px_y： 條件機率分佈
：param testDataArr： 測試集資料
：param testLabelArr： 測試集標記
：return： 準確率
’‘’
#錯誤值計數
errorCnt
=
0
#迴圈遍歷測試集中的每一個樣本
for
i
in
range
（
len
（
testDataArr
））：
#獲取預測值
presict
=
NaiveBayes
（
Py
，
Px_y
，
testDataArr
［
i
］）
#與答案進行比較
if
presict
！=
testLabelArr
［
i
］：
#若錯誤 錯誤值計數加1
errorCnt
+=
1
#返回準確率
return
1
-
（
errorCnt
/
len
（
testDataArr
））
def
getAllProbability
（
trainDataArr
，
trainLabelArr
）：
‘’‘
透過訓練集計算先驗機率分佈和條件機率分佈
：param trainDataArr： 訓練資料集
：param trainLabelArr： 訓練標記集
：return： 先驗機率分佈和條件機率分佈
’‘’
#設定樣本特診數目，資料集中手寫圖片為28*28，轉換為向量是784維。
# （我們的資料集已經從影象轉換成784維的形式了，CSV格式內就是）
featureNum
=
784
#設定類別數目，0-9共十個類別
classNum
=
10
#初始化先驗機率分佈存放陣列，後續計算得到的P（Y = 0）放在Py［0］中，以此類推
#資料長度為10行1列
Py
=
np
。
zeros
（（
classNum
，
1
））
#對每個類別進行一次迴圈，分別計算它們的先驗機率分佈
#計算公式為書中“4。2節 樸素貝葉斯法的引數估計 公式4。8”
for
i
in
range
（
classNum
）：
#下方式子拆開分析
#np。mat（trainLabelArr） == i：將標籤轉換為矩陣形式，裡面的每一位與i比較，若相等，該位變為Ture，反之False
#np。sum（np。mat（trainLabelArr） == i）：計算上一步得到的矩陣中Ture的個數，進行求和（直觀上就是找所有label中有多少個
#為i的標記，求得4。8式P（Y = Ck）中的分子）
#np。sum（np。mat（trainLabelArr） == i）） + 1：參考“4。2。3節 貝葉斯估計”，例如若資料集總不存在y=1的標記，也就是說
#手寫資料集中沒有1這張圖，那麼如果不加1，由於沒有y=1，所以分子就會變成0，那麼在最後求後驗機率時這一項就變成了0，再
#和條件機率乘，結果同樣為0，不允許存在這種情況，所以分子加1，分母加上K（K為標籤可取的值數量，這裡有10個數，取值為10）
#參考公式4。11
#（len（trainLabelArr） + 10）：標籤集的總長度+10。
#（（np。sum（np。mat（trainLabelArr） == i）） + 1） / （len（trainLabelArr） + 10）：最後求得的先驗機率
Py
［
i
］
=
（（
np
。
sum
（
np
。
mat
（
trainLabelArr
）
==
i
））
+
1
）
/
（
len
（
trainLabelArr
）
+
10
）
#轉換為log對數形式
#log書中沒有寫到，但是實際中需要考慮到，原因是這樣：
#最後求後驗機率估計的時候，形式是各項的相乘（“4。1 樸素貝葉斯法的學習” 式4。7），這裡存在兩個問題：1。某一項為0時，結果為0。
#這個問題透過分子和分母加上一個相應的數可以排除，前面已經做好了處理。2。如果特診特別多（例如在這裡，需要連乘的專案有784個特徵
#加一個先驗機率分佈一共795項相乘，所有數都是0-1之間，結果一定是一個很小的接近0的數。）理論上可以透過結果的大小值判斷， 但在
#程式執行中很可能會向下溢位無法比較，因為值太小了。所以人為把值進行log處理。log在定義域內是一個遞增函式，也就是說log（x）中，
#x越大，log也就越大，單調性和原資料保持一致。所以加上log對結果沒有影響。此外連乘項透過log以後，可以變成各項累加，簡化了計算。
#在似然函式中通常會使用log的方式進行處理
Py
=
np
。
log
（
Py
）
#計算條件機率 Px_y=P（X=x|Y = y）
#計算條件機率分成了兩個步驟，下方第一個大for迴圈用於累加，參考書中“4。2。3 貝葉斯估計 式4。10”，下方第一個大for迴圈內部是
#用於計算式4。10的分子，至於分子的+1以及分母的計算在下方第二個大For內
#初始化為全0矩陣，用於存放所有情況下的條件機率
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，計算式4。10的分母，以及分子和分母之間的除法
#迴圈每一個標記（共10個）
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
（）
# 獲取訓練集
（
‘start read transSet’
）
trainDataArr
，
trainLabelArr
=
loadData
（
‘。。/Mnist/mnist_train。csv’
）
# 獲取測試集
（
‘start read testSet’
）
testDataArr
，
testLabelArr
=
loadData
（
‘。。/Mnist/mnist_test。csv’
）
#開始訓練，學習先驗機率分佈和條件機率分佈
（
‘start to train’
）
Py
，
Px_y
=
getAllProbability
（
trainDataArr
，
trainLabelArr
）
#使用習得的先驗機率分佈和條件機率分佈對測試集進行測試
（
‘start to test’
）
accuracy
=
test
（
Py
，
Px_y
，
testDataArr
，
testLabelArr
）
#列印準確率
（
‘the accuracy is：’
，
accuracy
）
#列印時間
（
‘time span：’
，
time
。
time
（）
-
start
）
