Level3 — PART 4 机器学习算法

贝叶斯定理

朴素贝叶斯模型（Naive Bayes Model）

估计

离散估计

极大似然估计

案例

朴素贝叶斯扩展

高斯贝叶斯分类器

原理

应用

源码分析

伯努利贝叶斯分类器

原理

源码分析

多项朴素贝叶斯分类器

半朴素贝叶斯分类器

模拟题

CDA LEVEL III 模拟题（一）

CDA LEVEL III 模拟题（二）

贝叶斯定理

贝叶斯定理由英国数学家贝叶斯（Thomas Bayes 1702-1761）发展，用来描述两个条件概率之间的关系，比如 $P(A|B)$ 和 $P(B|A)$ ， $A$ 和 $B$ 是两个随机事件。

按照乘法法则，可以立刻导出：

$P(AB) = P(A)*P(B|A)=P(B)*P(A|B)$

其中 $P(AB)$ 表示随机事件 $A$ 和随机事件 $B$ 同时发生的概率， $P(B|A)$ 表示在事件 $A$ 发生的情况下，事件 $B$ 发生的概率，类似地， $P(A|B)$ 表示在事件 $B$ 发生的情况下，事件 $A$ 发生的概率。

如上公式也可变形为：

$P(A|B)=P(B|A)*\frac{P(A)}{P(B)}$

朴素贝叶斯模型（Naive Bayes Model）

假设有个数据，每一个数据由个特征构成，如下：

	特征1	特征2	...	特征n	类别
1	$x_{11}$	$x_{12}$	...	$x_{1n}$
2	$x_{21}$	$x_{22}$	...	$x_{2n}$
. . .	. . .	. . .		. . .	. . .
m	$x_{m1}$	$x_{m2}$	...	$x_{mn}$

记 $\mathbf{x}_i=(x_{i1},x_{i2},...,x_{in})$ 表示第个数据， $x_{ij}$ 为第i个数据的第j个特征的值， $y_i\in\{c_1,c_2,...,c_p\}$ 表示数据 $\mathbf{x}_i$ 的类别值，它是 $\{c_1,c_2,...,c_p\}$ 中某个类别值。

现在的任务是：给定一个由个特征构成的新数据 $\mathbf{x}=(x_1,x_2,...x_n)$ ，如何基于贝叶斯定理来对数据 $\mathbf{x}$ 进行分类呢？

对于每个类别 c_k $(k=1,2,...,p)$ ，在已知数据 $\mathbf{x}$ 的情况下，计算 $\mathbf{x}$ 属于类别 c_k 的概率，即条件概率 $P(c_k|\mathbf{x})$ 。这样会得到 $p$ 个条件概率 $P(c_k|\mathbf{x})$ 的最。其中条件概率最大值对应的类别就作为数据 $\mathbf{x}$ 所属的分类。这是合理的，因为概率 $P(c_k|\mathbf{x})$ 表示数据 $\mathbf{x}$ 属于类别 c_k 的可能性， $P(c_k|\mathbf{x})$ 值越大，说明数据 $\mathbf{x}$ 属于类别 c_k 的可能性也就越大。

那么如上的分类问题就转换为：在已知数据 $\mathbf{x}$ 和上面表格中m个数据以及它们的分类的情况下，求如下最大条件概率对应的类别的问题：

$c_k=argmax_{1\leq k\leq p}p(c_k|\mathbf{x})$

回顾上面的贝叶斯定理：

$P(A|B)=P(B|A)*\frac{P(A)}{P(B)}$

令 $A=c_k, B=\mathbf{x}$ ，贝叶斯公式转化为：

$P(c_k|\mathbf{x})=P(\mathbf{x}|c_k)*\frac{P(c_k)}{P(\mathbf{x})}$

使用上面已知的数据计算出公式右边的三个概率 $P(\mathbf{x}|c_k),P(c_k),P(\mathbf{x})$ 就可以计算出 $P(c_k|\mathbf{x})$ 。如何计算这些概率呢？对于条件概率 $P(\mathbf{x}|c_k)$ ，为了可以计算出该概率，我们需要做独立性假设，即假设每个特征之间是相互独立的，也即

$P(\mathbf{x})=P(x_1,x_2,...,x_n)=P(x_1)P(x_2)...P(x_n)=\prod_{i=1}^{n}P(x_i)$

这样

$P(\mathbf{x}|c_k)=P(x_1,x_2,...,x_n|c_k)=P(x_1|c_k)P(x_2|c_k)...P(x_n|c_k)$

$=\prod_{i=1}^{n}P(x_i|c_k)$

式中 $P(x_i|c_k)$ 表示类别为 $c_k$ 的情况下，第 $i$ 个特征值为 $x_i$ 的概率。因此条件概率可以如下计算：

$P(c_k|\mathbf{x})=\prod_{i=1}^{n}P(x_i|c_k)*\frac{P(c_k)}{P(\mathbf{x})}$

上面的问题转化为如下计算问题：

$c_k=argmax_{1\leq k\leq p}\prod_{i=1}^{n}P(x_i|c_k)\frac{P(c_k)}{P(\mathbf{x})}$

以上就是朴素贝叶斯模型的原理。之所以称其“朴素”，“Naive”，是因为我们对特征做了独立性假设：特征之间相互独立。

事实上，读者已经看到。对于 $k=1,2,...,p$ ， $P(\mathbf{x})$ 始终是固定的，因此 $P(\mathbf{x})$ 不影响条件概率 $P(c_k|\mathbf{x})$ 大小的比较，在具体计算时，无需计算该值。

估计

有了上面的贝叶斯模型，如何对 $P(\mathbf{x}|c_k),P(c_k),P(\mathbf{x})$ 进行估计，从而计算出 $P(c_k|\mathbf{x})$ 。这里介绍了两种估计方法：离散估计和极大似然估计。

离散估计

假设数据的每个特征取值范围都是有限的几个值，也就是说特征是类别数据，那么可以如下估计 P(c_k) 、 $P(\mathbf{x})$ 和 P(x_i|c_k) ：

（1）估计 P(c_k)

