决策论——朴素贝叶斯分类算法的R实现（三）

朴素贝叶斯算法（Naive Bayes, NB) 是应用最为广泛的分类算法之一，它是基于贝叶斯定义和特征条件独立假设的分类器方法。朴素贝叶斯法基于贝叶斯公式计算得到，有着坚实的数学基础，以及稳定的分类效率；NB模型所需估计的参数很少，对缺失数据不太敏感，算法也比较简单，当年的垃圾邮件分类都是基于朴素贝叶斯分类器识别的。朴素贝叶斯名字中的“朴素”二字代表着该算法对概率事件做了很大的简化，简化内容就是各个要素之间是相互独立的。比如今天刮风和气温低，两个要素导致了不下雨的结果。实际上刮风可能导致气温低，而且刮风对于天晴的影响会更大，朴素贝叶斯认为刮风和气温之间相互独立，且对于是否下雨这个结果的影响没有轻重之分。

一、贝叶斯决策理论

贝叶斯决策基于概率和误判损失来选择最优的类别标记。本部分内容将从以下问题导出：对于样本$x$，有$N$种可能的类别标记，即类别空间$y=\{c_1,c_2,\cdots ,c_N\}​$。$\lambda_{ij}$​表示将一个真实标记为$c_j$​的样本误分类为$c_i$所产生的损失。基于后验概率$(c_i|x)$（即对于给定样本，判断样本属于哪个分类，概率最大的那个类别是最可能正确的类别）可获得将样本$x$分类为$c_i$​所产生的期望损失，即在样本$x$上的“条件风险”：

$$R(c_i|x)=\sum_{j=1}^{N}\lambda_{ij}P(c_j|x)$$

我们训练模型的目的就是为了寻找一个映射函数$h:x\rightarrow y$以最小化总体风险：

$$R(h)=E_x[R(h(x)|x)]$$

那么对于每个样本$x$，若$h(x)$能最小化条件风险$(h(x)|x)$,则总体风险$R(h)$也将被最小化。所以为了得到最小的总体风险，只需在每个样本上选择哪个能使条件风险$R(c∣x)$最小的类别标记，即

$$h^*(x)=\underset{c\in y}{arg\ min}R(c|x)$$

其中$h^*$称为贝叶斯最优分类器，与之对应的总体风险$R(h^*)$称为贝叶斯风险。$1-R(h^*)$反映了分类器所能达到的最好性能，即通过机器学习所产生的模型精度的理论上限。
若$R(c|x)$是最小化分类错误率，则误判损失$\lambda_{ij}$可表示为：

$$\lambda_{ij}=\begin{cases} 0, \quad i =j \\ 1 , \quad otherwise \end{cases}$$

此时条件风险表示为： $$R(c|x)=1-P(c|x)$$
于是，最小化分类错误率的贝叶斯最优分类器为：

$$h^*(x)=\underset{c\in y}{argmax}P(c|x)$$

即对每个样本$x$，选择能使后验概率$P(c∣x)$最大的类别标记。

【贝叶斯模型的解释】

以上内容已经对问题解释的很清楚了，如果想要对样本进行分类，我们需要得到最大的$P(c|x)$，但是对于一般的情况下，直接求解，很难求出，所以一般借助贝叶斯定理，从侧面进行求解：

$$P(c|x)=\frac{P(c)P(x|c)}{P(x)}$$

$P(c|x)$叫后验概率，也就是我们要计算的后验概率，知道样本，计算这个样本属于某个类别的概率，概率最大的那个类别是最可能正确的类别。
$P(c)$是类“先验”概率。
$P(x∣c)$是条件概率，也就是在类别$c$的条件下，出现样本$x$的可能性。
对于每个样本$x$，其$P(x)= \frac{x的样本数}{总样本数}$​为恒定的数值(在不考虑特征属性的情况下,如此表示)，所以计算$P(c∣x)$的难点在于求解$P(x∣c)$或称之为“似然函数”。
求解类别$c$的类条件概率$P(x∣c)$（似然函数）有两种方法实现：极大似然估计；朴素贝叶斯分类器。

二、朴素贝叶斯算例

2.1 朴素贝叶斯分类案例分析

