7.4 多项式朴素贝叶斯原理与实现#
在上一节内容中，我们详细介绍了一种常见的朴素贝叶斯算法，也被称之为Categorical Naive Bayes。但实际上，”朴素贝叶斯“算法远不止这一种，而它们之间的主要区别在于对条件概率的处理上[3]，即式(7-10)中的$\prod{P}({{X}^{(i)}}={{x}^{(i)}}|Y={{c}_{k}})$部分。因此在接下来的这节内容中，我们将会介绍第2种基于朴素贝叶斯思想的分类模型，多项朴素贝叶斯（Multinomial Naive Bayes, MNB）。
7.4.1 算法思想#
在通过Categorical NB来进行文本分类的场景中，在计算条件概率时都是将词表中的每个词以是否出现为标准进行类别化（Categorization）处理，因此如果将词频作为特征维度的取值类别，那么将会出现在测试集中特征维度的取值情况数大于训练集中的情况。
例如在训练集中“客栈”这个词出现的最大次数为10，那么模型在拟合过程中就会认为“客栈”这个维度的特征取值有10种情况，并以此进行建模；但是当测试集中的某个样本里“客栈”这个词出现的频次为11时，那么模型便会认为该维度多了一种取值情况，进而无法取到对应的条件概率。
同时，在利用词袋模型对文本进行向量化表示时词频也是一个重要的考量因素，而多项朴素贝叶斯算法在处理这一问题时则是将每个维度的词频在总词频中的占比来作为条件概率进行建模[4]。
7.4.2 算法原理#
在MNB中，我们可以将类别$c_k$下条件概率的分布参数化为$\theta_{c_k}=(\theta_{c_k1},\theta_{c_k2}...,\theta_{c_kn})$这样的形式，其中$n$表示训练集中的特征维度（例如在文本分类中则是词表的长度），$\theta_{c_ki}$则是类别$c_k$下特征$i$的条件概率。进一步，参数$\theta_{c_k}$可以通过极大似然估计来计算得到[3]：
$$ \hat{\theta}_{c_ki}=\frac{N_{c_ki}+\alpha_i}{N_{c_k}+\alpha }\tag{7-26} $$其中$N_{c_ki}$表示在整个训练集$T$中样本属于$c_k$这个类别下特征$x^{(i)}$出现的频次，即$N_{c_ki}=\sum_{x\in T}x^{(i)}$；$N_{c_k}$表示类别$c_k$下所有维度特征的总频次，即$N_{c_k}=\sum_{i=1}^nN_{c_ki}$；$\alpha_i$表示每个特征维度对应的平滑系数，$\alpha$表示所有平滑系数的总和[4]。但是在实际处理时通常会将每个维度的平滑系数设为相等，因此式(7-26)可以改写为：
$$ \hat{\theta}_{c_ki}=\frac{N_{c_ki}+\alpha}{N_{c_k}+\alpha n }\tag{7-27} $$其中$\alpha$表示每个特征维度的平滑系数。
在根据式(7-27)估计得到每个类别下各个特征的条件概率（词频占比）后，便可以通过式(7-10)来最大化后验概率以此确定样本的分类类别。但是这里存在一个问题，那就是通过式(7-10)可以知道，在最大化后验概率时各个特征维度的条件概率是进行的累乘操作，而在动则上千维的文本向量中，这样累乘计算得到的结果将会出现下溢的情况。
因此，常见的一种做法是在式(7-10)的两边同时取自然对数$\log$，且由于$\log$函数是单调的因此这并不影响最终的预测结果[5]。所以式(7-10)可以改写为如下形式
$$ \hat{y} = \arg\max_{c_k} \log{\left(P(Y=c_k)\prod\limits_{i=1}^{n}{P}({{X}^{(i)}}|Y={{c}_{k}})\right)}\tag{7-28} $$其中${P}({{X}^{(i)}}|Y={{c}_{k}})$表示在类别$Y=c_k$下特征$X^{(i)}$对应的条件概率。
进一步，根据$\log(xy)=\log(x)+\log(y)$，式(7-28)可以改写为
$$ \hat{y} = \arg\max_{c_k} \left(\log P(Y=c_k) +\sum_{i=1}^n\log P({{X}^{(i)}}|Y={{c}_{k}})\right)\tag{7-29} $$同时对于MNB算法来说，从式(7-27)可以看出，此时条件概率计算的是训练集中特征维度的词频占比（相当于模型参数），因此在最终计算后验概率时需要同时考虑到每个维度的词频情况，即[4]
$$ \hat{y} = \arg\max_{c_k} \left(\log P(Y=c_k) +\sum_{i=1}^nf_i\log P({{X}^{(i)}}|Y={{c}_{k}})\right)\tag{7-30} $$其中$f_i$表示特征维度$i$对应的词频。
此时根据式(7-30)的形式来看，还可以将$\log P({{X}^{(i)}}|Y={{c}_{k}})$理解为特征$x^{(i)}$在对应类别下重要性大小的权重，而先验概率$\log P(Y=c_k)$则可以理解为数据集中各类别的相对频次（偏置），频次越大则当前样本越可能归属于该类别。因此，从这个角度看还可以将多项贝叶斯理解一个简单的线性模型。也正是因为这样的特性，MNB算法在处理TFIDF这类文本特征表示时依旧有着很好的效果[3]。
到此，对于多项贝叶斯算法的基本原理就介绍完了，下面我们再来通过一个实际的计算示例来帮助大家更加清晰地理解。
7.4.3 计算示例#
假设现在有一批基于词袋模型表示的文本数据，其一共包含有$X^{(1)},X^{(2)},X^{(3)}$这3个特征维度，每个维度表示词表中相应词的词频，$Y$表示样本对应的所属类别，如表7-2所示。现需要预测$x=[17,25,39]$这个样本的所属类别。
根据式(7-30)，由表7-2易知，各个类别先验概率的$\log$取值为
$$ \begin{aligned} \log P(Y=0)&=\log(\frac{5}{15})\approx-1.098\\[1ex] \log P(Y=1)&=\log(\frac{7}{15})\approx-0.762\\[1ex] \log P(Y=2)&=\log(\frac{3}{15})\approx-1.609 \end{aligned}\tag{7-31} $$各类别下不同特征取值的条件概率的$\log$取值为（设$\alpha=1$，即拉普拉斯平滑）
$$ \log P(X^{(1)}|Y=0)=\log(\frac{2+1+3+0+0+1}{17+11+24+5+16+1\cdot3})\approx-2.384\tag{7-32} $$其中分子中的$2,1,3,0,0$表示在$Y=0$这个类别下，特征$X^{(1)}$出现的总频次；分母$17,11,24,5,16$分别表示在$Y=0$这个类别下所有特征维度的总频次，例如$17=2+6+9$。同时，从这里也可以再次看出MNB算法在计算特征维度时考虑的是类别$Y=c_k$中所有样本当前各个特征维度的频次总数在类别$Y=c_k$中所有样本的所有特征维度总频次中的占比，这样便能够计算得到对于类别$Y=c_k$来说哪个特征维度的重要性更高。
同理，对于其它情况来说，对应的条件概率的$\log$取值为
$$ \begin{aligned} \log P(X^{(2)}|Y=0)=\log(\frac{31+1}{73+1\cdot3})\approx-0.864\\[1ex] \log P(X^{(3)}|Y=0)=\log(\frac{36+1}{73+1\cdot3})\approx-0.719\\[1ex] \log P(X^{(1)}|Y=1)=\log(\frac{44+1}{94+1\cdot3})\approx-0.768\\[1ex] \log P(X^{(2)}|Y=1)=\log(\frac{22+1}{94+1\cdot3})\approx-1.439\\[1ex] \log P(X^{(3)}|Y=1)=\log(\frac{28+1}{94+1\cdot3})\approx-1.207\\[1ex] \log P(X^{(1)}|Y=2)=\log(\frac{19+1}{43+1\cdot3})\approx-0.832\\[1ex] \log P(X^{(2)}|Y=2)=\log(\frac{10+1}{43+1\cdot3})\approx-1.430\\[1ex] \log P(X^{(3)}|Y=2)=\log(\frac{14+1}{43+1\cdot3})\approx-1.120\\[1ex] \end{aligned}\tag{7-33} $$根据计算得到的先验概率、条件概率和式(7-30)，样本$x=[17,25,39]$的后验概率置信度为
$$ \begin{aligned} P(Y=0|x)&=\log P(Y=0)+X^{(1)}\cdot \log P(X^{(1)}|Y=0)\\[2ex] &+X^{(2)}\cdot\log P(X^{(2)}|Y=0)+X^{(3)}\cdot \log P(X^{(3)}|Y=0)\\[2ex] &=-(1.098+17\times2.384+25\times0.864+39\times0.719)\approx-91.267\\[1ex] \end{aligned}\tag{7-34} $$同理，其它两个类别的后验概率置信度分别为
$$ \begin{aligned} P(Y=1|x)&\approx-96.888\\[1ex] P(Y=2|x) &\approx-95.240\\[1ex] \end{aligned}\tag{7-35} $$根据上述结果可知，样本$x=[17,25,39]$应该属于类别$Y=0$。
7.4.4 多项式贝叶斯实现#
在有了前面Categorical贝叶斯算法的实现经验后，Multinomial贝叶斯的实现过程就非常容易理解了。下面，我们依旧分步进行讲解。以下实现代码均参考自sklearn 0.24.0 中的MultinomialNB模块，只是对部分处理逻辑进行了修改与简化，完整代码可参见AllBooKCode/Chapter07/C02_naive_bayes_multinomial.py 文件。
1. 特征计数实现
通过7.4.2节内容可知，不管是计算先验概率还是条件概率，在这之前都需要先统计训练集中各个样本及样本特征取值的分布情况。因此，这里首先需要初始化相关的计数器，然后再对样本和特征取值的分布情况进行统计。
对于两个重要计数器的初始化过程，示例代码如下：
1 class MyMultinomialNB(object):
2 def __init__(self, alpha=0):
3 self.alpha = alpha
4 self._ALPHA_MIN = 1e-10
5 def _init_counters(self, ):
6 self.class_count_ = np.zeros(self.n_classes, dtype=np.float64)
7 self.feature_count_ = np.zeros((self.n_classes, self.n_features_),dtype=np.float64)在上述代码中，第3~4行是初始化平滑项系数alpha。第6行class_count_被初始化成了一个形状为[n_classes,]的全零向量，其中n_classes表示分类的类别数量，而每个维度分别表示每个类别的样本数量（例如[2,2,3]表示0,1,2这三个类别的样本数分别是2,2,3），其目的是后续用于计算每个类别的先验概率。第7行feature_count_被初始化成了一个形状为[n_classes,n_features]的全零矩阵，其中feature_count_[i][j]表示，对于整个训练集来说，在第i个类别中第j个特征的出现频次。
例如对于表7-2中的示例来说，其计算后的feature_count_为
feature_count_ = [[ 6. 31. 36.],[44. 22. 28.],[19. 10. 14.]]在上述结果中，第1行6的含义便是在训练样本中$Y=0$这个类别下特征$X^{(1)}$的总频次为 $2+1+3+0+0=6$ ；第1行22的含义便是在训练样本中$Y=1$这个类别下特征$X^{(2)}$的总频次为 $11+10+0+0+0+1+0=22$。
在初始化完两个计数器之后，下面就可以来对样本类别和特征分布进行计数了，示例代码如下：
1 def _count(self, X, Y):
2 self.class_count_ += Y.sum(axis=0) # shape [n_classes,]
3 # 计算得到每个类别下的样本数量
4 self.feature_count_ += np.dot(Y.T, X) # [n_classes,n] @ [n,n_features_]在上述代码中，第1行参数Y是原始标签经过one-hot编码后的形式，例如3分类问题中类别1会被编码成[0,1,0]的形式，因此Y的形状为[n,n_classes]。第2行代码是计算得到每个类别对应的样本数量，形状为[n_classes,]。第4行则是统计每个类别下所有样本在各个特征维度下的取值频次。
进一步，在计算得到训练集中样本及特征的分布情况后，便可以计算先验概率和条件概率。
2. 先验概率实现
先验概率实现较为简单，整体实现代码如下：
1 def _update_class_prior(self):
2 log_class_count = np.log(self.class_count_)
3 self.class_prior_ = log_class_count - np.log(self.class_count_.sum())在上述代码中，第2~3便是用来计算各个类别的先验概率，并同时进行了取自然对数处理。
3. 条件概率实现
根据式(7-27)可知，计算条件概率的实现过程如下：
1 def _update_feature_prob(self, ):
2 smoothed_fc = self.feature_count_ + self.alpha
3 smoothed_cc = smoothed_fc.sum(axis=1)
4 self.feature_prob_ = (np.log(smoothed_fc) -
5 np.log(smoothed_cc.reshape(-1, 1)))从上面第1步最后的示例结果可知，feature_count_中每一行求和便是第i个类别下所有特征频次的总和，这对应的便是上述代码第3行。第4~5行代码则是根据式(7-27)来计算条件概率并同时取了对数。
4. 模型拟合实现
上述先验概率和条件概率的计算过程对应的便是整个模型的拟合过程，换句话说对于贝叶斯算法来说，所谓模型拟合就是计算先验概率和条件概率。在实现这部分代码之后，再通过一个函数将整个过程结合来即可，示例代码如下：
1 def fit(self, X, y):
2 self.n_features_ = X.shape[1]
3 labelbin = LabelBinarizer() # 将标签转化为one-hot形式
4 Y = labelbin.fit_transform(y) # one-hot 形式标签 shape: [n,n_classes]
5 self.classes_ = labelbin.classes_ # 原始标签类别 shape: [n_classes,]
6 if Y.shape[1] == 1: # 当数据集为二分类时fit_transform处理后的结果并不是one-hot形式
7 Y = np.concatenate((1 - Y, Y), axis=1) # 改变为one-hot形式
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~13行则分别是上面7.3节内容介绍到的初始化计数器、特征取值情况统计、计算先验概率和计算条件概率。
5. 后验概率实现
在完成模型的拟合过程后，对于新输入的样本来说其最终的预测结果则取决于对应的极大后验概率。根据式(7-30)可知，后验概率实现代码如下所示：
1 def _joint_likelihood(self, X):
2 return np.dot(X, self.feature_prob_.T) + self.class_prior_在上述代码中，第2行便是根据拟合得到的条件概率（权重）和先验概率（偏置）来计算得到样本的后验概率。
在实现每个样本后验概率的计算过程后，最后一步需要完成的便是极大化操作，即从所有后验概率中选择最大的概率值对应的类别作为该样本的预测类标即可，示例代码如下：
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样则是根据对应的参数来返回预测后的结果。
6. 使用示例
下面，我们先以表7-2中的模拟数据来通过上述实现代码进行示例。首先是构造模拟数据集，示例代码如下：
1 def load_simple_data():
2 import numpy as np
3 x = np.array([[5, 3, 2, 1, 0, 5, 12, 12, 10, 7, 8, 3, 0, 0, 1],
4 [11, 10, 6, 8, 7, 0, 0, 0, 0, 3, 1, 9, 1, 7, 0],
5 [7, 1, 9, 2, 12, 2, 15, 2, 2, 0, 0, 12, 4, 9, 1]]).transpose()
6 y = np.array([1, 1, 0, 0, 2, 1, 1, 2, 1, 2, 1, 0, 0, 0, 1])
7 return x, y进一步通过MyMultinomialNB来实现分类任务，实现代码如下所示：
1 def test_naive_bayes():
2 x, y = load_simple_data()
3 logging.info(f"MyMultinomialNB运行结果：")
4 model = MyMultinomialNB(alpha=1.)
5 model.fit(x, y)
6 logging.info(model.predict(np.array([[17, 25, 39]]), with_prob=True))
7 logging.info(f"MultinomialNB 运行结果：")
8 model = MultinomialNB(alpha=1.)
9 model.fit(x, y)
10 logging.info(model.predict(np.array([[17, 25, 39]])))
11 logging.info(model.predict_proba(np.array([[17, 25, 39]])))
12
13 if __name__ == '__main__':
14 test_naive_bayes()在上述代码中，第4~6行是根据实现的MyMultinomialNB来进行模型训练并预测。第8~11行则是利用sklearn中的MultinomialNB模块来完成模型的训练与预测工作。
上述代码运行结束后便可得到类似如下结果：
1 MyMultinomialNB运行结果：
2 各个类别下特征维度的频次为： [[ 6. 31. 36.],[44. 22. 28.], [19. 10. 14.]]
3 每个类别下的样本数class_count_:[5. 7. 3.]
4 计算每个类别的先验概率class_prior_: [-1.09861229 -0.76214005 -1.60943791]
5 各个类别下特征维度的占比为：
6 [[-2.38482319 -0.86499744 -0.71981543]
7 [-0.76804849 -1.43921676 -1.20741515]
8 [-0.83290912 -1.43074612 -1.1205912 ]]
9 样本预测原始概率为：[[-91.33834415 -96.88857422 -95.24060271]]
10 预测结果为：[0], [0.97648353, 0.00379516, 0.0197213 ]
11 MultinomialNB 运行结果：[0], [0.97648353 0.00379516 0.0197213 ]7.4.5 基于Multinomial的垃圾邮件分类#
下面，再通过一个垃圾邮件分类实例来测试上述代码的实现过程。首先同样需要构造数据集，这里采用考虑的词频的TFIDF来进行文本表示，示例代码如下：
1 def load_data():
2 x, y = load_cut_spam()
3 x_train, x_test, y_train, y_test = \
4 train_test_split(x, y, test_size=0.3, random_state=2020)
5 vect = TfidfVectorizer(max_features=1000)
6 x_train = vect.fit_transform(x_train)
7 x_test = vect.transform(x_test)
8 return x_train, x_test, y_train, y_test在上述代码中，第2行是载入分词后的原始样本。第5~7行则是先在训练集上得到词表并将训练样本转换为TFIDF表示，然后再用同一个词表将测试集上的训练样本转换为TFIDF表示，其中max_features取出现频率最高的1000个词。
进一步，完成整个垃圾邮件分类的建模任务，示例代码如下：
1 def test_spam_classification():
2 x_train, x_test, y_train, y_test = load_data()
3 logging.info(f"MyMultinomialNB 运行结果：")
4 model = MyMultinomialNB(alpha=1.)
5 model.fit(x_train.toarray(), y_train)
6 y_pred = model.predict(x_test.toarray())
7 logging.info(classification_report(y_pred, y_test))
8
9 logging.info(f"MultinomialNB 运行结果：")
10 model = MultinomialNB(alpha=1.)
11 model.fit(x_train, y_train)
12 y_pred = model.predict(x_test)
13 logging.info(classification_report(y_pred, y_test))在上述代码中，第4~5行是先实例化类MyMultinomialNB，然后再进行模型拟合，需要注意的是MyMultinomialNB并不支持稀疏表示，所以第5行里需要使用toarray()方法将x_train转换为稠密表示。
上述代码运行后便可以看到类似如下结果：
1 MyMultinomialNB 运行结果：
2 precision recall f1-score support
3 0 0.93 0.99 0.96 1418
4 1 0.99 0.94 0.96 1583
5 accuracy 0.96 3001
6 macro avg 0.96 0.96 0.96 3001
7 weighted avg 0.96 0.96 0.96 3001
8 MultinomialNB 运行结果：
9 precision recall f1-score support
10 0 0.93 0.99 0.96 1418
11 1 0.99 0.94 0.96 1583
12 accuracy 0.96 3001
13 macro avg 0.96 0.96 0.96 3001
14 weighted avg 0.96 0.96 0.96 3001从上述结果可以看到，两者在各项评估指标上并没有任何差异。
7.4.6 小结#
在本节内容中，我们首先介绍了多项式朴素贝叶斯的基本原理；然后通过一个实际的示例对其原理进行了进一步解释与阐述；接着一步一步详细地介绍了多项式朴素贝叶斯的实现方法；最后，我们分别以模拟数据和真实的文本数据来测试了实现代码的正确性。在下一节内容中，我们讲继续介绍第3种贝叶斯算法。