P(c_k) 表示的是类别 c_k 发生的概率，可以使用表格中类别 c_k 的频率与数据总个数的比值来近似估计。假设表中类别为 c_k 的数据个数为 r_k ，则可以如下估计出 P(c_k) ：

$P(c_k)=\frac{r_k}{m}$

（2）估计 $P(\mathbf{x})$

$P(\mathbf{x})$ 表示的是数据 $\mathbf{x}$ 发生的概率，由于数据的每个特征之间是相互独立的，所以

$P(\mathbf{x})=\prod_{i=1}^{n}P(x_i)$ ,

其中 P(x_i) 表示数据 $\mathbf{x}$ 的第个特征取值为 x_i 的概率，它可以使用 x_i 的频率与数据总个数的比值来代替，假设第个特征值为 x_i 的数据个数为 m_i ，则

$P(x_i)=\frac{m_i}{m}$ ，

则 $P(\mathbf{x})$ 可以如下估计：

$P(\mathbf{x})=\frac{1}{m^n}\prod_{i=1}^{n}m_i$

（3）估计 P(x_i|c_k)

P(x_i|c_k) 表示当类别为 c_k 的情况下，数据 $\mathbf{x}$ 的第个特征值为 x_i 的概率，假设在类别为 c_k 的数据中，第个特征值为 x_i 的数据个数为 q_i ，则可以如下估计：

$P(x_i|c_k)=\frac{q_i}{r_k}$

注：离散估计适用于数据特征为类别型的数据。

极大似然估计

所谓极大似然估计（Maximum Likelihood Estimation，简称MLE），不严谨地说，就是已经知道了随机变量 $x$ 服从的概率密度函数形式为 $f(x,\theta)$ ，但参数 $\theta$ 是未知的。然后通过实验或者其他方式获得了关于随机变量 $x$ 的 $m$ 个数据，比如 $x_1,x_2,...,x_m$ 。如何使用这些数据估计出参数 $\theta$ 的值呢？

首先需要回答的是为什么可以使用这些数据估计参数 $\theta$ ？这点应该不难理解，因为 $\theta$ 的取值不同，会影响获得的数据，也就是说 $\theta$ 在一定程度上决定了会产生这 $m$ 个数据，所以当然可以使用这些数据估计参数 $\theta$ 。那么，什么样的 $\theta$ 才会获得像 $x_1,x_2,...,x_m$ 这样形式的数据呢？从概率的角度来说，也就是 $\theta$ 为何值时，会使获得像 $x_1,x_2,...,x_m$ 这样的数据可能性更大呢？这就是极大似然估计的思想。使 $x_1,x_2,...,x_m$ 出现的可能性最大的 $\hat{\theta}$ ，就是 $\theta$ 的一个极大似然估计。

下面举个例子来理解下，比如最常见的“抛硬币”。对于一枚不均匀的硬币，假设抛掷后只会出现正面和反面，出现正面记作1，反面记作0，正面出现的概率为 $p$ ， $p$ 是未知的，它就是上面我们提到的 $\theta$ 。进一步假设每次抛硬币相互独立。虽然不知道概率 $p$ 的具体值，但是我们知道抛硬币出现的结果是一个随机变量，假设为 $x$ ，它服从如下的二项分布

$f(x,p)=p^x(1-p)^{1-x}$

通过抛掷10次硬币，获得了如下10个数据：

1, 0, 0, 1, 1, 1, 1, 0, 1, 1

下面使用这10个数据估计出概率 $p$ 。用 $A$ 表示这组数据， $P(A|p)$ 表示在概率为 $p$ 的情况下，出现数据组 $A$ 的概率。根据极大似然估计思想，需要求解的就是使 $P(A|p)$ 最大的 $p$ 值。但是首先需要给出 $P(A|p)$ 的具体表达式，因为每次抛掷硬币是相互独立的，所以 $P(A|p)$ 可以如下表示：

$P(A|p)=\prod_{i=1}^{10}P(x_i|p)=\prod_{i=1}^{10}f(x_i,p)=\prod_{i=1}^{10}p^{x_i}(1-p)^{1-x_i}$ ，

这就是所谓的极大似然函数。

因为 $P(A|p)$ 的最大值和对数 $lnP(A|p)$ 的最大值是等价的，所以上式两边取对数，得到：

$lnP(A|p)=\sum_{i=1}^{10}ln[p^{x_i}(1-p)^{1-x_i}]=\sum_{i=1}^{10}[x_ilnp+(1-x_i)ln(1-p)]$

两边对 $p$ 求导，得到

$\frac{\partial lnP(A|p)}{\partial p}=\sum_{i=1}^{10}(\frac{x_i}{p}-\frac{1-x_i}{1-p})$

令导数为0，得到：

$p=\frac{1}{10}\sum_{i=1}^{10}x_i=0.7$ .

也就是说当概率 $p=0.7$ 时，出现上面10个数据的可能性最大。这与实际得到的10个数据是相符的，从上面数据就可以看到，10个数据中，有7次是正面朝上。其实概率 $p$ 的极大似然估计就是获得的10个数据的平均值。概率 $p=0.7$ 有从侧面验证了硬币确实是不均匀的。

注：极大似然估计适用于已知特征数据服从某个形式的概率密度函数。

案例

表中是飞机延误相关的训练数据，总计14条数据。

Season	Weather	A_Control	Airline	Delay
Summer	Sunny	no	CZ	no
Summer	Sunny	no	CA	no
Autumn	Sunny	no	CZ	yes
WinterSpring	RainyOrSnowy	no	CZ	yes
WinterSpring	Cloudy	yes	CZ	yes
WinterSpring	Cloudy	yes	CA	no
Autumn	Cloudy	yes	CA	yes
Summer	RainyOrSnowy	no	CZ	no
Summer	Cloudy	yes	CZ	yes
WinterSpring	RainyOrSnowy	yes	CZ	yes
Summer	RainyOrSnowy	yes	CA	yes
Autumn	RainyOrSnowy	no	CA	yes
Autumn	Sunny	yes	CZ	yes
WinterSpring	RainyOrSnowy	no	CA	no