案例:采用西瓜数据集3.0进行贝叶斯分类

色泽	根蒂	敲声	纹理	脐部	触感	密度	含糖率	好瓜
青绿	蜷缩	浊响	清晰	凹陷	硬滑	0.697	0.46	是
乌黑	蜷缩	沉闷	清晰	凹陷	硬滑	0.774	0.376	是
乌黑	蜷缩	浊响	清晰	凹陷	硬滑	0.634	0.264	是
青绿	蜷缩	沉闷	清晰	凹陷	硬滑	0.608	0.318	是
浅白	蜷缩	浊响	清晰	凹陷	硬滑	0.556	0.215	是
青绿	稍蜷	浊响	清晰	稍凹	软粘	0.403	0.237	是
乌黑	稍蜷	浊响	稍糊	稍凹	软粘	0.481	0.149	是
乌黑	稍蜷	浊响	清晰	稍凹	硬滑	0.437	0.211	是
乌黑	稍蜷	沉闷	稍糊	稍凹	硬滑	0.666	0.091	否
青绿	硬挺	清脆	清晰	平坦	软粘	0.243	0.267	否
浅白	硬挺	清脆	模糊	平坦	硬滑	0.245	0.057	否
浅白	蜷缩	浊响	模糊	平坦	软粘	0.343	0.099	否
青绿	稍蜷	浊响	稍糊	凹陷	硬滑	0.639	0.161	否
浅白	稍蜷	沉闷	稍糊	凹陷	硬滑	0.657	0.198	否
乌黑	稍蜷	浊响	清晰	稍凹	软粘	0.36	0.37	否
浅白	蜷缩	浊响	模糊	平坦	硬滑	0.593	0.042	否
青绿	蜷缩	沉闷	稍糊	稍凹	硬滑	0.719	0.103	否

其中类先验概率$P(c)$为:

$$P(\text{好瓜=是})=\frac{8}{17}\approx0.471 \quad \quad P(\text{好瓜=否})=\frac{9}{17}\approx0.529$$

然后，每个属性的似然概率（条件概率）$P(x_i|c)$为:

是(好瓜)	否(好瓜)
$P_{青绿|是}=P(色泽=青绿|好瓜=是)=\frac{3}{8}=0.375$	$P_{青绿|否}=P(色泽=青绿|好瓜=否)=\frac{3}{9}=0.333$
$P_{蜷缩|是}=P(根蒂=蜷缩|好瓜=是)=\frac{5}{8}=0.625$	$P_{蜷缩|否}=P(根蒂=蜷缩|好瓜=否)=\frac{3}{9}=0.333$
$P_{浊响|是}=P(敲声=浊响|好瓜=是)=\frac{6}{8}=0.750$	$P_{浊响|否}=P(敲声=浊响|好瓜=否)=\frac{4}{9}=0.444$
$P_{清晰|是}=P(纹理=清晰|好瓜=是)=\frac{7}{8}=0.875$	$P_{清晰|否}=P(纹理=清晰|好瓜=否)=\frac{2}{9}=0.222$
$P_{凹陷|是}=P(脐部=凹陷|好瓜=是)=\frac{6}{8}=0.750$	$P_{凹陷|否}=P(脐部=凹陷|好瓜=否)=\frac{2}{9}=0.222$
$P_{硬滑|是}=P(触感=硬滑|好瓜=是)=\frac{6}{8}=0.750$	$P_{硬滑|否}=P(触感=硬滑|好瓜=否)=\frac{6}{9}=0.667$
$p(密度=0.697|好瓜=是)=\frac{1}{\sqrt{2\pi}0.129}\exp\left(-\frac{(0.697-0.574)^2}{2\cdot 0.129^2}\right)\approx1.959$	$p(密度=0.697|好瓜=否)=\frac{1}{\sqrt{2\pi}0.195}\exp\left(-\frac{(0.697-0.496)^2}{2\cdot 0.195^2}\right)\approx1.203$
$p(含糖率=0.460|好瓜=是)=\frac{1}{\sqrt{2\pi}0.101}\exp\left(-\frac{(0.460-0.279)^2}{2\cdot 0.101^2}\right)\approx0.788$	$p(含糖率=0.460|好瓜=否)=\frac{1}{\sqrt{2\pi}0.108}\exp\left(-\frac{(0.460-0.154)^2}{2\cdot 0.108^2}\right)\approx0.066$

