朴素贝叶斯
1 朴素叶斯简介及基本思想
贝叶斯定理是由英国数学家贝叶斯(Thomas Bayes 1702-1761)创建的，被用来描述两个条件概率之间的关系，也是当前机器学习领域中的一个常见的有监督分类算法。在本小节中我们着重介绍朴素贝叶斯分类器(Naive Bayes Classifiers)的内容。
贝叶斯的基本思想就在于通过考虑特征概率来预测分类,朴素贝叶斯分类器是一个基于贝叶斯定理的比较简单的概率分类器，其中朴素(naive)二字是指的对于模型中各个特征有强独立性的假设，并未将特征间的相关性纳入考虑中。
在实际进行分类的过程中首先利用训练数据并根据特征的取值来计算每个类别被观察的概率。当分类器用于无标签数据时，分类器就会根据观测到的概率来预测新的特征最有可能属于哪个类。
需要特别注意的是，朴素贝叶斯算法是在以下述两个假设为前提的情况下建立的：
- 假设所有特征都具有相同的重要性和独立性(条件特征独立)
- 假设连续变量符合正态分布
这为后续的计算提供了很大的便利，但是这样的假设在大部分情况下不符合现实场景，因此会在一定程度上影响算法的精度，同时也带来了一些不适用性。
2 朴素贝叶斯原理及解释
朴素贝叶斯算法是基于概率的分类算法，其中涉及两个变量，分别是：
- 自变量(\(X\))：表示发生事件或样本的特征
- 因变量(\(Y\))：表示事件或样本所对应的类别
朴素贝叶斯的思想在于根据某个事件的先验概率来计算因变量属于某个类别的后验概率。朴素贝叶斯定理的公式展示如下： \[P(Y|X)=\frac{P(YX)}{P(X)}=\frac{P(X|Y)P(Y)}{P(X)}\] 从公式1中可以看出，如果要想得到事件\(X\)发生的情况下\(Y\)所属某个类别的概率，只需要计算以下几点：
- \(P(X)\)：\(X\)事件的概率，即\(X\)的先验概率；
- \(P(Y)\)：\(Y\)属于某类的概率，即\(Y\)的先验概率；
- \(P(X|Y)\)：以及已知\(Y\)的某个分类下，事件\(X\)的概率，即后验概率。
观察上述需求可知，不同水平的\(P(C)\)比较容易获得，只需计算数据中各个类别所对应的概率即可；所以贝叶斯中的难点就在于后验概率的计算，因此在这里我们便可以运用到上述的假设一，即假设所有特征都具有相同的重要性和独立性(条件特征独立)，在这种前提下，通过概率论的知识我们便可以用如下的公式很容易地计算出\(P(X|C)\)： \[P(X|C_i)=P(x_1|C_i)*P(x_2|C_i)*…*P(x_n|C_i)\]
在自变量为离散性变量的情况下，我们可以直接通过计算频率来得出先验概率，并通过公式1.2计算出后验概率。
但是对于连续性变量来说情况要稍微复杂一些，这时我们需要运用到假设二的内容，即假设连续变量符合正态分布。正态分布的概率密度函数如下所示： \[f(x)=\frac{1}{\sqrt{2\pi}\sigma}e^-\frac{(x-\mu)^2}{2\sigma^2}\]
其中\(\mu\)代表数学期望(均值)，\(\sigma\)代表标准差。
从该函数中可以看出只要变量的均值和标准差就能得出其对应的概率值。总之，朴素贝叶斯分类的过程就是通过上述过程计算出某事物归属于不同类的概率，再从中挑选出最大的概率，并将其归为相关类别。
3 案例
本部分主要介绍了朴素贝叶斯分类在R语言中的实现过程。本次实验采用了经典的Mushrooms数据集，目的通过数据集中给出的蘑菇的物理特性，结合朴素贝叶斯分类算法，来判断蘑菇属于有毒的还是可食用的类别。
案例主要涉及以下六个小部分：
- 数据集介绍
- 导入数据
- 数据探索与预处理
- 特征工程
- 朴素贝叶斯建模
- 模型评估
3.1 数据集介绍
数据集名称为Mushroom Data Set，数据集中包含从奥杜邦的野外指南中得到的描述蘑菇的物理特性。本数据集记录了有关伞菌属和环柄菇属蘑菇的对应假想样本的23种描述，每一个物种都被认为是可以食用的，绝对有毒的，或者是未知的可食性，不推荐食用。而且该指南明确指出，没有简单的规则来确定蘑菇的可食性。
数据有23列，8124行，由于变量数量太多，选择部分特征变量进行解释，展示如下：
| 变量名 | 变量解释 |
|---|---|
| class | 蘑菇类别：edible=e, poisonous=p |
| cap.shape | 菌盖形状：bell=b, conical=c… |
| cap.surface | 菌盖表面： fibrous=f, grooves=g… |
| cap-color | 菌盖颜色： brown=n, buff=b… |
| odor | 气味：almond=a, anise=l… |
| …… | …… |
| stalk.shape | 茎的形状：enlarging=e, tapering=t |
| stalk.surface.above.ring | 菌环上方茎表面：fibrous=f scaly=y… |
| ring.number | 菌环数量：none=n, one=o, two=t |
3.2 导入数据
本次实验使用kaggle的mushrooms数据集，通过read.csv()函数获取数据集：
setwd("E:/") # 设置工作路径为E盘
# 使用read.csv()函数导入mushroom数据集，并命名为mydata
mydata <- read.csv('mushrooms.csv')3.3 数据探索与预处理
在本步骤中进行数据的展示和简单的数据探索，并根据具体情况进行数据预处理。
注意：如果没有下载plyr和dplyr程序包，请先使用install.packages("plyr")命令和install.packages("dplyr")下载程序包之后再加载。
# install.packages("plyr") # 下载plyr程序包
# install.packages("dplyr") # 下载dplyr程序包
library(plyr) # 加载plyr程序包
library(dplyr) # 加载dplyr程序包，使用glimpse()函数
glimpse(mydata) # 查看数据的基本结构## Observations: 8,124
## Variables: 23
## $ class <fctr> p, e, e, p, e, e, e, e, p, e, e, e, ...
## $ cap.shape <fctr> x, x, b, x, x, x, b, b, x, b, x, x, ...
## $ cap.surface <fctr> s, s, s, y, s, y, s, y, y, s, y, y, ...
## $ cap.color <fctr> n, y, w, w, g, y, w, w, w, y, y, y, ...
## $ bruises <fctr> t, t, t, t, f, t, t, t, t, t, t, t, ...
## $ odor <fctr> p, a, l, p, n, a, a, l, p, a, l, a, ...
## $ gill.attachment <fctr> f, f, f, f, f, f, f, f, f, f, f, f, ...
## $ gill.spacing <fctr> c, c, c, c, w, c, c, c, c, c, c, c, ...
## $ gill.size <fctr> n, b, b, n, b, b, b, b, n, b, b, b, ...
## $ gill.color <fctr> k, k, n, n, k, n, g, n, p, g, g, n, ...
## $ stalk.shape <fctr> e, e, e, e, t, e, e, e, e, e, e, e, ...
## $ stalk.root <fctr> e, c, c, e, e, c, c, c, e, c, c, c, ...
## $ stalk.surface.above.ring <fctr> s, s, s, s, s, s, s, s, s, s, s, s, ...
## $ stalk.surface.below.ring <fctr> s, s, s, s, s, s, s, s, s, s, s, s, ...
## $ stalk.color.above.ring <fctr> w, w, w, w, w, w, w, w, w, w, w, w, ...
## $ stalk.color.below.ring <fctr> w, w, w, w, w, w, w, w, w, w, w, w, ...
## $ veil.type <fctr> p, p, p, p, p, p, p, p, p, p, p, p, ...
## $ veil.color <fctr> w, w, w, w, w, w, w, w, w, w, w, w, ...
## $ ring.number <fctr> o, o, o, o, o, o, o, o, o, o, o, o, ...
## $ ring.type <fctr> p, p, p, p, e, p, p, p, p, p, p, p, ...
## $ spore.print.color <fctr> k, n, n, k, n, k, k, n, k, k, n, k, ...
## $ population <fctr> s, n, n, s, a, n, n, s, v, s, n, s, ...
## $ habitat <fctr> u, g, m, u, g, g, m, m, g, m, g, m, ...summary(mydata) # 查看变量的主要描述性统计量## class cap.shape cap.surface cap.color bruises odor
## e:4208 b: 452 f:2320 n :2284 f:4748 n :3528
## p:3916 c: 4 g: 4 g :1840 t:3376 f :2160
## f:3152 s:2556 e :1500 s : 576
## k: 828 y:3244 y :1072 y : 576
## s: 32 w :1040 a : 400
## x:3656 b : 168 l : 400
## (Other): 220 (Other): 484
## gill.attachment gill.spacing gill.size gill.color stalk.shape
## a: 210 c:6812 b:5612 b :1728 e:3516
## f:7914 w:1312 n:2512 p :1492 t:4608
## w :1202
## n :1048
## g : 752
## h : 732
## (Other):1170
## stalk.root stalk.surface.above.ring stalk.surface.below.ring
## ?:2480 f: 552 f: 600
## b:3776 k:2372 k:2304
## c: 556 s:5176 s:4936
## e:1120 y: 24 y: 284
## r: 192
##
##
## stalk.color.above.ring stalk.color.below.ring veil.type veil.color
## w :4464 w :4384 p:8124 n: 96
## p :1872 p :1872 o: 96
## g : 576 g : 576 w:7924
## n : 448 n : 512 y: 8
## b : 432 b : 432
## o : 192 o : 192
## (Other): 140 (Other): 156
## ring.number ring.type spore.print.color population habitat
## n: 36 e:2776 w :2388 a: 384 d:3148
## o:7488 f: 48 n :1968 c: 340 g:2148
## t: 600 l:1296 k :1872 n: 400 l: 832
## n: 36 h :1632 s:1248 m: 292
## p:3968 r : 72 v:4040 p:1144
## b : 48 y:1712 u: 368
## (Other): 144 w: 192glimpse()函数展示了数据的观察长度(Observations)，值为8124；变量数(Variables)为23；及各个变量的名称、属性和其对应的值。
通过观察summary()的输出结果可知，数据中的第17列(veil.type)只存在一个值。恒定变量在数据中是没有任何作用的，所以接下来要剔除此变量。
mydata$veil.type <- NULL # 剔除此变量3.4 分割训练集与测试集
使用creat包中的createDataPartition()函数对数据进行分层抽样，以3:7的比例划分测试集和训练集，并分别保存在train和test两个变量中，设置随机种子以保证实验结果相同。
# install.packages("caret")
library(caret)
set.seed(12) # 设置随机种子，保持随机抽样一致性
# 对数据集进行分层抽样，设置训练集占比为70%
index <- createDataPartition(1:nrow(mydata), p = 0.7, list = FALSE)
train <- mydata[index,] # 抽取70%的数据作为训练集
test <- mydata[-index,] # 抽取30%的数据作为测试集使用prop.table()函数查看抽样与总体之间是否吻合。prop.table()函数在此处计算出了每个类别所对应的频率。
prop.table(table(mydata$class)) # 统计每个类别所对应的频率##
## e p
## 0.5179714 0.4820286prop.table(table(train$class)) # 统计每个类别所对应的频率##
## e p
## 0.5179325 0.4820675prop.table(table(test$class)) # 统计每个类别所对应的频率##
## e p
## 0.5180624 0.4819376可以看出mydata，train，test返回的值大致相同，因此抽样与总体之间相较吻合。
3.5 特征工程
原数据为23列8124行的数据，进行预处理删除一列，除去因变量class后依旧剩下21列特征，特征数量依旧过多，我们需要选择有意义的特征输入机器学习的算法和模型进行训练。在这里我们使用caret包中方法进行特征选择。
rfeControls.nb <- rfeControl( # 设置rfe函数的控制参数
functions = nbFuncs, # 使用朴素贝叶斯算法
method = 'cv', # 使用交叉验证的抽样方法
number = 5) # 指定5折交叉验证
fs.nb <- rfe(x = train[,-1], # 指定输入变量
y = train[,1], # 指定输出变量
sizes = seq(4,21,2), # 指定保留特征的数量
rfeControl = rfeControls.nb) # 指定控制参数在这里我们使用封装发的特征递减发来选择变量，函数及参数说明如下：
- 使用rfeControl()函数来设置rfe()函数的控制参数，其中采用贝叶斯函数(nbFuncs)和5次10则交叉验证的方法进行设置。
- 通过rfe()函数进行特征选择，其中x指定输入变量，y指定输出变量，sizes指定需要保留的特征项，并待入上述设置的控制参数。
plot(fs.nb, class = c('g','o')) # 通过可视化查看选择效果
fs.nb$optVariables # 展示筛选变量## [1] "gill.color" "bruises"
## [3] "gill.size" "ring.type"
## [5] "stalk.root" "stalk.surface.above.ring"
## [7] "stalk.surface.below.ring" "gill.spacing"
## [9] "habitat" "population"
## [11] "stalk.color.above.ring" "stalk.color.below.ring"
## [13] "cap.surface" "odor"
## [15] "ring.number" "stalk.shape"特征选择后通过plot()函数可视化可以看出当选择16个变量时模型效果达到最优，并通过fs.nb$optVariables展示出选择的16个特征变量。
3.6 朴素贝叶斯建模
本部分主要介绍和朴素贝叶斯核心函数的使用方法以及对蘑菇(mushrooms)数据的建模预测，细化介绍如下：
3.6.1 NaiveByaes()函数简介
在进行建模先对贝叶斯的核心函数进行介绍，在本次实验中使用的是klaR包中的NaiveByaes()函数，其主要作用是使用贝叶斯规则计算给定独立变量的分类变量的条件A后验概率。该函数总共有两种构建格式，分别是公式formula格式和默认格式。
NaiveBayes(formula, data, …, subset, na.action = na.pass)
NaiveBayes(x, grouping, prior, usekernel = FALSE, fL = 0, …)
下面列举出NaiveByaes()函数的常用参数及其对应的解释：
- formula：指定参与模型计算的变量，以公式形式给出，类似于y=x1+x2+x3；
- data：用于指定需要分析的数据对象；
- na.action：指定缺失值的处理方法，默认情况下不将缺失值纳入模型计算，也不会发生报错信息，当设为“na.omit”时则会删除含有缺失值的样本；
- x：指定需要处理的数据，可以是数据框形式，也可以是矩阵形式；
- grouping：为每个观测样本指定所属类别；
- prior：可为各个类别指定先验概率，默认情况下用各个类别的样本比例作为先验概率；
- usekernel：指定密度估计的方法(在无法判断数据的分布时，采用密度密度估计方法)，默认情况下使用正态分布密度估计，设为TRUE时，则使用核密度估计方法；
- fL：指定是否进行拉普拉斯修正，默认情况下不对数据进行修正，当数据量较小时，可以设置该参数为1，即进行拉普拉斯修正。
3.6.2 数据建模
上面我们已经在特征工程中筛选出了重要的变量，下面我们就利用上面筛选的数据进行建模：
- 将经过特征选择的变量和因变量class放入vars中，用于选择变量；
- 使用NaiveBayes()的函数形式进行建模，并将模型赋值给fit变量
# install.packages("klaR")
library(klaR)
vars <- c('class',fs.nb$optVariables) # 选择特征变量
fit <- NaiveBayes(class ~ ., data = train[,vars]) # 拟合朴素贝叶斯模型3.6.3 模型输出量解释
在上一步中，模型被记为fit，在这一步骤中我们将对模型常用的输出项及其对应的意义做说明：
3.6.3.1 names()
显示出fit中所包含的输出项的名称
names(fit) # 展示输出项的名称## [1] "apriori" "tables" "levels" "call" "x" "usekernel"
## [7] "varnames"3.6.3.2 model$apriori
输出中的apriori记录了该次执行过程中使用的先验概率
fit$apriori # 展示先验概率## grouping
## e p
## 0.5179325 0.4820675从apriori项中可以看出，标记为e(即可食用)的蘑菇占51.78%;标记为p(即有毒)的蘑菇占48.21%。
3.6.3.3 model$tables
输出中的tables记录了建立判别规则的所有变量在各类别下的条件概率(后验概率)
fit$tables # 展示后验概率## $gill.color
## var
## grouping b e g h k n
## e 0.0000000 0.02172437 0.05940258 0.04786151 0.08146640 0.21554650
## p 0.4427425 0.00000000 0.13092633 0.13785558 0.01458789 0.02917578
## var
## grouping o p r u w y
## e 0.01629328 0.1999321 0.00000000 0.10488798 0.23693143 0.015953836
## p 0.00000000 0.1630197 0.00547046 0.01130562 0.06017505 0.004741065
##
## $bruises
## var
## grouping f t
## e 0.3452138 0.6547862
## p 0.8409920 0.1590080
##
## $gill.size
## var
## grouping b n
## e 0.9270197 0.07298031
## p 0.4361780 0.56382203
##
## $ring.type
## var
## grouping e f l n p
## e 0.2348948 0.01221996 0.0000000 0.000000000 0.7528853
## p 0.4511306 0.00000000 0.3351568 0.008388038 0.2053246
##
## $stalk.root
## var
## grouping ? b c e r
## e 0.1704005 0.4582485 0.121860149 0.20298710 0.04650373
## p 0.4496718 0.4762947 0.009846827 0.06418673 0.00000000
##
## $stalk.surface.above.ring
## var
## grouping f k s y
## e 0.09877800 0.03360489 0.8642227 0.003394433
## p 0.03501094 0.57111597 0.3924143 0.001458789
##
## $stalk.surface.below.ring
## var
## grouping f k s y
## e 0.10828242 0.03462322 0.8071962 0.04989817
## p 0.03719912 0.55981036 0.3862144 0.01677608
##
## $gill.spacing
## var
## grouping c w
## e 0.7145282 0.28547183
## p 0.9755653 0.02443472
##
## $habitat
## var
## grouping d g l m p u
## e 0.4490835 0.3285811 0.05804481 0.063815343 0.03360489 0.02342159
## p 0.3129103 0.1936543 0.15426696 0.009117433 0.26258206 0.06746900
## var
## grouping w
## e 0.04344874
## p 0.00000000
##
## $population
## var
## grouping a c n s v y
## e 0.08791582 0.06653089 0.09368635 0.21215207 0.2875085 0.2522064
## p 0.00000000 0.01203501 0.00000000 0.09153902 0.7279358 0.1684902
##
## $stalk.color.above.ring
## var
## grouping b c e g n o
## e 0.0000000 0.000000000 0.02104549 0.140869 0.003394433 0.04616429
## p 0.1126915 0.008388038 0.00000000 0.000000 0.109409190 0.00000000
## var
## grouping p w y
## e 0.1350984 0.6534284 0.000000000
## p 0.3347921 0.4332604 0.001458789
##
## $stalk.color.below.ring
## var
## grouping b c e g n o
## e 0.0000000 0.000000000 0.02104549 0.1337407 0.01527495 0.04616429
## p 0.1119621 0.008388038 0.00000000 0.0000000 0.12071481 0.00000000
## var
## grouping p w y
## e 0.1347590 0.6490156 0.000000000
## p 0.3249453 0.4296134 0.004376368
##
## $cap.surface
## var
## grouping f g s y
## e 0.3740665 0.0000000000 0.2671419 0.3587916
## p 0.1907367 0.0007293946 0.3650620 0.4434719
##
## $odor
## var
## grouping a c f l m n
## e 0.09911745 0.00000000 0.0000000 0.09470468 0.000000000 0.80617787
## p 0.00000000 0.04631656 0.5576222 0.00000000 0.008388038 0.03026988
## var
## grouping p s y
## e 0.00000000 0.0000000 0.0000000
## p 0.06418673 0.1487965 0.1444201
##
## $ring.number
## var
## grouping n o t
## e 0.000000000 0.8761032 0.12389681
## p 0.008388038 0.9719183 0.01969365
##
## $stalk.shape
## var
## grouping e t
## e 0.385947 0.614053
## p 0.484318 0.515682上述输出中展示了在不同特征中的蘑菇分别为有毒或可食用的概率大小，比如当菌褶颜色(gill.color)为浅黄色(b)时，它可食用的概率为0，有毒的概率为0.44。
3.6.3.4 others
本项中包含判别变量等级项,判别命令项，是否使用标准密度估计，特征变量名等输出项
fit$levels # 判别变量等级项## [1] "e" "p"输出显示存在两个levels，分别为e和p。
fit$call # 判别命令项## NaiveBayes.default(x = X, grouping = Y)输出模型的判别命令项
fit$usekernel # 是否使用标准密度估计## [1] FALSE输出表示模型没有使用标准密度估计
fit$varnames # 特征变量名## [1] "gill.color" "bruises"
## [3] "gill.size" "ring.type"
## [5] "stalk.root" "stalk.surface.above.ring"
## [7] "stalk.surface.below.ring" "gill.spacing"
## [9] "habitat" "population"
## [11] "stalk.color.above.ring" "stalk.color.below.ring"
## [13] "cap.surface" "odor"
## [15] "ring.number" "stalk.shape"输出展示出了用于构建模型的16个特征变量名称。
3.6.4 模型预测
对测试集使用predict()函数，进行模型预测:
- 在newdata中带入先前划分好的测试集，即test数据，同时依旧在test中选择用vars保留的特征列，并删除位于第一列的因变量class列,将模型的预测结果赋值给pred变量。
- 输出结果包括各样本的预测类别class和预测过程中各样本属于每一个类别的后验概率posterior。
pred <- predict(fit, newdata = test[,vars][,-1]) # 对测试集进行预测pred的输出结果包括样本的预测类别和预测样本中属于每一类别的后验概率(posterior),其中后验概率大得一项被判定为相应类别。
3.7 模型评估
有很多标准可以用来衡量模型的好坏，本次实验采用了两种方式对朴素贝叶斯分类器进行评估，分别是混淆矩阵(Confusion Matrix)和ROC曲线(Receiver Operating Characteristic)，具体步骤展示如下：
3.7.1 Confusion Matrix
混淆矩阵也称误差矩阵，是表示模型精度的一种格式，其结果通过将每个实测像元的位置和分类与分类图像中的相应位置和分类相比较计算得出。
通过table()函数建立混淆矩阵来展示预测现象，其中第一个元素放置模型预测出的类别(pred$class)；第二个元素放置测试集真实的标签列(test[,1])。
freq <- table(pred$class, test[,1]) # 构建混淆矩阵
freq # 展示混淆矩阵##
## e p
## e 1262 134
## p 0 1040通过混淆矩阵可以看出，可食用得蘑菇被全部预测正确，但是有92种有毒得蘑菇被错误预测为可食用的蘑菇。下一步通过混淆矩阵的内容计算正确率。
accuracy <- sum(diag(freq))/sum(freq) # 计算正确率
accuracy # 展示正确率## [1] 0.9449918可见模型正确率为95.4%，具有较高的分类准确率。
3.7.2 ROC curve
ROC曲线和AUC值经常被用来作为二分类模型的评价标准，本部分通过pROC包中的roc()计算出模型的AUC值，并使用plot()函数绘制ROC曲线，下面简单介绍一下ROC曲线和AUC值：
- ROC曲线全称为受试者工作特征曲线(receiver operating characteristic curve)，是根据一系列二分类器，以真阳性率为纵坐标，假阳性率为横坐标绘制的曲线。它能够容易地查出在不同界限值时的模型的识别能力和模型的准确率。
- AUC(Area Under Curve)值为ROC曲线下的面积，AUC值越大，当前的分类算法越有可能将正样本排在负样本前面，即能够更好的分类。
# install.packages("pROC")
library(pROC)
# 模型的AUC值
modelroc <- roc(as.integer(test[,1]),
as.integer(factor(pred$class)))
# 绘制ROC曲线
plot(modelroc, print.auc = TRUE, auc.polygon = TRUE,
grid = c(0.1,0.2), grid.col = c('green','red'),
max.auc.polygon = TRUE, auc.polygon.col = 'steelblue')
由以上展示的ROC曲线可以看出不同界限下模型识别能力的变化；以及模型的AUC值为0.953，分类效果较佳。
4 总结
朴素贝叶斯是一个高效且易于实现的算法，对小规模的数据表现很好，适合多分类任务，对缺失值也不是那么敏感；但是朴素贝叶斯算法也存在本身的限制，它对输入数据的表达形式很敏感，而且由于其对变量之间相互独立的假设在现实场景中很难实现，所以也会对分类的准确率带来一定影响。因此在使用朴素贝叶斯分类器时也要考虑是否是其适应的场景。
5 参考文献
[1]黄文, 王正林. 数据挖掘 : R语言实战[M]. 电子工业出版社, 2014.