它包含4个特征列：Season、Weather、A_Control、Airline，1个类别列（或标签列）Delay。其中

Season有三种不同类型的取值：Summer、Autumn和WinterSpring

Weather有三种不同类型的取值：Sunny、Cloudy和RainyOrSnowy

A_Control有两种不同类型的取值：no和yes，

Airline有两种不同类型的取值：CA和CZ。

基于上述训练数据，预测新样本

x=<Season=Summer，Weather=RainyOrSnowy，A_Control=no，Airline=CA>

是否会延误，也就是预测Delay的值是no还是yes？

从样本数据，可以看到这些特征数据是类别型数据，可以使用离散估计来计算。回顾上面的贝叶斯模型

$P(c_k|\mathbf{x})=\prod_{i=1}^{n}P(x_i|c_k)*\frac{P(c_k)}{P(\mathbf{x})}$

对于当前案列， $c_k=Delay$ ，

$\mathbf{x}$ =<Season=Summer,Weather=RainyOrSnowy,A_Control=no,Airline=CA>

（1）估算 $P(\mathbf{x})$

也就是估算

P(Season=Summer,Weather=RainyOrSnowy,A_Control=no,Airline=CA)

=P(Season=Summer) P(Weather=RainyOrSnowy) P(A_Control=no) P(Airline=CA)

= $\frac{5}{14}*\frac{6}{14}*\frac{7}{14}*\frac{6}{14}$

（2）估算 P(c_k)

也就是估算 $P(Delay)$ ，

$P(Delay=yes)=\frac{9}{14}$ ， $P(Delay=no)=\frac{5}{14}$

（3）估算 $\prod_{i=1}^{n}P(x_i|c_k)$

Delay取值yes的情况:

P(Season=Summer|Delay=yes)= $\frac{2}{9}$

P(Weather=RainyOrSnowy|Delay=yes)= $\frac{4}{9}$

P(A_Control=no|Delay=yes)= $\frac{3}{9}$

P(Airline=CA|Delay=yes)= $\frac{3}{9}$

于是得到

$\prod_{i=1}^{n}P(x_i|c_k)=\frac{2}{9}*\frac{4}{9}*\frac{3}{9}*\frac{3}{9}$

Delay取值no的情况:

P(Season=Summer|Delay=yes)= $\frac{3}{5}$

P(Weather=RainyOrSnowy|Delay=yes)= $\frac{2}{5}$

P(A_Control=no|Delay=yes)= $\frac{4}{5}$

P(Airline=CA|Delay=yes)= $\frac{3}{5}$

于是得到

$\prod_{i=1}^{n}P(x_i|c_k)=\frac{3}{5}*\frac{2}{5}*\frac{4}{5}*\frac{3}{5}$

（4）计算 $P(c_k|\mathbf{x})$

$P(c_k|\mathbf{x})=\prod_{i=1}^{n}P(x_i|c_k)*\frac{P(c_k)}{P(\mathbf{x})}$

其中 $c_k=Delay$ ，

$\mathbf{x}$ =<Season=Summer,Weather=RainyOrSnowy,A_Control=no,Airline=CA>

Delay取值yes：

$\prod_{i=1}^{n}P(x_i|c_k)=\frac{2}{9}*\frac{4}{9}*\frac{3}{9}*\frac{3}{9}$

$P(c_k)=P(Delay=yes)=\frac{9}{14}$

$P(\mathbf{x})=\frac{5}{14}*\frac{6}{14}*\frac{7}{14}*\frac{6}{14}$

$P(c_k|\mathbf{x})=\prod_{i=1}^{n}P(x_i|c_k)*\frac{P(c_k)}{P(\mathbf{x})}=0.2151$

Delay取值no：

$\prod_{i=1}^{n}P(x_i|c_k)=\frac{3}{5}*\frac{2}{5}*\frac{4}{5}*\frac{3}{5}$

$P(c_k)=P(Delay=no)=\frac{5}{14}$

$P(\mathbf{x})=\frac{5}{14}*\frac{6}{14}*\frac{7}{14}*\frac{6}{14}$

$P(c_k|\mathbf{x})=\prod_{i=1}^{n}P(x_i|c_k)*\frac{P(c_k)}{P(\mathbf{x})}=1.2544$

有读者可能会怀疑计算有误，概率怎么会大于1。原因就在于我们做了独立性假设，而该案例中特征之间并不都是相互独立的，例如Season和Weather存在关联性，通常夏天的天气可能会是Sunny，而春季和冬季的天气可能会是RainyOrSnowy。

从计算结果来看，新样本

x=<Season=Summer，Weather=RainyOrSnowy，A_Control=no，Airline=CA>

的Delay预测值为no，也就是不会延误。

手工计算总是不现实的，下面使用Python代码方式来实现。

训练数据集如下，其存储在csv格式的delay.csv文件中：

Season	Weather	A_Control	Airline	Delay
Summer	Sunny	no	CZ	no
Summer	Sunny	no	CA	no
Autumn	Sunny	no	CZ	yes
WinterSpring	RainyOrSnowy	no	CZ	yes
WinterSpring	Cloudy	yes	CZ	yes
WinterSpring	Cloudy	yes	CA	no
Autumn	Cloudy	yes	CA	yes
Summer	RainyOrSnowy	no	CZ	no
Summer	Cloudy	yes	CZ	yes
WinterSpring	RainyOrSnowy	yes	CZ	yes
Summer	RainyOrSnowy	yes	CA	yes
Autumn	RainyOrSnowy	no	CA	yes
Autumn	Sunny	yes	CZ	yes
WinterSpring	RainyOrSnowy	no	CA	no
Summer	RainyOrSnowy	no	CA

最后一行数据是需要预测的数据。

sklearn库的naive_bayes模块提供了CategoricalNB朴素贝叶斯分类器来处理特征为类别型数据的分类问题，代码如下：

#!/usr/bin/env python
# -*- encoding: utf-8 -*-
'''
@File    :   Bayes.py    
@Contact :   raogx.vip@hotmail.com
@License :   (C)Copyright 2017-2018, Liugroup-NLPR-CASIA

@Modify Time      @Author    @Version    @Desciption
------------      -------    --------    -----------
2024/8/24 12:19   gxrao      1.0         None
'''

