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

传送门

《统计学习方法》读书笔记——机器学习常用评价指标

《统计学习方法》读书笔记——感知机（原理+代码实现）

《统计学习方法》读书笔记——K近邻法（原理+代码实现）

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

朴素贝叶斯法

传送门写在前面朴素贝叶斯法代码实现参考

写在前面

朴素贝叶斯法与贝叶斯估计是不同的概念。 损失函数与风险函数损失函数用于度量一次预测的好坏； 风险函数用于度量平均意义下模型的好坏。 全概率公式与逆概率公式设A1,A2,...,AnA_1,A_2,...,A_nA1​,A2​,...,An​为一组完备事件组，则对任一事件BBB，有如下全概率公式：

P(B)=∑i=1nP(Ai)P(B∣Ai)P(B) =\sum_{i=1}^nP(A_i)P(B|A_i) P(B)=i=1∑n​P(Ai​)P(B∣Ai​)

若P(B)>0P(B)>0P(B)>0，则有如下贝叶斯公式，或称逆概率公式：

P(Ai∣B)=P(AiB)P(B)=P(Ai)P(B∣Ai)∑j=1nP(Aj)P(B∣Aj)P(A_i|B)=\frac{P(A_iB)}{P(B)}=\frac{P(A_i)P(B|A_i)}{\sum\limits_{j=1}^nP(A_j)P(B|A_j)} P(Ai​∣B)=P(B)P(Ai​B)​=j=1∑n​P(Aj​)P(B∣Aj​)P(Ai​)P(B∣Ai​)​ 先验概率与后验概率常称P(Ai)P(A_i)P(Ai​)为事件AiA_iAi​发生的先验概率，而称P(Ai∣B)P(A_i|B)P(Ai​∣B)为事件AiA_iAi​发生的后验概率。(概率论教程P25)

朴素贝叶斯法

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

要搞明白朴素贝叶斯法的原理，要先知道“球门方向”在哪里。

设输入空间X⊆Rn\mathcal{X}\subseteq{\bm{R}^n}X⊆Rn为nnn维向量的集合，输入空间为类标记集合Y={c1,c2,...,cK}\mathcal{Y} = \{c_1,c_2,...,c_K\}Y={c1​,c2​,...,cK​}。输入为特征向量x∈Xx\in\mathcal{X}x∈X，输入为类别y∈Yy\in\mathcal{Y}y∈Y。XXX是定义在输入空间上的随机变量，YYY是定义在输出空间Y\mathcal{Y}Y上的随机变量。对于一个输入的数据xxx，要预测xxx所属的类别，相当于求以下概率的值：

P(Y=ck∣X=x)P(Y=c_k|X=x) P(Y=ck​∣X=x)

根据贝叶斯公式（逆概率公式），可得：

P(Y=ck∣X=x)=P(X=x∣Y=ck)P(Y=ck)∑kP(X=x∣Y=ck)P(Y=ck)P(Y=c_k|X=x) = \frac{P(X=x|Y=c_k)P(Y=c_k)}{\sum\limits_k P(X=x|Y=c_k)P(Y=c_k)} P(Y=ck​∣X=x)=k∑​P(X=x∣Y=ck​)P(Y=ck​)P(X=x∣Y=ck​)P(Y=ck​)​

插入部分

上述的贝叶斯公式中，P(Y=ck)P(Y=c_k)P(Y=ck​)容易求得，统计数据集中各个类别样本的数量就可以计算出来。

因为xxx是向量，即x=(x(1),x(2),...,x(n))x = (x^{(1)},x^{(2)},...,x^{(n)})x=(x(1),x(2),...,x(n))，因此将P(X=x∣Y=ck)P(X=x|Y=c_k)P(X=x∣Y=ck​)部分展开来就是P(X(1)=x(1),X(2)=x(2),...,X(n)=x(n)∣Y=ck)P(X^{(1)}=x^{(1)},X^{(2)}=x^{(2)},...,X^{(n)}=x^{(n)} |Y=c_k)P(X(1)=x(1),X(2)=x(2),...,X(n)=x(n)∣Y=ck​)。

这样看来，假设对mnist数据集来说，每一张图片有K=10K=10K=10种可能的预测结果（即ck,k=1,2,...,10c_k,k=1,2,...,10ck​,k=1,2,...,10），每一张图片由28*28=784个像素点组成(即x(j)x^{(j)}x(j),j=1,2,...,784j=1,2,...,784j=1,2,...,784)，每一个像素点的取值范围为0<c(j)<2550<c^{(j)}<2550<c(j)<255，即x(j)x^{(j)}x(j)的可值为Sj=255S_j=255Sj​=255个。那么要计算的数据量为K∏j=1n=10×255784K\prod\limits_{j=1}^n=10\times255^{784}Kj=1∏n​=10×255784，无法计算。

朴素贝叶斯法对条件概率分布作了条件独立性的假设，朴素贝叶斯法也由此得名，即

P(X=x∣Y=ck)=P(X(1)=x(1),X(2)=x(2),...,X(n)=x(n)∣Y=ck)=P(X(1)=x(1)∣Y=ck)P(X(2)=x(2)∣Y=ck)...P(X(n)=x(n)∣Y=ck)=∏j=1nP(X(j)=x(j)∣Y=ck)\begin{aligned} 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\limits_{j=1}^nP(X^{(j)}=x^{(j)}|Y=c_k) \end{aligned} P(X=x∣Y=ck​)​=P(X(1)=x(1),X(2)=x(2),...,X(n)=x(n)∣Y=ck​)=P(X(1)=x(1)∣Y=ck​)P(X(2)=x(2)∣Y=ck​)...P(X(n)=x(n)∣Y=ck​)=j=1∏n​P(X(j)=x(j)∣Y=ck​)​

