在《機器學習---朴素貝葉斯分類器（Machine Learning Naive Bayes Classifier）》一文中，我們介紹了朴素貝葉斯分類器的原理。現在，讓我們來實踐一下。
在這里，我們使用一份皮馬印第安女性的醫學數據，用來預測其是否會得糖尿病。文件一共有768個樣本，我們先剔除缺失值，然后選出20%的樣本作為測試樣本。
文件下載地址：https://raw.githubusercontent.com/jbrownlee/Datasets/master/pima-indians-diabetes.data.csv
特征分別是：
懷孕次數
口服葡萄糖耐量試驗中血漿葡萄糖濃度
舒張壓（mm Hg）
三頭肌組織褶厚度（mm）
2小時血清胰島素（μU/ ml）
體重指數（kg/（身高(m)）^ 2）
糖尿病系統功能
年齡（歲）
（注：特征值為0表示缺失值）
標簽是：
是否患有糖尿病，0代表沒有糖尿病，1代表患有糖尿病
在開始之前，先回顧一下朴素貝葉斯算法：
假設訓練集特征數為i，記作xi，目標為k個類別，記為Ck，樣本數為n，新樣本特征記作xnew_sample。
1，計算出先驗概率P(Ck)，即每個類別在訓練集中的概率；
2，分別計算出訓練集中每個特征在每個類別下的條件概率P(xi|Ck)，具體分為以下三種情況；
a）如果特征數據是離散值，那么我們假設其符合多項式分布，代入公式
；
b）如果特征數據是布爾類型的離散值，那么我們假設其符合伯努利分布，代入公式
；
c）如果特征數據是連續值，那么我們假設其符合高斯分布。只需要求出這個特征在每個類別下的平均值μk和方差σk2，然后代入公式
；
3，對於需要預測類別的新樣本，分別計算每個類別下的極大后驗概率：
（我們在實際計算中為防止數據下溢，將連乘運算取對數變成相加運算）；
4，對上述極大后驗概率進行比較，最大的那個對應的類別即為結果；
首先導入pandas和numpy庫，讀取csv文件，然后去除缺失值：
import pandas as pd import numpy as np data=pd.read_csv(r"C:\Users\ccav\pima-indians-diabetes.data.csv",header=None,\ names=["懷孕次數","口服葡萄糖耐量試驗中血漿葡萄糖濃度",\ "舒張壓（mm Hg）","三頭肌組織褶厚度（mm）",\ "2小時血清胰島素（μU/ ml）","體重指數（kg/（身高(m)）^ 2）",\ "糖尿病系統功能","年齡（歲）","是否患有糖尿病"]) data.iloc[:,0:8]=data.iloc[:,0:8].applymap(lambda x:np.NaN if x==0 else x) #把屬性值為0的地方轉換成NaN data=data.dropna(how="any",axis=0) #去除有缺失值的行
經過處理還只剩336個有效樣本：
print(len(data)) 336
通過查看，我們發現在這336個有效樣本中，225個是沒有糖尿病人的樣本，111個是患有糖尿病人的樣本，顯然原數據里兩個類別的樣本是不均衡的。這個問題可以通過不同的采樣方法解決，但是鑒於簡便起見，我們隨機選取大約20%的樣本用於測試：
#隨機選取80%的樣本作為訓練樣本 data_train=data.sample(frac=0.8,random_state=4,axis=0) #剩下的作為測試樣本 test_idx=[i for i in data.index.values if i not in data_train.index.values] data_test=data.loc[test_idx,:]
由於我們一共有8個特征，2個類別，因此，我們一共需要計算16個條件概率以及2個先驗概率，然后根據這些數據計算出2個極大后驗概率。該如何計算這些數據，又該如何合並，如果不事先想清楚怎么做，很容易被搞暈。此外，考慮到數據是連續值，我們假設其符合正態分布，因此在訓練時需要計算出每個特征的平均值和方差。
讓我們先來理清一下思路：
1，按類別分隔數據
2，獲取類別總數和類別名稱
3，訓練數據：計算每個類別的先驗概率，計算每個類別每個特征的平均值和方差
4，預測：計算每個類別每個特征的條件概率，計算每個類別的極大后驗概率並合並，得出最大可能的類別名稱
按照這個步驟，我們發現以下幾項都是需要復用的：按類別分隔數據，計算每個類別的先驗概率，計算每個類別每個特征的平均值和方差，計算每個類別每個特征的條件概率。因此，我們把這幾項分別做成function。
按類別分隔數據：
def SepByClass(X, y): ###按類別分隔數據### ###輸入未分類的特征和目標，輸出分類完成的數據（字典形式）### num_of_samples=len(y) #總樣本數 y=y.reshape(X.shape[0],1) data=np.hstack((X,y)) #把特征和目標合並成完整數據 data_byclass={} #初始化分類數據，為一個空字典 #提取各類別數據，字典的鍵為類別名，值為對應的分類數據 for i in range(len(data[:,-1])): if i in data[:,-1]: data_byclass[i]=data[data[:,-1]==i] class_name=list(data_byclass.keys()) #類別名 num_of_class=len(data_byclass.keys()) #類別總數 return data_byclass
計算每個類別的先驗概率：
def CalPriorProb(y_byclass): ###計算y的先驗概率（使用拉普拉斯平滑）### ###輸入當前類別下的目標，輸出該目標的先驗概率### #計算公式：（當前類別下的樣本數+1）/（總樣本數+類別總數） return (len(y_byclass)+1)/(num_of_samples+num_of_class)
計算每個類別每個特征的平均值和方差：
def CalXMean(X_byclass): ###計算各類別特征各維度的平均值### ###輸入當前類別下的特征，輸出該特征各個維度的平均值### X_mean=[] for i in range(X_byclass.shape[1]): X_mean.append(np.mean(X_byclass[:,i])) return X_mean def CalXVar(X_byclass): ###計算各類別特征各維度的方差### ###輸入當前類別下的特征，輸出該特征各個維度的方差### X_var=[] for i in range(X_byclass.shape[1]): X_var.append(np.var(X_byclass[:,i])) return X_var
計算每個類別每個特征的條件概率：
def CalGaussianProb(X_new, mean, var): ###計算訓練集特征（符合正態分布）在各類別下的條件概率### ###輸入新樣本的特征，訓練集特征的平均值和方差，輸出新樣本的特征在相應訓練集中的分布概率### #計算公式：(np.exp(-(X_new-mean)**2/(2*var)))*(1/np.sqrt(2*np.pi*var)) gaussian_prob=[] for a,b,c in zip(X_new, mean, var): formula1=np.exp(-(a-b)**2/(2*c)) formula2=1/np.sqrt(2*np.pi*c) gaussian_prob.append(formula2*formula1) return gaussian_prob
接下來，我們開始訓練數據。首先按類別分隔數據，然后遍歷每一個類別，分別計算每個類別的先驗概率和每個類別每個特征的平均值和方差，並儲存在列表中。
def fit(X, y): ###訓練數據### ###輸入訓練集特征和目標，輸出目標的先驗概率，特征的平均值和方差### #將輸入的X,y轉換為numpy數組 X, y = np.asarray(X, np.float32), np.asarray(y, np.float32) data_byclass=Gaussian_NB.SepByClass(X,y) #將數據分類 #計算各類別數據的目標先驗概率，特征平均值和方差 for data in data_byclass.values(): X_byclass=data[:,:-1] y_byclass=data[:,-1] prior_prob.append(Gaussian_NB.CalPriorProb(y_byclass)) X_mean.append(Gaussian_NB.CalXMean(X_byclass)) X_var.append(Gaussian_NB.CalXVar(X_byclass)) return prior_prob, X_mean, X_var
最后，輸入一個新樣本的特征，預測其所屬的類別。首先，遍歷之前在訓練數據時計算出的每個類別的先驗概率和每個類別每個特征的平均值和方差，再以此計算條件概率和極大后驗概率。選出最大的極大后驗概率所對應的索引，最后提取其對應的類別名稱。
def predict(X_new): ###預測數據### ###輸入新樣本的特征，輸出新樣本最有可能的目標### #將輸入的x_new轉換為numpy數組 X_new=np.asarray(X_new, np.float32) posteriori_prob=[] #初始化極大后驗概率 for i,j,o in zip(prior_prob, X_mean, X_var): gaussian=Gaussian_NB.CalGaussianProb(X_new,j,o) posteriori_prob.append(np.log(i)+sum(np.log(gaussian))) idx=np.argmax(posteriori_prob) return class_name[idx]
整理一下以上代碼，我們把這個高斯朴素貝葉斯分類器做成一個類，完整代碼如下：
import pandas as pd import numpy as np data=pd.read_csv(r"C:\Users\ccav\pima-indians-diabetes.data.csv",header=None,\ names=["懷孕次數","口服葡萄糖耐量試驗中血漿葡萄糖濃度",\ "舒張壓（mm Hg）","三頭肌組織褶厚度（mm）",\ "2小時血清胰島素（μU/ ml）","體重指數（kg/（身高(m)）^ 2）",\ "糖尿病系統功能","年齡（歲）","是否患有糖尿病"]) data.iloc[:,0:8]=data.iloc[:,0:8].applymap(lambda x:np.NaN if x==0 else x) #把屬性值為0的地方轉換成NaN data=data.dropna(how="any",axis=0) #去除有缺失值的行 #隨機選取80%的樣本作為訓練樣本 data_train=data.sample(frac=0.8,random_state=4,axis=0) #剩下的作為測試樣本 test_idx=[i for i in data.index.values if i not in data_train.index.values] data_test=data.loc[test_idx,:] #提取訓練集和測試集的特征和目標 X_train=data_train.iloc[:,:-1] y_train=data_train.iloc[:,-1] X_test=data_test.iloc[:,:-1] y_test=data_test.iloc[:,-1] class Gaussian_NB: def __init__(self): self.num_of_samples = None self.num_of_class = None self.class_name = [] self.prior_prob = [] self.X_mean = [] self.X_var = [] def SepByClass(self, X, y): ###按類別分隔數據### ###輸入未分類的特征和目標，輸出分類完成的數據（字典形式）### self.num_of_samples=len(y) #總樣本數 y=y.reshape(X.shape[0],1) data=np.hstack((X,y)) #把特征和目標合並成完整數據 data_byclass={} #初始化分類數據，為一個空字典 #提取各類別數據，字典的鍵為類別名，值為對應的分類數據 for i in range(len(data[:,-1])): if i in data[:,-1]: data_byclass[i]=data[data[:,-1]==i] self.class_name=list(data_byclass.keys()) #類別名 self.num_of_class=len(data_byclass.keys()) #類別總數 return data_byclass def CalPriorProb(self, y_byclass): ###計算y的先驗概率（使用拉普拉斯平滑）### ###輸入當前類別下的目標，輸出該目標的先驗概率### #計算公式：（當前類別下的樣本數+1）/（總樣本數+類別總數） return (len(y_byclass)+1)/(self.num_of_samples+self.num_of_class) def CalXMean(self, X_byclass): ###計算各類別特征各維度的平均值### ###輸入當前類別下的特征，輸出該特征各個維度的平均值### X_mean=[] for i in range(X_byclass.shape[1]): X_mean.append(np.mean(X_byclass[:,i])) return X_mean def CalXVar(self, X_byclass): ###計算各類別特征各維度的方差### ###輸入當前類別下的特征，輸出該特征各個維度的方差### X_var=[] for i in range(X_byclass.shape[1]): X_var.append(np.var(X_byclass[:,i])) return X_var def CalGaussianProb(self, X_new, mean, var): ###計算訓練集特征（符合正態分布）在各類別下的條件概率### ###輸入新樣本的特征，訓練集特征的平均值和方差，輸出新樣本的特征在相應訓練集中的分布概率### #計算公式：(np.exp(-(X_new-mean)**2/(2*var)))*(1/np.sqrt(2*np.pi*var)) gaussian_prob=[] for a,b,c in zip(X_new, mean, var): formula1=np.exp(-(a-b)**2/(2*c)) formula2=1/np.sqrt(2*np.pi*c) gaussian_prob.append(formula2*formula1) return gaussian_prob def fit(self, X, y): ###訓練數據### ###輸入訓練集特征和目標，輸出目標的先驗概率，特征的平均值和方差### #將輸入的X,y轉換為numpy數組 X, y = np.asarray(X, np.float32), np.asarray(y, np.float32) data_byclass=Gaussian_NB.SepByClass(X,y) #將數據分類 #計算各類別數據的目標先驗概率，特征平均值和方差 for data in data_byclass.values(): X_byclass=data[:,:-1] y_byclass=data[:,-1] self.prior_prob.append(Gaussian_NB.CalPriorProb(y_byclass)) self.X_mean.append(Gaussian_NB.CalXMean(X_byclass)) self.X_var.append(Gaussian_NB.CalXVar(X_byclass)) return self.prior_prob, self.X_mean, self.X_var def predict(self,X_new): ###預測數據### ###輸入新樣本的特征，輸出新樣本最有可能的目標### #將輸入的x_new轉換為numpy數組 X_new=np.asarray(X_new, np.float32) posteriori_prob=[] #初始化極大后驗概率 for i,j,o in zip(self.prior_prob, self.X_mean, self.X_var): gaussian=Gaussian_NB.CalGaussianProb(X_new,j,o) posteriori_prob.append(np.log(i)+sum(np.log(gaussian))) idx=np.argmax(posteriori_prob) return self.class_name[idx] if __name__=="__main__": Gaussian_NB=Gaussian_NB() #實例化Gaussian_NB Gaussian_NB.fit(X_train,y_train) #使用Gaussian_NB模型訓練數據 acc=0 TP=0 FP=0 FN=0 for i in range(len(X_test)): predict=Gaussian_NB.predict(X_test.iloc[i,:]) target=np.array(y_test)[i] if predict==1 and target==1: TP+=1 if predict==0 and target==1: FP+=1 if predict==target: acc+=1 if predict==1 and target==0: FN+=1 print("准確率:",acc/len(X_test)) print("查准率:",TP/(TP+FP)) print("查全率:",TP/(TP+FN)) print("F1:",2*TP/(2*TP+FP+FN))
結果如下：
准確率: 0.7611940298507462 查准率: 0.7391304347826086 查全率: 0.6296296296296297 F1: 0.68
由於樣本不均衡，因此光看准確率是不夠的，在這里一並計算了查准率和查全率以及F1值。可以看出這個分類器的效果不是很好，一方面是因為樣本數量較少，另一方面朴素貝葉斯的預測能力確實比不上復雜模型，只能提供一個粗略的判斷。