import pandas as pd
from sklearn.naive_bayes import CategoricalNB
from sklearn.preprocessing import LabelEncoder

if __name__ == '__main__':
    # 读取数据.
    df = pd.read_csv(filepath_or_buffer='G:\pycharm_workspace\machine_learning\cda\data\delay.csv',
                     sep="\t",
                     header=0)
    print('原始数据df=', '\n', df)

    print("==========================================================")
    # 对数据进行编码，字符串特征转换成数字型特征
    label_encoder = LabelEncoder()
    label_df = pd.DataFrame()
    for col in df.columns:
        label_df[col] = label_encoder.fit_transform(df[col])
    print('转换后的数据label_df=', '\n', label_df)

    print("==========================================================")
    # 定义离散型的朴素贝叶斯分类器
    bayes_classifier = CategoricalNB()

    # 数据划分为训练集和测试集.
    train_df = label_df[:-1]
    test_df = label_df[-1:]

    X_train = train_df.drop('Delay', axis=1)
    print('训练数据特征列X_train=', '\n', X_train)
    y_train = train_df['Delay']
    print('训练数据标签列y_train=', '\n', y_train)

    # 使用训练数据拟合分类器
    bayes_classifier.fit(X_train, y_train)

    print("==========================================================")
    # 预测测试数据.
    X_test = test_df.drop('Delay', axis=1)
    predictions = bayes_classifier.predict(X_test)

    # 将预测结果数值转换成原始字符串.
    prediction_ray = label_encoder.inverse_transform(predictions)
    print('预测结果predictions', '\n', prediction_ray)

（1）pd.read_csv()从csv格式的文件中读取样本数据

（2）LabelEncoder()用于将类别字符串转换成数字，因为Python无法基于字符串直接计算，需要将其转换成数字。可以理解成对类别字符串进行编码。

（3）label_encoder.inverse_transform()将预测的数字型结果再转回编码前的字符串形式。

程序运行结果如下：

原始数据df= 
           Season       Weather A_Control Airline Delay
0         Summer         Sunny        no      CZ    no
1         Summer         Sunny        no      CA    no
2         Autumn         Sunny        no      CZ   yes
3   WinterSpring  RainyOrSnowy        no      CZ   yes
4   WinterSpring        Cloudy       yes      CZ   yes
5   WinterSpring        Cloudy       yes      CA    no
6         Autumn        Cloudy       yes      CA   yes
7         Summer  RainyOrSnowy        no      CZ    no
8         Summer        Cloudy       yes      CZ   yes
9   WinterSpring  RainyOrSnowy       yes      CZ   yes
10        Summer  RainyOrSnowy       yes      CA   yes
11        Autumn  RainyOrSnowy        no      CA   yes
12        Autumn         Sunny       yes      CZ   yes
13  WinterSpring  RainyOrSnowy        no      CA    no
14        Summer  RainyOrSnowy        no      CA   NaN
==========================================================
转换后的数据label_df= 
     Season  Weather  A_Control  Airline  Delay
0        1        2          0        1      0
1        1        2          0        0      0
2        0        2          0        1      1
3        2        1          0        1      1
4        2        0          1        1      1
5        2        0          1        0      0
6        0        0          1        0      1
7        1        1          0        1      0
8        1        0          1        1      1
9        2        1          1        1      1
10       1        1          1        0      1
11       0        1          0        0      1
12       0        2          1        1      1
13       2        1          0        0      0
14       1        1          0        0      2
==========================================================
训练数据特征列X_train= 
     Season  Weather  A_Control  Airline
0        1        2          0        1
1        1        2          0        0
2        0        2          0        1
3        2        1          0        1
4        2        0          1        1
5        2        0          1        0
6        0        0          1        0
7        1        1          0        1
8        1        0          1        1
9        2        1          1        1
10       1        1          1        0
11       0        1          0        0
12       0        2          1        1
13       2        1          0        0
训练数据标签列y_train= 
 0     0
1     0
2     1
3     1
4     1
5     0
6     1
7     0
8     1
9     1
10    1
11    1
12    1
13    0
Name: Delay, dtype: int32
==========================================================
预测结果predictions 
 ['no']

Process finished with exit code 0

结果最终预测Delay为no，与手工计算结果一致。

朴素贝叶斯扩展

naive_bayes模块中，朴素贝叶斯分类器根据特征服从的分布不同有不同的分类器实现，常见的有

（1）高斯贝叶斯分类器（GaussianNB），

（2）伯努利贝叶斯分类器（BernoulliNB），

（3）多项贝叶斯分类器（MultinomialNB），

（4）类别贝叶斯分类器（CategoricalNB）。

如图：

类之间的继承关系如图：

下面逐个介绍它们的原理以及Python源码实现。

高斯贝叶斯分类器

高斯贝叶斯分类器就是假设特征服从的分布是高斯分布。

原理

现在我们回到贝叶斯公式的估计上来，假设已经知道了在分类为 c_k 的情况下，第个特征 $x_i$ 服从如下的正态分布：

$p(x_i|c_k)=\frac{1}{\sqrt{2\pi}\sigma_{k}}e^{-\frac{(x_i-u_k)^2}{2\sigma_{k}^2}}$

这里 $u_k,\sigma_k^2$ 就是我们前面极大似然估计里面说的未知参数 $\theta$ ，只不过这里是两个未知参数，但是道理是一样的。

事实上，对于每个分类 c_k ，都有一个 $u_k,\sigma_k^2$ 与之对应。从表中找到类别为 c_k 的所有数据，然后将它们的第列特征的值提取出来组成一个集合，记为 $A_k=\{t_1,t_2,...,t_q\}$ 。事实上 A_k 就是在分类为 c_k ，第个特征 $x_i$ 服从的正态分布的情况下产生的数据。使 $P(A_k|c_k)$ 达到最大的 $\hat{u_k},\hat{\sigma_k^2}$ 就是 $u_k,\sigma_k^2$ 的极大似然估计。

下面构造极大似然函数具体推导，因为各数据之间相互独立，所以