这一假设使朴素贝叶斯法变得简单，但会牺牲一定的分类准确率。

（因为向量的特征之间大概率是不独立地，如果我们独立了，会无法避免地抛弃一些前后连贯的信息，比方说我说“三人成？”，后面大概率就是个”虎“，这个虎明显依赖于前面的三个字。）

那么上面的贝叶斯公式就可以转换为以下形式：

P(Y=ck∣X=x)=P(X=x∣Y=ck)P(Y=ck)∑kP(X=x∣Y=ck)P(Y=ck)=P(Y=ck)∏j=1nP(X(j)=x(j)∣Y=ck))∑kP(Y=ck)∏j=1nP(X(j)=x(j)∣Y=ck)\begin{aligned} P(Y=c_k|X=x) &= \frac{P(X=x|Y=c_k)P(Y=c_k)}{\sum\limits_k P(X=x|Y=c_k)P(Y=c_k)} \\ &= \frac{P(Y=c_k)\prod\limits_{j=1}^nP(X^{(j)}=x^{(j)}|Y=c_k))}{\sum\limits_k P(Y=c_k)\prod\limits_{j=1}^nP(X^{(j)}=x^{(j)}|Y=c_k)} \end{aligned} P(Y=ck​∣X=x)​=k∑​P(X=x∣Y=ck​)P(Y=ck​)P(X=x∣Y=ck​)P(Y=ck​)​=k∑​P(Y=ck​)j=1∏n​P(X(j)=x(j)∣Y=ck​)P(Y=ck​)j=1∏n​P(X(j)=x(j)∣Y=ck​))​​

这是朴素贝叶斯法分类的基本公式。于是，朴素贝叶斯分类器可表示为

y=arg max⁡ckP(Y=ck)∏j=1nP(X(j)=x(j)∣Y=ck))∑kP(Y=ck)∏j=1nP(X(j)=x(j)∣Y=ck)y= \argmax_{c_k}\frac{P(Y=c_k)\prod\limits_{j=1}^nP(X^{(j)}=x^{(j)}|Y=c_k))}{\sum\limits_k P(Y=c_k)\prod\limits_{j=1}^nP(X^{(j)}=x^{(j)}|Y=c_k)} y=ck​argmax​k∑​P(Y=ck​)j=1∏n​P(X(j)=x(j)∣Y=ck​)P(Y=ck​)j=1∏n​P(X(j)=x(j)∣Y=ck​))​

其中分母对所有ckc_kck​都是相同的，所以

y=arg max⁡ckP(Y=ck)∏j=1nP(X(j)=x(j)∣Y=ck))y= \argmax_{c_k}P(Y=c_k)\prod\limits_{j=1}^nP(X^{(j)}=x^{(j)}|Y=c_k)) y=ck​argmax​P(Y=ck​)j=1∏n​P(X(j)=x(j)∣Y=ck​))

该公式即为朴素贝叶斯分类器。

至于式子里面的两项具体怎么求，我们首先看第一项。

P(Y=ck)=∑i=1NI(yi=ck)N,k=1,2,⋯,KP\left(Y=c_{k}\right)=\frac{\sum_{i=1}^{N} I\left(y_{i}=c_{k}\right)}{N}, \quad k=1,2, \cdots, K P(Y=ck​)=N∑i=1N​I(yi​=ck​)​,k=1,2,⋯,K

N为训练样本的数目，假设我们手里现在有100个样本，那N就是100。

分子中I(...)I(...)I(...)是指示函数，括号内条件为真时指示函数为1，反之为0。分子的意思求在这100个样本里有多少是属于ckc_kck​分类的。

再看第二项

P(X(j)=ajl∣Y=ck)=∑i=1NI(xi(j)=ajl,yi=ck)∑i=1NI(yi=ck)P(X^{(j)}=a_{j l} |Y=c_{k})=\frac{\sum_{i=1}^{N} I(x_{i}^{(j)}=a_{j l},y_{i}=c_{k})}{\sum_{i=1}^{N} I(y_{i}=c_{k})} \\ P(X(j)=ajl​∣Y=ck​)=∑i=1N​I(yi​=ck​)∑i=1N​I(xi(j)​=ajl​,yi​=ck​)​

其中j=1,2,...,n;l=1,2,...,Sj;k=1,2,...,Kj=1,2, ..., n ; \quad l=1,2, ..., S_{j};\quad k=1,2,..., Kj=1,2,...,n;l=1,2,...,Sj​;k=1,2,...,K

构造原理和上面第一项的一样，也是通过指示函数来计数得出概率。那么两项都能求了，给出朴素贝叶斯算法。

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

代码实现

# coding=utf-8import numpy as npimport timedef 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和labelreturn dataArr, labelArrdef NaiveBayes(Py, Px_y, x):featrueNum = 784classNum = 10P = [0] * classNumfor i in range(classNum):sum = 0for 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 = 0for i in range(len(testDataArr)):predict = NaiveBayes(Py, Px_y, testDataArr[i])if predict != testLabelArr[i]:errorCnt += 1return 1 - (errorCnt / len(testDataArr))def getAllProbability(trainDataArr, trainLabelArr):featureNum = 784classNum = 10lambbaVal = 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_yif __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.csvstart to read data：./mnist/mnist_test.csvstart to trainstart to testthe accuracy is: 0.8432999999999999time span: 129.84226727485657

参考

原理：《统计学习方法》

代码： /Dod-o/Statistical-Learning-Method_Code