条件（似然）概率	是（好瓜）	否(好瓜）
青绿(色泽)	$P_{青绿|是}=P({色泽=青绿|好瓜=是})=\frac{3}{8}=0.375$	$P_{青绿|否}=P({色泽=青绿|好瓜=否})=\frac{3}{9}=0.333$
乌黑(色泽)	$P_{乌黑|是}=P({色泽=乌黑|好瓜=是})=\frac{4}{8}=0.500$	$P_{乌黑|否}=P({色泽=乌黑|好瓜=否})=\frac{2}{9}=0.222$
浅白(色泽)	$P_{浅白|是}=P({色泽=浅白|好瓜=是})=\frac{1}{8}=0.125$	$P_{浅白|否}=P({色泽=浅白|好瓜=否})=\frac{4}{9}=0.444$

于是，根据朴素贝叶斯公式可得：

$$P(\text{好瓜=是})\times P_\text{青绿|是} \times P_\text{蜷缩|是}\times P_\text{浊响|是}\times P_\text{清晰|是}\times P_\text{凹陷|是}\times P_\text{硬滑|是}\times p_{\text{密度:0.697|是}}\times p_{\text{含糖:0.460|是}}\approx0.063$$

$$P(\text{好瓜=否})\times P_\text{青绿|否} \times P_\text{蜷缩|否}\times \text{浊响|否}\times P_\text{清晰|否}\times P_\text{凹陷|否}\times P_\text{硬滑|否}\times p_{\text{密度:0.697|否}}\times p_{\text{含糖:0.460|否}}\approx 6.80 \times10^{-5}$$

由于$0.063>6.80\times10^{-5}$，因此，预测结果为“好瓜”。

2.2 拉普拉斯平滑

若某个属性值在训练集中没有与某个类同时出现过，则直接基于朴素贝叶斯进行估计，很有很能出现错误的估计：例如：对一个“敲声=清脆”的测试例，有

$$P(\text{清脆|是})=P(\text{敲声=清脆|好瓜=是})=\frac{0}{8}=0$$

为了避免其他属性携带的信息被训练集中未出现的属性值“抹去”，在估计概率值时通常要进行“平滑”，常用“拉普拉斯修正”，令$N$表示训练集$D$中可能的类别数，$N_i$​表示第$i$个属性可能的取值数，则修正后的$P(c)$和$P(x_i∣c)$如下

$$\hat{P(c)}=\frac{|D_c|+1}{|D|+N} \quad \hat{P(x_i|c)}=\frac{|D_{c,x_i}|+1}{|D_c|+N_i}​$$

也可采用下面扩展公式平滑修正计算，此公式添加了调整参数$\lambda$，当$\lambda=1$时和上式一致。

$$\begin{aligned} & P_\lambda\left(y=c_n\right)=\frac{\sum_{i=1}^N I\left(y_i=c_n\right)+\lambda}{N+K \lambda},\quad n=1,2, \ldots, K \\ & P_\lambda\left(x^j=a_j \mid y=c_n\right)=\frac{\sum_{i=1}^N I\left(x_i^j=a_j \mid y=c_n\right)+\lambda}{\sum_{i=1}^N I\left(y_i=c_n\right)+S_j \lambda} \end{aligned}$$

其中$x_i^j$是第i个样本的第j个特征；$a_j$是第j个特征可能取的某个值；$I$指示函数或示性函数（indicator function）。
考虑下表15个数据的贝叶斯分类，采用拉普拉斯平滑修正后的先验概率和条件概率(似然概率)，计算过程如下：

Number	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15
$x^1$	1	1	1	1	1	2	2	2	2	2	3	3	3	3	3
$x^2$	S	M	M	S	S	S	M	M	L	L	L	M	M	L	L
y	-1	-1	1	1	-1	-1	1	1	1	1	1	1	1	1	-1

$$P(x=(2, S) \mid y=1)=P(y=1) * P\left(x^1=2 \mid y=1\right) * P\left(x^2=S \mid y=1\right)$$

因为 $P(y=1)=\frac{9+1}{15+\left[2]*[1]\right.}[2$ 表示 $y$ 有 2 个取值, $\lambda为1]$

$$P\left(x^1=2 \mid y=1\right)=\frac{3+1}{9+[3]^* 1} \quad\left[3\right. 表示 x^1 有 3 个取值 ]$$

$$P\left(x^2=S \mid y=1\right)=\frac{1+1}{9+[3]^* 1} \quad\left[3\right. 表示 x^2 有 3 个取值 ]$$

因此 $P(x=(2, S) \mid y=1)=\frac{10}{17} * \frac{4}{12} * \frac{2}{12}=\frac{5}{153}$

$$P(x=(2, S) \mid y=-1)=P(y=-1) * P\left(x^1=2 \mid y=-1\right) * P\left(x^2=S \mid y=-1\right)$$

因为 $P(y=-1)=\frac{6+[1]}{15+[2]^*[1]}[2$ 表示 $y$ 有 2 个取值, $\lambda为1]$

$$P\left(x^1=2 \mid y=-1\right)=\frac{2+1}{6+[3]^* 1} \quad\left[3\right. 表示 x^1 有 3 个取值 ]$$

$$P\left(x^2=S \mid y=-1\right)=\frac{3+1}{6+[3]^* 1} \quad\left[3\right. 表示 x^2 有 3 个取值]$$

因此 $P(x=(2, S) \mid y=-1)=\frac{7}{17} * \frac{3}{9} * \frac{4}{9}=\frac{28}{459}$

三、朴素贝叶斯的R实现

3.1 似然R函数

#构造训练集
data <- matrix(c("sunny","hot","high","weak","no",
                 "sunny","hot","high","strong","no",
                 "overcast","hot","high","weak","yes",
                 "rain","mild","high","weak","yes",
                 "rain","cool","normal","weak","yes",
                 "rain","cool","normal","strong","no",
                 "overcast","cool","normal","strong","yes",
                 "sunny","mild","high","weak","no",
                 "sunny","cool","normal","weak","yes",
                 "rain","mild","normal","weak","yes",
                 "sunny","mild","normal","strong","yes",
                 "overcast","mild","high","strong","yes",
                 "overcast","hot","normal","weak","yes",
                 "rain","mild","high","strong","no"), byrow = TRUE,
               dimnames = list(day = c(),
                               condition = c("outlook","temperature",
                                             "humidity","wind","playtennis")), nrow=14, ncol=5);

#计算先验概率
prior.yes = sum(data[,5] == "yes") / length(data[,5]);
prior.no  = sum(data[,5] == "no")  / length(data[,5]);

#模型
naive.bayes.prediction <- function(condition.vec) {
  # Calculate unnormlized posterior probability for playtennis = yes.
  playtennis.yes <-
    sum((data[,1] == condition.vec[1]) & (data[,5] == "yes")) / sum(data[,5] == "yes") * # P(outlook = f_1 | playtennis = yes)
    sum((data[,2] == condition.vec[2]) & (data[,5] == "yes")) / sum(data[,5] == "yes") * # P(temperature = f_2 | playtennis = yes)
    sum((data[,3] == condition.vec[3]) & (data[,5] == "yes")) / sum(data[,5] == "yes") * # P(humidity = f_3 | playtennis = yes)
    sum((data[,4] == condition.vec[4]) & (data[,5] == "yes")) / sum(data[,5] == "yes") * # P(wind = f_4 | playtennis = yes)
    prior.yes; # P(playtennis = yes)
  
  # Calculate unnormlized posterior probability for playtennis = no.
  playtennis.no <-
    sum((data[,1] == condition.vec[1]) & (data[,5] == "no"))  / sum(data[,5] == "no")  * # P(outlook = f_1 | playtennis = no)
    sum((data[,2] == condition.vec[2]) & (data[,5] == "no"))  / sum(data[,5] == "no")  * # P(temperature = f_2 | playtennis = no)
    sum((data[,3] == condition.vec[3]) & (data[,5] == "no"))  / sum(data[,5] == "no")  * # P(humidity = f_3 | playtennis = no)
    sum((data[,4] == condition.vec[4]) & (data[,5] == "no"))  / sum(data[,5] == "no")  * # P(wind = f_4 | playtennis = no)
    prior.no; # P(playtennis = no)
  
  return(list(post.pr.yes = playtennis.yes,
              post.pr.no  = playtennis.no,
              prediction  = ifelse(playtennis.yes >= playtennis.no, "yes", "no")));
}

#预测
naive.bayes.prediction(c("rain",     "hot",  "high",   "strong"));
naive.bayes.prediction(c("sunny",    "mild", "normal", "weak"));
naive.bayes.prediction(c("overcast", "mild", "normal", "weak"));

3.2 调用e1071包naiveBayes函数

naiveBayes(formula, data, laplace = 0, ..., subset, na.action = na.pass)

参数	解释
formula	类似一般线性回归表达式，不含常数项
data	需要分析的训练数据对象
laplace	拉普拉斯估计值，默认为0
subset	抽取要分析的训练数据子集
na.action	缺失值的处理方法。默认情况下不将缺失值纳入模型计算，如果设定为na.omit则会删除缺失值进行计算

library(e1071)  #包中有naiveBayes函数
##离散数据
data(HouseVotes84, package = "mlbench")
# 构建模型
bayes_model <- naiveBayes(Class ~ ., data = HouseVotes84)
# 进行预测判别
predict(bayes_model, HouseVotes84[1:5,])
# 获取各样本的概率值
predict(bayes_model, HouseVotes84[1:5,], type = "raw")
# 获取模型的预测精度
pred <- predict(bayes_model, HouseVotes84)
table(pred, HouseVotes84$Class)
# 使用laplace校准
bayes_model_laplace <- naiveBayes(Class ~ ., data = HouseVotes84, laplace = 3)

##连续数据
bayes_iris <- naiveBayes(Species ~ ., data = iris)
iris_pred <- predict(bayes_iris,iris)
table(iris_pred,iris$Species)  

3.3 调用klaR包NaiveBayes函数

NaiveBayes（）
默认情况：Naivebayes(x,grouping,prior,usekernd=FALSE,fL=0,…)
数据对象为公式：Naivebayes(formula,data,…,subset,na.nation=na.pass)

参数	解释
x	要处理的数据库data.frame或者数据矩阵matrix
formula	放置生成判别规则的公式
data/subset	以formula为对象的函数格式中，分别用于指明该formula中变量所来自的数据集名称和纳入规则建立规程的样本
grouping	指明每个观测样本属于的类别
prior	设置各类别的先验概率
na.action	默认na.pass不会将缺失值纳入计算，不会以你选哪个函数运行，取值为ma.omit时便是删除相应的含有缺失值的观测变量样本。

library(klaR)
#数据集分割
index <- sample(2,size = nrow(iris),replace = TRUE,prob = c(0.7,0.3))
train_data <- iris[index == 1,]
test_data <- iris[index == 2,]
# 构建贝叶斯模型
Bayes_model <- NaiveBayes(Species ~ ., data = train_data)
# 进行预测
Bayes_model_pre <- predict(Bayes_model, newdata = test_data[,1:4])
# 生成实际与预判交叉表
table(test_data$Species,Bayes_model_pre$class)
#生成预判精度
sum(diag(table(test_data$Species,Bayes_model_pre$class)))/sum(table(test_data$Species,Bayes_model_pre$class))

总结

最为广泛的两种分类模型是决策树模型(Decision Tree Model)和朴素贝叶斯模型（Naive Bayesian Model，NBM）。和决策树模型相比，朴素贝叶斯分类器(Naive Bayes Classifier 或 NBC)有着坚实的数学基础，以及稳定的分类效率；NBC模型所需估计的参数很少，对缺失数据不太敏感，算法也比较简单。理论上NBC模型与其他分类方法相比具有最小的误差率。但是NBC模型假设属性之间相互独立，在一定程度上降低了贝叶斯分类算法的分类效果，但在实际的应用场景中极大地简化了贝叶斯方法的复杂性。

参考文献

1. 朴素贝叶斯算法（带例题解释）
2. 模式识别课程(一):贝叶斯决策及python实现 莺尾花分类
3. 第三节 贝叶斯决策
4. 贝叶斯判别简要原理及其实例