$P(A_k|c_k)=\prod_{i=1}^{q}P(t_i|c_k)$ ，

两边取对数，得到

$lnP(A_k|c_k)=\sum_{i=1}^{q}lnP(t_i|c_k)\\ =\sum_{i=1}^{q}ln\frac{1}{\sqrt{2\pi}\sigma_k}e^{-\frac{(t_i-u_k)^2}{2\sigma_k^2}} =-\sum_{i=1}^{q}[\frac{1}{2}ln\sigma_k^2+\frac{1}{2}ln2\pi+\frac{(t_i-u_k)^2}{2\sigma_k^2}]$

两边对 $u_k,\sigma_k^2$ 求偏导数，得到：

$\frac{\partial lnP(A_k|c_k)}{\partial u_k}=\sum_{i=1}^{q}\frac{t_i-u_k}{\sigma_k^2}$

$\frac{\partial lnP(A_k|c_k)}{\partial \sigma_k^2}=-\sum_{i=1}^{q}[\frac{1}{2\sigma_k^2}-\frac{(t_i-u_k)^2}{2(\sigma_k^2)^2}]$

令偏导数为0，得到 $u_k,\sigma_k^2$ 的极大似然估计：

$u_k=\frac{1}{q}\sum_{i=1}^{q}t_i$

$\sigma_k^2=\frac{1}{q}\sum_{i=1}^{q}(t_i-u_k)^2$

从结果来看， $u_k,\sigma_k^2$ 分别是集合 A_k 中所有数据的平均值和方差。

以上估计出了 $P(x_i|c_k)$ ，同理使用极大似然估计很容易估计出 $P(\mathbf{x})$ ，而 $P(c_k)$ 可以使用离散估计。

应用

下面使用Python中实现的高斯贝叶斯分类器来训练数据，并对数据进行分类，具体代码如下：

from sklearn import datasets
from sklearn.naive_bayes import GaussianNB

iris = datasets.load_iris()
gnb = GaussianNB()
y_pred = gnb.fit(iris.data, iris.target).predict(iris.data)
print("Number of mislabeled points out of a total %d points : %d"
      % (iris.data.shape[0],(iris.target != y_pred).sum()))

代码主要分为以下几个步骤：

（1）加载scikit-learn中自带了标准数据集iris，对应代码如下：

iris = datasets.load_iris()

scikit-learn中自带了一些标准数据集（datasets）, 例如在分类中有iris和digits，在回归中有boston house prices dataset。我们这里使用的就是iris数据集，关于iris数据集，Wikipedia中介绍如下：

但是截至目前，iris数据集中已经变成150个样本，每个样本有4个特征，如下：

from sklearn import datasets

iris = datasets.load_iris()
print(iris.data.shape)

程序输出结果：

(150, 4)

我们可以查看iris数据集的data和target。如下：

print('前5条data:',iris.data[0:5])
print('target:',iris.target)

前5条data: [[5.1 3.5 1.4 0.2]
 [4.9 3.  1.4 0.2]
 [4.7 3.2 1.3 0.2]
 [4.6 3.1 1.5 0.2]
 [5.  3.6 1.4 0.2]]
target: [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2
 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
 2 2]

（2）初始化高斯贝叶斯分类器GaussianNB，这里使用的是默认初始化，对应代码如下：

gnb = GaussianNB()

（3）使用标准数据集iris训练出一个高斯贝叶斯分类器，其实这一步使用极大似然估计完成了正态分布中的未知参数估计，对应代码如下：

gnb.fit(iris.data, iris.target)

（4）使用训练出来的高斯贝叶斯分类器对数据进行预测，这里预测的是标准数据集iris，对应代码如下：

y_pred = gnb.fit(iris.data, iris.target).predict(iris.data)

程序运行结果：

Number of mislabeled points out of a total 150 points : 6

结果表明：使用iris数据集训练得到的高斯贝叶斯分类器，只对6个样本分类错误，准确度可以接受。

源码分析

首先来看三个贝叶斯分类器的基类BaseNB，具体代码如下：

class BaseNB(six.with_metaclass(ABCMeta, BaseEstimator, ClassifierMixin)):
    """贝叶斯估计（naive Bayes estimators）的抽象类"""

    @abstractmethod
    def _joint_log_likelihood(self, X):

    def predict(self, X):
        jll = self._joint_log_likelihood(X)
        return self.classes_[np.argmax(jll, axis=1)]

    def predict_log_proba(self, X):
        jll = self._joint_log_likelihood(X)
        # normalize by P(x) = P(f_1, ..., f_n)
        log_prob_x = logsumexp(jll, axis=1)
        return jll - np.atleast_2d(log_prob_x).T

    def predict_proba(self, X):
        return np.exp(self.predict_log_proba(X))

该类包含了一个抽象方法和三个已经实现的方法：

（1）_joint_log_likelihood(self, X)

该方法是是一个抽象方法，具体实现由子类完成。它用于计算联合概率 $P(c_k,\mathbf{x})$ 的对数，也就是 $lnP(c_k,\mathbf{x})=lnP(c_k)+lnP(\mathbf{x}|c_k)$ 。返回如下 $m\times p$ 阶联合概率对数矩阵：

$jll=\begin{bmatrix} lnP(c_1,\mathbf{x_1}) & lnP(c_2,\mathbf{x_1}) & ... &lnP(c_p,\mathbf{x_1}) \\ lnP(c_1,\mathbf{x_2}) & lnP(c_2,\mathbf{x_2}) & ... &lnP(c_p,\mathbf{x_2}) \\ \vdots&\vdots & &\vdots \\ lnP(c_1,\mathbf{x_m}) & lnP(c_2,\mathbf{x_m}) & ... &lnP(c_p,\mathbf{x_m}) \\ \end{bmatrix}$

（2）predict(self, X)

该方法通过调用_joint_log_likelihood(self, X)获得联合概率对数矩阵 $jll$ ，并选择矩阵每行中最大值对应的下标，该下标也是类的下标。

（3）predict_log_proba(self, X)

该方法就是标准化联合概率对数矩阵 $jll$ 。具体先计算第行：

