《统计学习方法》读书笔记——朴素贝叶斯法（公式推导+代码实现）

传送门

《统计学习方法》读书笔记——机器学习常用评价指标
《统计学习方法》读书笔记——感知机（原理+代码实现）
《统计学习方法》读书笔记——K近邻法（原理+代码实现）
《统计学习方法》读书笔记——朴素贝叶斯法（公式推导+代码实现）

---

写在前面

朴素贝叶斯法与贝叶斯估计是不同的概念。

 

损失函数与风险函数

损失函数用于度量一次预测的好坏；

风险函数用于度量平均意义下模型的好坏。

全概率公式与逆概率公式

设 $A_1,A_2,...,A_n$为一组完备事件组，则对任一事件 $B$，有如下 全概率公式：
$$P(B)=\sum _{i=1}^{n}P(A_i)P(B|A_i)$$
若 $P(B)>0$，则有如下 贝叶斯公式，或称 逆概率公式：
$$P(A_i|B)=\frac{P(A_{i}B)}{P(B)}=\frac{P(A_i)P(B|A_i)}{\sum _{j=1}^{n}P(A_j)P(B|A_j)}$$

先验概率与后验概率

常称 $P(A_i)$为事件 $A_i$发生的 先验概率，而称 $P(A_i|B)$为事件 $A_i$发生的 后验概率。(概率论教程P25)

朴素贝叶斯法

“如果你不知道怎样踢球，就往球门方向踢 ” ——施拉普纳

要搞明白朴素贝叶斯法的原理，要先知道“球门方向”在哪里。
设输入空间$\mathcal{X}\subseteq {\mathbf{R}}^n$为$n$维向量的集合，输入空间为类标记集合 Y = { c 1 , c 2 , . . . , c K } \mathcal{Y} = \{c_1,c_2,...,c_K\} Y={c1​,c2​,...,cK​}。输入为特征向量$x\in \mathcal{X}$，输入为类别$y\in \mathcal{Y}$。$X$是定义在输入空间上的随机变量，$Y$是定义在输出空间$\mathcal{Y}$上的随机变量。对于一个输入的数据$x$，要预测$x$所属的类别，相当于求以下概率的值：
$$P(Y=c_k|X=x)$$
根据贝叶斯公式（逆概率公式），可得：
$$P(Y=c_k|X=x)=\frac{P(X=x|Y=c_k)P(Y=c_k)}{\sum _{k}P(X=x|Y=c_k)P(Y=c_k)}$$

---

插入部分
上述的贝叶斯公式中，$P(Y=c_k)$容易求得，统计数据集中各个类别样本的数量就可以计算出来。
因为$x$是向量，即$x=(x^{(1)},x^{(2)},...,x^{(n)})$，因此将$P(X=x|Y=c_k)$部分展开来就是$P(X^{(1)}=x^{(1)},X^{(2)}=x^{(2)},...,X^{(n)}=x^{(n)}|Y=c_k)$。
这样看来，假设对mnist数据集来说，每一张图片有$K=10$种可能的预测结果（即$c_k,k=1,2,...,10$），每一张图片由28*28=784个像素点组成(即$x^{(j)}$,$j=1,2,...,784$)，每一个像素点的取值范围为$0<c^{(j)}<255$，即$x^{(j)}$的可值为$S_j=255$个。那么要计算的数据量为$K\prod _{j=1}^n=10\times 255^{784}$，无法计算。
朴素贝叶斯法对条件概率分布作了条件独立性的假设，朴素贝叶斯法也由此得名，即
$$\begin{array}{rl}P(X=x|Y=c_k)& =P(X^{(1)}=x^{(1)},X^{(2)}=x^{(2)},...,X^{(n)}=x^{(n)}|Y=c_k)\\ & =P(X^{(1)}=x^{(1)}|Y=c_k)P(X^{(2)}=x^{(2)}|Y=c_k)...P(X^{(n)}=x^{(n)}|Y=c_k)\\ & =\prod _{j=1}^{n}P(X^{(j)}=x^{(j)}|Y=c_k)\end{array}$$
这一假设使朴素贝叶斯法变得简单，但会牺牲一定的分类准确率。
（因为向量的特征之间大概率是不独立地，如果我们独立了，会无法避免地抛弃一些前后连贯的信息，比方说我说“三人成？”，后面大概率就是个”虎“，这个虎明显依赖于前面的三个字。）

---

