7.3 朴素贝叶斯实现#
经过前面两个小节内容的介绍，对于朴素贝叶斯算法的原理我们已经有了清晰的认识。在本节内容中，我们将开始分步对各个部分的实现进行详细地介绍。同时，需要说明的是以下实现代码均参考自sklearn 0.24.0 中的CategoricalNB模块，只是对部分处理逻辑进行了修改与简化，完整代码见 AllBooKCode/Chapter07/C01_naive_bayes_category.py 文件。
7.3.1 特征计数实现#
通过7.1节的内容可知，不管是计算先验概率还是条件概率，在这之前都需要先统计训练集中各个样本及样本特征取值的分布情况。因此，这里首先需要初始化相关的计数器，然后再对样本和特征取值的分布情况进行统计。
具体地，对于计数器的初始化工作实现过程，示例代码如下：
1 class MyCategoricalNB(object):
2
3 def __init__(self, alpha=1.0):
4 self.alpha = alpha
5
6 def _init_counters(self):
7 self.class_count_ = np.zeros(self.n_classes, dtype=np.float64)
8 self.category_count_ = [np.zeros((self.n_classes, 0))
9 for _ in range(self.n_features_)]在上述代码中，第3~4行是初始化平滑项系数alpha。第7行class_count_被初始化成了一个形状为[n_classes,]的全零向量，其中n_classes表示分类的类别数量，而每个维度分别表示每个类别的样本数量（例如[2,2,3]表示0、1、2这3个类别的样本数分别是2、2、3），其目的是后续用于计算每个类别的先验概率。第8行category_count_被初始化成了一个包含有n_features_个元素的列表，其中n_features_表示数据集的特征维度数量，同时category_count_中每个元素的形状是[n_classes,0]（后续每个元素将会更新为[n_classes,len(X_i)]的形状， len(X_i)表示X_i这个特征的取值情况数量）；而category_count_的作用是记录在各个类别下每个特征变量中各种取值情况的数量，例如category_count_[i][j][k]为10表示含义就是特征i在类别j下特征取值为k的样本数量为10个。
在初始化两个计数器之后，进一步便可以实现各个类别及特征分布的统计，示例代码如下：
1 def _count(self, X, Y):
2 def _update_cat_count(X_feature, Y, cat_count, n_classes):
3 for j in range(n_classes): # 遍历每个类别
4 mask = Y[:, j].astype(bool) # 取每个类别下对应样本的索引
5 counts = np.bincount(X_feature[mask]) # 统计当前类别下，特征X_feature中各个取值下的数量
6 indices = np.nonzero(counts)[0]
7 cat_count[j, indices] += counts[indices]
8
9 self.class_count_ += Y.sum(axis=0) # Y: shape(n,n_classes) Y.sum(): shape(n_classes,)
10 self.n_categories_ = X.max(axis=0) + 1
11 for i in range(self.n_features_): # 遍历每个特征
12 X_feature = X[:, i] # 取每一列的特征
13 self.category_count_[i] = np.pad(self.category_count_[i],
14 [(0, 0), (0, self.n_categories_[i])], 'constant')
15 _update_cat_count(X_feature, Y,self.category_count_[i],self.n_classes)在上述代码中，第1行参数Y是原始标签经过one-hot编码后的形式，例如3分类问题中类别1会被编码成[0,1,0]的形式，因此Y的形状为[n,n_classes]。第9行代码是计算得到每个类别对应的样本数量。第10行则是统计每个特征维度的取值数量（因为特征取值是从0开始的所以后面加了1）,例如[3 3 3 3]表示四个特征维度的取值均有3种情况。第11~12行开始遍历每个特征并取对应的特征列。第13~14行是对category_count_中的每个元素填充self.n_categories_[i]列全0向量，此时category_count_中每个元素将会变成形状为[n_classes,len(X_i)]的全零矩阵。第15行则是根据输入的每一列特征等相关参数来更新category_count_计数器。
第3~5行为遍历每一个样本类别，并取每个类别下对应样本的索引，同时统计当前类别下特征列X_feature中各个取值下的数量。同时，第5行中np.bincount的作用的是统计每个值出现的次数，例如：
1 counts = np.bincount(np.array([0, 3, 5, 1, 4, 4]))
2 print(counts) # [1 1 0 1 2 1]
3 # 表示[0, 3, 5, 1, 4, 4]中0,1,2,3,4,5这个6个值的出现的频次分别是1,1,0,1,2,1第6~7行则是用来更新cat_count中当前输入特征每种取值情况的分布数量，例如cat_count[i,k]表示的是第i个类别下，特征X_feature的第k个取值的数量。
例如对于表7-1中的数据集来说，其对应的category_count_计数器的结果为：
1 [[[4., 1.],[3., 7.]],
2 [[4., 1.],[4., 6.]],
3 [[1., 1., 3.],[2., 3., 5.]]]其中[[4., 1.],[3., 7.]]分别表示的含义就是：对于第1个特征$X^{(1)}$来说，在$Y=0$这个类别下，$X^{(1)}$取值为0的情况有4种（表中第4~7号样本），取值为1的情况有1种（第14号样本）；在$Y=1$这个类别下，$X^{(1)}$取值为0的情况有3种（第1~3号样本），取值为1的情况有7种（第8~13和第15号样本）。同理，对于第2个特征$X^{(2)}$来说，在$Y=0$这个类别下，$X^{(2)}$取值为0的情况有4种，取值为1的情况有1种；在$Y=1$这个类别下，$X^{(2)}$取值为0的情况有4种，取值为1的情况有6种。
到此，对于数据集中样本及特征分布情况的计数就算是统计完了，接下来开始实现计算先验概率和条件概率。
7.3.2 先验概率实现#
有了数据集中各个样本的分布情况后，计算先验概率就变得十分简单了，示例代码如下：
1 def _update_class_prior(self):
2 self.class_prior_ = (self.class_count_ + self.alpha) / # shape: [n_classes, ]
3 (self.class_count_.sum() + self.n_classes * self.alpha)在上述代码中，第2~3行便是用来计算各个类别的先验概率，其中带有alpha的地方便是第7.2.1节内容中的平滑处理项。
7.3.3 条件概率实现#
进一步，可以根据category_count_中的统计结果来实现各特征取值的条件概率，示例代码如下：
1 def _update_feature_prob(self):
2 feature_prob = []
3 for i in range(self.n_features_): # 遍历 每一个特征
4 smoothed_cat_count = self.category_count_[i] + self.alpha
5 smoothed_class_count = self.category_count_[i].sum(axis=1) +
6 self.category_count_[i].shape[1] * self.alpha
7 cond_prob = smoothed_cat_count / smoothed_class_count.reshape(-1, 1)
8 feature_prob.append(cond_prob)
9 self.feature_prob_ = feature_prob在上述代码中，第4行是取当前特征在不同类别下各个特征的取值的统计结果，并加上平滑项系数。第5~6行是计算在当前特征下各类别样本的总数，并进行平滑处理。第7行是计算每个特征取值对应的条件概率。第9行中feature_prob_ 为一个包含有n_features_个元素的列表，每个元素的shape为 [self.n_classes, 特征取值数]，因此feature_prob_[i][j][k]表示的就是特征i在类别j下取值为k时的概率。
例如对于表7-1中的训练数据集来说，计算得到feature_prob_的结果如下所示：
1 [[0.8 0.2] [[0.8 0.2] [[0.2 0.2 0.6]
2 [0.3 0.7]] [0.4 0.6]] [0.2 0.3 0.5]]上述结果中，第1个元素表示：对于第1个特征$X^{(1)}$来说，在$Y=0$这个类别下，$X^{(1)}=0$的概率为$0.8$，$X^{(1)}=1$的概率为$0.2$；对于第1个特征$X^{(1)}$来说，在$Y=1$这个类别下，$X^{(1)}=0$的概率为$0.3$，$X^{(1)}=1$的概率为$0.7$。这一结果也可以同7.1.3节中的计算结果进行对比。
7.3.4 模型拟合实现#
上述先验概率和条件概率的计算过程对应便是整个模型的拟合过程，换句话说对于贝叶斯算法来说，所谓模型拟合就是计算先验概率和条件概率。在实现这部分代码之后，再通过一个函数将整个过程整合起来即可，示例代码如下：
1 def fit(self, X, y):
2 self.n_features_ = X.shape[1]
3 labelbin = LabelBinarizer()
4 Y = labelbin.fit_transform(y)
5 self.classes_ = labelbin.classes_
6 if Y.shape[1] == 1:
7 Y = np.concatenate((1 - Y, Y), axis=1)
8 self.n_classes = Y.shape[1]
9 self._init_counters() # 初始化计数器
10 self._count(X, Y) # 对各个特征的取值情况进行计数，以计算条件概率等
11 self._update_class_prior()
12 self._update_feature_prob()
13 return self在上述代码中，第2~4行用来将原始[0,1,2,3...]这样的标签转换为one-hot编码形式的标签值，形状为[n,n_classes]。第5行用来记录原始的标签类别，形状为[n_classes,]。第6~7行用来处理当数据集为二分类时fit_transform处理后的结果并不是one-hot形式，需要添加一列来转换成one-hot形式。第8行是获取数据集的类别数量。第9~12行则分别是上面3节内容介绍到的初始化计数器、特征取值分布情况统计、计算先验概率和计算条件概率。
7.3.5 后验概率实现#
在完成模型的拟合过程后，对于新输入的样本来说其最终的预测结果则取决于对应的极大后验概率，这部分示例代码如下：
1 def _joint_likelihood(self, X):
2 if not X.shape[1] == self.n_features_:
3 raise ValueError("期望输入特征个数为 %d , 实际输入为 %d "
4 % (self.n_features_, X.shape[1]))
5 jll = np.ones((X.shape[0], self.class_count_.shape[0]))
6 for i in range(self.n_features_):
7 indices = X[:, i] # 取对应的每一列特征
8 if self.feature_prob_[i].shape[1] <= indices.max():
9 raise IndexError(f"测试集中的第{i}个特征维度的取值情况"
10 f" {indices.max()} 超出了训练集中该维度的取值情况！")
11 jll *= self.feature_prob_[i][:, indices].T
12 total_ll = jll * self.class_prior_
13 return total_ll在上述代码中，第2~4行用来判断输入样本的特征维度是否等于训练集中样本的特征维度。第5行是初始化每个样本的联合似然估计值。第6~7行是遍历测试数据中对应的每一列特征。第11行是取每个特征取值下对应类别的条件概率，并进行累乘（因为self.feature_prob_[i][:, indices]的形状为[n_classes,n]，而jll的初始形状为[n,n_classes]所以前者要进行转置）。第12行则是通过条件概率乘以先验概率来计算每个样本对应的后验概率估计。
同时值得注意的一点是，在文本向量化中对于考虑词频的词袋模型来说，其向量表示并不能直接使用于这里实现的贝叶斯模型。因为此时训练集词表中每个词出现的次数并不是该特征维度对应的取值情况数，而是表示的频次。例如在训练集中“客栈”这个词出现的最大次数为10，那么模型在拟合过程中就会认为“客栈”这个维度的特征取值有10种情况，并以此进行建模；但是当测试集中的某个样本里“客栈”这个词出现的频次为11时，那么模型便会认为该维度多了一种取值情况，进而无法取到对应的条件概率。所以我们在上面第8~10行中加了对应的判断条件。
在实现各样本后验概率的计算过程后，最后一步需要完成的便是极大化操作，即从所有后验概率中选择最大的概率值对应的类别作为该样本的预测类别即可。这部分示例代码如下：
1 def predict(self, X, with_prob=False):
2 from scipy.special import softmax
3 jll = self._joint_likelihood(X)
4 y_pred = self.classes_[np.argmax(jll, axis=1)]
5 if with_prob:
6 prob = softmax(jll)
7 return y_pred, prob
8 return y_pred在上述代码中，第4行便是极大化后验概率的操作。第5~8行则是根据对应的参数来返回预测后的结果是否包括概率值。
7.3.6 使用示例#
下面，我们先以表7-1中的模拟数据来通过上述代码进行示例。首先需要构建表7-1中的模拟数据，示例代码如下：
1 def load_simple_data():
2 x = np.array([[0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1],
3 [1, 1, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0],
4 [2, 1, 1, 2, 2, 2, 0, 2, 2, 0, 0, 2, 2, 1, 1]]).transpose()
5 y = np.array([1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 1])
6 return x, y其中特征$X^{(3)}$的取值这里使用了0,1,2来进行代替。
接着来实现模型的构建部分，代码如下：
1 def test_naive_bayes():
2 x, y = load_simple_data()
3 logging.info(f"My Bayes 运行结果：")
4 model = MyCategoricalNB(alpha=0)
5 model.fit(x, y)
6 logging.info(model.predict(np.array([[0, 1, 0]]), with_prob=True))
7 logging.info(f"CategoricalNB 运行结果：")
8 model = CategoricalNB(alpha=0)
9 model.fit(x, y)
10 logging.info(model.predict(np.array([[0, 1, 0]])))
11 logging.info(model.predict_proba(np.array([[0, 1, 0]])))
12 if __name__ == '__main__':
13 test_naive_bayes()在上述代码中，分别使用了我们自己动手实现的贝叶斯模型和sklearn中的CategoricalNB模块进行对比，其中logging为日志打印模块。
最终，上述代码运行以后将会输出类似如下结果：
1 MyCategoricalNB 运行结果：
2 每个类别下的样本数class_count_(n_classes,): [ 5. 10.]
3 每个特征的取值种数n_categories_:[2 2 3]
4 各个特征每个取值的数量分布（未平滑处理） category_count_:
5 [[[4., 1.],[3., 7.]], [[4., 1.],[4., 6.]],
6 [[1., 1., 3.], [2., 3., 5.]]]
7 n_classes:2 class_count_:[ 5. 10.]
8 计算每个类别的先验概率class_prior_:[0.333 0.667]
9 第0个特征在各类别下各特征取值的条件概率为: [[0.8 0.2],[0.3 0.7]]
10 第1个特征在各类别下各特征取值的条件概率为: [[0.8 0.2],[0.4 0.6]]
11 第2个特征在各类别下各特征取值的条件概率为: [[0.2 0.2 0.6],[0.2 0.3 0.5]]
12 样本预测原始概率为：[[0.011 0.024]]
13 [1] [[0.497, 0.503]]
14 CategoricalNB 运行结果：[1] [[0.307 0.692]]在上述结果中，第12行输出的便是样本[0, 1, 0]的原始预测概率，也就是7.1.3节中的计算结果。第13行则是最终的预测类别以及经过归一化后的概率输出。第14行则是sklearn中CategoricalNB模块输出的预测结果和概率值，其之所以与第13行中的概率值不一样是因为CategoricalNB在计算条件概率以及先验概率时都先取了$\log$对数，但是这并不影响最终的预测结果。不过从结果来看，取对数后两个类别预测的概率有更大的区分度。
下面，还可以通过7.2.3节中的垃圾邮件数据集来对上述实现进行验证，整个示例代码如下：
1 def test_spam_classification():
2 x_train, x_test, y_train, y_test = load_data()
3 model = MyCategoricalNB(alpha=1.0)
4 model.fit(x_train, y_train)
5 y_pred = model.predict(x_test)
6 logging.info(f"MyCategoricalNB 运行结果：")
7 logging.info(classification_report(y_test, y_pred))
8
9 model = CategoricalNB()
10 model.fit(x_train, y_train)
11 y_pred = model.predict(x_test)
12 logging.info(f"CategoricalNB 运行结果：")
13 logging.info(classification_report(y_test, y_pred))上述代码便是一个基于朴素贝叶斯模型的垃圾邮件分类示例，代码运行结束后便可以看到类似如下结果：
1 MyCategoricalNB 运行结果：
2 precision recall f1-score support
3 0 0.97 0.96 0.97 1504
4 1 0.96 0.97 0.97 1497
5 accuracy 0.97 3001
6 macro avg 0.97 0.97 0.97 3001
7 weighted avg 0.97 0.97 0.97 3001
8
9 CategoricalNB 运行结果：
10 precision recall f1-score support
11 0 0.97 0.96 0.97 1504
12 1 0.96 0.97 0.97 1497
13 accuracy 0.97 3001
14 macro avg 0.97 0.97 0.97 3001
15 weighted avg 0.97 0.97 0.97 3001可以看到，本文实现的MyCategoricalNB和CategoricalNB在各项评估指标上并没有任何差异。
7.3.7 小结#
在这篇文章中，我们首先详细介绍了朴素贝叶斯中样本分布及特征取值的统计、先验概率和条件概率的计算、极大化后验概率的代码实现等；然后一步一步介绍了如何从零实现不考虑词频的词袋模型；最后，再次以垃圾邮件这一分类任务来测试了实现代码的正确性，以及同第6.2节中K近邻的分类结果进行了对比。