$total=ln\sum_{j=1}^{p}exp(lnP(c_j,\mathbf{x}_i))$ ，

然后第行的每个元素除以 $total$ ，这样就得到一个标准化后的联合概率对数矩阵。

（4）predict_proba(self, X)

该方法是通过消除标准化后的联合概率对数矩阵中的对数，计算真正的概率。

下面来看它的一个实现类GaussianNB的源码，使用GaussianNB分类器分类，主要分为5个阶段：

（1）数据加载

查看上面的应用。

（2）分类器初始化

当数据加载之后，接下来是初始化GaussianNB分类器，实际上调用的是__init__()方法，源码如下：

    def __init__(self, priors=None):
        self.priors = priors

该方法中初始化了priors，其实是 P(c_k) 的先验概率，可以由程序员事先给定，如果选择默认实现，就使用离散估计中的 P(c_k) 。

（3）训练或拟合

当GaussianNB分类器初始化结束后，调用了fit()方法，具体代码如下：

    def fit(self, X, y, sample_weight=None):
        X, y = check_X_y(X, y)
        return self._partial_fit(X, y, np.unique(y), _refit=True,
                                 sample_weight=sample_weight)

该方法很简单，首先检查数据X和分类y是否有效，然后调用_partial_fit()训练数据，注意参数_refit=True，表示如果之前已经拟合过，全部抛弃，重新开始拟合。事实上，_partial_fit()方法才是GaussianNB分类器真正核心，具体代码如下：

    def _partial_fit(self, X, y, classes=None, _refit=False,
                     sample_weight=None):
        #检查数据有效性
        X, y = check_X_y(X, y)
        
        if sample_weight is not None:
            sample_weight = check_array(sample_weight, ensure_2d=False)
            check_consistent_length(y, sample_weight)

        epsilon = 1e-9 * np.var(X, axis=0).max()
        #如果重新拟合，将类别设置为None
        if _refit:
            self.classes_ = None
        #检查是否第一次拟合
        if _check_partial_fit_first_call(self, classes):
            #初始化特征数，类别数，u_k和sigma_k
            n_features = X.shape[1]
            n_classes = len(self.classes_)
            self.theta_ = np.zeros((n_classes, n_features))
            self.sigma_ = np.zeros((n_classes, n_features))

            self.class_count_ = np.zeros(n_classes, dtype=np.float64)

            # Initialise the class prior
            n_classes = len(self.classes_)
            # 当P(c_k)不为空时，对其进行有效性检验
            if self.priors is not None:
                priors = np.asarray(self.priors)
                # Check that the provide prior match the number of classes
                if len(priors) != n_classes:
                    raise ValueError('Number of priors must match number of'
                                     ' classes.')
                # Check that the sum is 1
                if priors.sum() != 1.0:
                    raise ValueError('The sum of the priors should be 1.')
                # Check that the prior are non-negative
                if (priors < 0).any():
                    raise ValueError('Priors must be non-negative.')
                self.class_prior_ = priors
            #第一次拟合时，如果P(C_K)为None，全部初始化为0
            else:
                # Initialize the priors to zeros for each class
                self.class_prior_ = np.zeros(len(self.classes_),
                                             dtype=np.float64)
        else:
            if X.shape[1] != self.theta_.shape[1]:
                msg = "Number of features %d does not match previous data %d."
                raise ValueError(msg % (X.shape[1], self.theta_.shape[1]))
            # Put epsilon back in each time
            self.sigma_[:, :] -= epsilon

        classes = self.classes_

        unique_y = np.unique(y)
        unique_y_in_classes = np.in1d(unique_y, classes)

        if not np.all(unique_y_in_classes):
            raise ValueError("The target label(s) %s in y do not exist in the "
                             "initial classes %s" %
                             (unique_y[~unique_y_in_classes], classes))
        #对于每个类c_k，估计u_k和sigma_k
        for y_i in unique_y:
            i = classes.searchsorted(y_i)
            X_i = X[y == y_i, :]

            if sample_weight is not None:
                sw_i = sample_weight[y == y_i]
                N_i = sw_i.sum()
            else:
                sw_i = None
                N_i = X_i.shape[0]

            new_theta, new_sigma = self._update_mean_variance(
                self.class_count_[i], self.theta_[i, :], self.sigma_[i, :],
                X_i, sw_i)

            self.theta_[i, :] = new_theta
            self.sigma_[i, :] = new_sigma
            self.class_count_[i] += N_i

        self.sigma_[:, :] += epsilon

        # p(c_k)的计算，它可通过实例化GaussianNB分类器时参数传入
        # 如果p(c_k)为None，使用频率代替
        if self.priors is None:
            # Empirical prior, with sample_weight taken into account
            self.class_prior_ = self.class_count_ / self.class_count_.sum()

        return self

已通过注释方式对代码的核心部分进行了解释。对照高斯贝叶斯分类器的极大似然估计推导，不难理解。

（4）预测

当拟合（训练）好高斯贝叶斯分类器后，就可以调用predict()方法使用它对数据进行预测分类了，实际上，GaussianNB类并未对predict()方法进行重写，而是直接调用基类BaseNB中已经实现的predict()方法，接着predict()调用GaussianNB实现的_joint_log_likelihood()方法，具体代码如下：

    def _joint_log_likelihood(self, X):
        check_is_fitted(self, "classes_")

        X = check_array(X)
        joint_log_likelihood = []
        for i in range(np.size(self.classes_)):
            jointi = np.log(self.class_prior_[i])
            n_ij = - 0.5 * np.sum(np.log(2. * np.pi * self.sigma_[i, :]))
            n_ij -= 0.5 * np.sum(((X - self.theta_[i, :]) ** 2) /
                                 (self.sigma_[i, :]), 1)
            joint_log_likelihood.append(jointi + n_ij)

        joint_log_likelihood = np.array(joint_log_likelihood).T
        return joint_log_likelihood

该方法具体完成 $jll$ 矩阵计算。当predict()方法拿到 $jll$ 矩阵后就可以完成预测分类了。

（5）评估
查看上面应用。

伯努利贝叶斯分类器

原理