那么上面的贝叶斯公式就可以转换为以下形式：
$$\begin{array}{rl}P(Y=c_k|X=x)& =\frac{P(X=x|Y=c_k)P(Y=c_k)}{\sum _{k}P(X=x|Y=c_k)P(Y=c_k)}\\ & =\frac{P(Y=c_k)\prod _{j=1}^{n}P(X^{(j)}=x^{(j)}|Y=c_k))}{\sum _{k}P(Y=c_k)\prod _{j=1}^{n}P(X^{(j)}=x^{(j)}|Y=c_k)}\end{array}$$
这是朴素贝叶斯法分类的基本公式。于是，朴素贝叶斯分类器可表示为
$$y=\underset{c_k}{\mathrm{argmax}}\frac{P(Y=c_k)\prod _{j=1}^{n}P(X^{(j)}=x^{(j)}|Y=c_k))}{\sum _{k}P(Y=c_k)\prod _{j=1}^{n}P(X^{(j)}=x^{(j)}|Y=c_k)}$$
其中分母对所有$c_k$都是相同的，所以
$$y=\underset{c_k}{\mathrm{argmax}}P(Y=c_k)\prod _{j=1}^{n}P(X^{(j)}=x^{(j)}|Y=c_k))$$
该公式即为朴素贝叶斯分类器。
至于式子里面的两项具体怎么求，我们首先看第一项。
$$P\left(Y=c_k\right)=\frac{{\sum }_{i=1}^{N}I\left(y_i=c_k\right)}{N},k=1,2,\cdots ,K$$
N为训练样本的数目，假设我们手里现在有100个样本，那N就是100。
分子中$I(...)$是指示函数，括号内条件为真时指示函数为1，反之为0。分子的意思求在这100个样本里有多少是属于$c_k$分类的。
再看第二项
$$P(X^{(j)}=a_{jl}|Y=c_k)=\frac{{\sum }_{i=1}^{N}I(x_i^{(j)}=a_{jl},y_i=c_k)}{{\sum }_{i=1}^{N}I(y_i=c_k)}$$
其中$j=1,2,...,n;l=1,2,...,S_j;k=1,2,...,K$
构造原理和上面第一项的一样，也是通过指示函数来计数得出概率。那么两项都能求了，给出朴素贝叶斯算法。

但是，步骤（2）中，那么多概率连乘，如果其中有一个概率为0怎么办？那整个式子直接就是0了，这样不对。所以我们连乘中的每一项都得想办法让它保证不是0。

代码实现

# coding=utf-8
import numpy as np
import time
def loadData(fileName):
    print('start to read data：' + fileName)
    dataArr = []
    labelArr = []
    fr = open(fileName, 'r')
    # 将文件按行读取
    for line in fr.readlines():
        # 对每一行数据按切割符','进行切割，返回字段列表
        curLine = line.strip().split(',')
        # 将mnist图像二值化
        dataArr.append([int(int(num) > 128) for num in curLine[1:]])
        labelArr.append(int(curLine[0]))
    # 返回data和label
    return dataArr, labelArr


def NaiveBayes(Py, Px_y, x):
    featrueNum = 784
    classNum = 10
    P = [0] * classNum
    for i in range(classNum):
        sum = 0
        for j in range(featrueNum):
            sum += Px_y[i][j][x[j]]
        P[i] = sum + Py[i]
    return P.index(max(P))


def model_test(Py, Px_y, testDataArr, testLabelArr):
    errorCnt = 0
    for i in range(len(testDataArr)):
        predict = NaiveBayes(Py, Px_y, testDataArr[i])
        if predict != testLabelArr[i]:
            errorCnt += 1
    return 1 - (errorCnt / len(testDataArr))


def getAllProbability(trainDataArr, trainLabelArr):
    featureNum = 784
    classNum = 10
    lambbaVal = 1

    #初始化先验概率分布存放数组，后续计算得到的P(Y = 0)放在Py[0]中，以此类推
    #数据长度为10行1列
    Py = np.zeros((classNum, 1))
    #对每个类别进行一次循环，分别计算它们的先验概率分布
    #计算公式为书中"4.2节 朴素贝叶斯法的参数估计 公式4.8"
    for i in range(classNum):
        # 计算先验概率
        Py[i] = ((np.sum(np.mat(trainLabelArr) == i)) + lambbaVal) / (len(trainLabelArr) + classNum * lambbaVal)
    # 防止数据过小向下溢出使用log()
    # 在似然函数中通常会使用log的方式进行处理
    Py = np.log(Py)

    # 计算条件概率 Px_y=P（X=x|Y = y）
    # 计算分子部分
    # 2表示一个像素点可取值的个数（因为对数据做了二值化处理）
    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 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()
    # 获取训练集、测试集
    trainData, trainLabel = loadData('./mnist/mnist_train.csv')
    testData, testLabel = loadData('./mnist/mnist_test.csv')

    #开始训练，学习先验概率分布和条件概率分布
    print('start to train')
    Py, Px_y = getAllProbability(trainData, trainLabel)

    #使用习得的先验概率分布和条件概率分布对测试集进行测试
    print('start to test')
    accuracy = model_test(Py, Px_y, testData, testLabel)

    print('the accuracy is:', accuracy)
    print('time span:', time.time() -start)

输出结果

start to read data：./mnist/mnist_train.csv
start to read data：./mnist/mnist_test.csv
start to train
start to test
the accuracy is: 0.8432999999999999
time span: 129.84226727485657

参考

原理：《统计学习方法》
代码： https://github.com/Dod-o/Statistical-Learning-Method_Code