首先看一个经典分类问题：文本分类。我们将一切可能的词组成一个字典。当然实际中“一切可能”的词都放入字典中是不可能的，但是我们先这么假设。我们如下来构造文本的一个特征向量：

首先按照字典中词的个数构造一个维向量，通常字典是比较大的，所以对应的也非常大。当字典的第个词在文本中出现过，维向量的第个分量就取值为1，否则取值为0。这样就得到一个取值为0或1的维向量，它就是文本对应的特征向量。

我们的任务是使用文本的特征向量对文本进行分类。每一个文本特征向量就是一个数据，只是每个数据的特征只能取值为0或1。这就是伯努利贝叶斯分类器（Bernoulli Naive Bayes）的模型。它是一种实现服从多元伯努利分布数据的朴素贝叶斯训练和分类算法。下面使用贝叶斯公式，对伯努利贝叶斯分类器模型进行推导。推导过程与上面的高斯朴素分类器一样，只不过这里的数据服从的不再是正态分布，而是伯努利分布。对于元随机向量 $(x_1,x_2,...,x_n)$ ，分量 $x_i$ 是一个随机变量，在分类为的情况下，服从如下伯努利分布：

$P(x_i|y)=P(x_i=1|y)x_i+P(x_i=0|y)(1-x_i)$

注意到 $P(x_i=0|y)=1-P(x_i=1|y)$ ，所以

$P(x_i|y)=P(x_i=1|y)x_i+(1-P(x_i=1|y))(1-x_i)$

P(x_i=1|y) 就是需要估计的参数，不失一般性，我们估计在分类为 c_k 和 $x_i$ 下的参数，简记为 $P_{ik}$ 。简化上式得到

$P(x_i|y)=P_{ik}[x_i-(1-x_i)]+(1-x_i)$ ，

后面代码实现也是采用的这种形式。

从表中找到类为 c_k 的所有数据，然后将它们的第列特征数据提取出来组成一个集合，记为 $A_k=\{t_1,t_2,...,t_q\}$ ，并设 A_k 中等于 $1$ 的数有 $n_k$ 个。同样使用极大似然估计，构造似然函数

$P(A_k|c_k)=\prod_{i=1}^{q}P(t_i|c_k)$

两边取对数，得到

$lnP(A_k|c_k)=\sum_{j=1}^{q}lnP(t_j|c_k)=\sum_{j=1}^{q}ln[P_{ik}t_j+(1-P_{ik})(1-t_j)]$

两边对 $P_{ik}$ 求导数，得到：

$\frac{\partial lnP(A_k|c_k)}{\partial P_{ik}}=\sum_{j=1}^{q}\frac{t_j-(1-t_j)}{P_{ik}t_j+(1-P_{ik})(1-t_j)}=\frac{n_k}{P_{ik}}-\frac{q-n_k}{1-P_{ik}}$

令导数为0，得到 $P_{ik}$ 的极大似然估计：

$P_{ik}=\frac{n_k}{q}$

从结果来看， $P_{ik}$ 是类别为 c_k 的数据中第个分量为1的数据的占比。

从结果来看， $P_{ik}$ 的极大似然估计似乎很完美，但是如果获取的数据中的第个分量全部为0，则使用这些数据估计出的概率 $P_{ik}=0$ 。这是不符合实际情况的，虽然获得的数据中没有出现第个分量为1的数据，但不代表第个分量为1的数据不会在以后预测中出现。为了规避这个问题，引入一个平滑系数 $\alpha$ ，如下：

$P_{ik}=\frac{n_k+\alpha}{q+\alpha}$

其中 $0<\alpha \leq 1$ 。

源码分析

从类继承关系图来看，伯努利贝叶斯分类器BernoulliNB，并没有直接继承BaseNB，而是继承BaseDiscreteNB。该类是离散贝叶斯分类器的具体抽象，包含BernoulliNB和MultinomialNB两个子类。其中MultinomialNB就是后面要详细介绍的多项式贝叶斯分类器。依然按照数据加载、分类器初始化、训练（拟合）、预测和评估五个阶段来分析BernoulliNB的源码：

（1）数据加载

略。

（2）分类器初始化

    def __init__(self, alpha=1.0, binarize=.0, fit_prior=True,
                 class_prior=None):
        self.alpha = alpha
        self.binarize = binarize
        self.fit_prior = fit_prior
        self.class_prior = class_prior

self.alpha：初始化平滑系数；

self.binarize：二值化的阈值，BernoulliNB分类器处理的只有0和1值的情况，如果输入是连续值，可以使用二值化阈值将连续值转换成0和1的值；

self.fit_prior：是否学习类的先验概率，默认为true，这样会在训练数据不断增加的同时来优化类的先验概率。

self.class_prior：类的先验概率，可以事先指定类的先验概率。如果未指定，使用类的频率代替。

（3）训练（拟合）

BernoulliNB中并没有fit()方法，而是调用的基类BaseDiscreteNB中的fit()方法，源码如下：

    def fit(self, X, y, sample_weight=None):
        #数据检查
        X, y = check_X_y(X, y, 'csr')
        #获取数据的特征数
        _, n_features = X.shape

        #标签处理，获取矩阵Y，Y的行表示数据，列表示类别，取值只有0和1
        labelbin = LabelBinarizer()
        Y = labelbin.fit_transform(y)
        self.classes_ = labelbin.classes_
        if Y.shape[1] == 1:
            Y = np.concatenate((1 - Y, Y), axis=1)

        Y = Y.astype(np.float64)
        if sample_weight is not None:
            sample_weight = np.atleast_2d(sample_weight)
            Y *= check_array(sample_weight).T

        class_prior = self.class_prior

        n_effective_classes = Y.shape[1]
        #初始化类别和特征频率
        self.class_count_ = np.zeros(n_effective_classes, dtype=np.float64)
        self.feature_count_ = np.zeros((n_effective_classes, n_features),
                                       dtype=np.float64)
        #调用实现类的_count/_update_feature_log_prob/_update_class_log_prior方法
        self._count(X, Y)
        alpha = self._check_alpha()
        self._update_feature_log_prob(alpha)
        self._update_class_log_prior(class_prior=class_prior)
        return self

该方法比较容易理解，已通过注释方式将核心部分进行了解析。fit()方法核心是调用了子类BernoulliNB的三个方法：_count()、_update_feature_log_prob()和_update_class_log_prior()，具体来看：

    def _count(self, X, Y):
        """Count and smooth feature occurrences."""
        if self.binarize is not None:
            X = binarize(X, threshold=self.binarize)
        self.feature_count_ += safe_sparse_dot(Y.T, X)
        self.class_count_ += Y.sum(axis=0)

_count()方法计算特征和类的频率。

   def _update_feature_log_prob(self, alpha):
        """Apply smoothing to raw counts and recompute log probabilities"""
        smoothed_fc = self.feature_count_ + alpha
        smoothed_cc = self.class_count_ + alpha * 2

        self.feature_log_prob_ = (np.log(smoothed_fc) -
                                  np.log(smoothed_cc.reshape(-1, 1)))

_update_feature_log_prob()方法计算加入平滑系数后的 $lnP_{ik}$ 。

    def _joint_log_likelihood(self, X):
        """Calculate the posterior log probability of the samples X"""
        check_is_fitted(self, "classes_")

        X = check_array(X, accept_sparse='csr')

        if self.binarize is not None:
            X = binarize(X, threshold=self.binarize)

        n_classes, n_features = self.feature_log_prob_.shape
        n_samples, n_features_X = X.shape

        if n_features_X != n_features:
            raise ValueError("Expected input with %d features, got %d instead"
                             % (n_features, n_features_X))

        neg_prob = np.log(1 - np.exp(self.feature_log_prob_))
        # Compute  neg_prob · (1 - X).T  as  ∑neg_prob - X · neg_prob
        jll = safe_sparse_dot(X, (self.feature_log_prob_ - neg_prob).T)
        jll += self.class_log_prior_ + neg_prob.sum(axis=1)

        return jll

_joint_log_likelihood方法计算联合概率对数矩阵 $jll$ 。

（4）预测

略。

（5）评估

略。

整个计算流程如图：

多项朴素贝叶斯分类器

多项朴素贝叶斯分类器（Multinomial Naive Bayes）是一种实现服从多项分布数据的分类算法。多项分布其实是伯努利分布的直接推广。在多项分布中，出现的结果不再是两个（0或1），而是 $k$ （ $k\geq 2$ ）个。所以，伯努利分布是多项分布中 $k=2$ 的情况。假设在分类为 c_k 的情况下，第个特征 $x_i$ 服从多项分布，多项分布出现的可能结果为 $\{a_1,a_2,...,a_k\}$ ， a_i 取值0或1， $\sum_{i=1}^{k}a_i=1$ ，且 $P(a_i=1)=p_i$ ， $\sum_{i=1}^{k}p_i=1$ 。

从表中找到类为 c_k 且第个特征取值为 $a_j$ 的所有数据，然后将它们的第列特征数据提取出来组成一个集合，记为 $A_k=\{t_1,t_2,...,t_q\}$ ，并设 A_k 中等于 $1$ 的数有 $n_k$ 个。同样使用极大似然估计，构造似然函数

$P(A_k|c_k)=\prod_{i=1}^{q}P(t_i|c_k)$

两边取对数，得到

$lnP(A_k|c_k)=\sum_{i=1}^{q}lnP(t_j|c_k)=\sum_{i=1}^{q}ln[p_ja_j+(1-p_j)(1-a_j)]$

两边对 $p_j$ 求导数，得到：

$\frac{\partial lnP(A_k|c_k)}{\partial p_j}=\sum_{i=1}^{q}\frac{a_j-(1-a_j)}{p_ja_j+(1-p_j)(1-a_j)}=\frac{n_k}{p_j}-\frac{q-n_k}{1-p_j}$

令导数为0，得到 $p_j$ 的极大似然估计：

$p_j=\frac{n_k}{q}$

这个结果和离散估计是一致的。

与伯努利分布类似，引入一个平滑系数 $\alpha$ ，如下：

$p_j=\frac{n_k+\alpha }{q+n\alpha }$

源码与伯努利朴素贝叶斯分类器类似，不再累述。

我们回到文本分类的问题上来，在伯努利朴素分类器中，提取了文本特征向量，每个向量取值为0或1。但是在文本分类和主题分析时，这种方法不是很好。它只度量了字典中的词是否在文本中出现，对于出现的频率并没有度量。现在换一种方法提取文本特征。

构建一个维向量，向量的第个分量表示字典中第个位置的词在文本中出现的频率。

此时，每个分量的取值是非负整数。这个就是文本的多项分布模型。文字对主题文本的影响分为两部分：

（1）某个文字在文本中出现的比例越高，它与文本主题也就越相关；

（2）一个文字出现的文本数占总文本数的比例越小，它与文本主题也就越相关；

第一条就是熟悉的TF，第二条是IDF，假设文本的特征向量为 $\mathbf{x}_i=(x_{i1},x_{i2},...,x_{in})$ ，则

$TF_{ij}=\frac{x_{ij}}{\sum_{j=1}^nx_{ij}}$

$IDF_j=ln\frac{m}{\sum_{i=1}^m1_{x_{ij}>0}}$

$TFIDF_{ij}=TF_{ij}*IDF_{j}$

半朴素贝叶斯分类器

在朴素贝叶斯分类器中，我们严格假设特征之间完全独立的。在实际情况中，几乎是不可能的。这就需要考虑特征之间的依赖关系。例如常见的半朴素贝叶斯分类器有SPODE、TAN、AODE。

TAN分类器是由Friedman 等人提出的一种树状贝叶斯网络, 是朴素贝叶斯分类器的一种改进模型, TAN 分类器的分类性能明显优于朴素贝叶斯分类器，其基本思路是放松朴素贝叶斯分类器中的独立性假设条件, 借鉴贝叶斯网络中表示依赖关系的方法, 扩展朴素贝叶斯的结构, 使其能容纳属性间存在的依赖关系, 但对其表示依赖关系的能力加以限制。

半朴素贝叶斯分类器的原理，这里就不再详述。