朴素贝叶斯算法

  朴素贝叶斯(naive bayes)算法是基于贝叶斯定理与特征条件独立假设的分类方法。对于给定的训练数据集，首先基于特征条件独立假设学习输入、输出的联合概率分布，然后基于此模型，对给定的输入$x$，利用贝叶斯定理求出后验概率最大的输出$y$。朴素贝叶斯法实现简单，学习与预测的效率都很高，是一种常用方法。

数学原理

  设输入空间$\mathcal X \subseteq R^n$为$n$维向量的集合，输出空间为类标记集合$\mathcal Y={c_1,c_2,…,c_K}$。输入为特征向量$x \in \mathcal X$，输出为类标记$y \in \mathcal Y$。$X$是定义在输入空间$\mathcal X$上的随机变量，$Y$是定义在输出空间$\mathcal Y$上的随机变量。$P(X,Y)$是$X$和$Y$的联合概率分布。训练数据集

$$T={(x_1,y_1),(x_2,y_2),\ldots,(x_N,y_N)}$$
由$P(X,Y)$独立同分布产生。

  朴素贝叶斯法通过训练数据集学习联合概率分布$P(X,Y)$。具体地，学习以下先验概率分布及条件概率分布。先验概率分布，

$$P(Y=c_K),\quad k=1,2,\ldots,K$$

条件概率分布

$$P(X=x|Y=c_k)=P(X^{(1)}=x^{(1)},…,X^{(n)}=x^{(n)}),\quad k=1,2,\ldots,K$$

于是学习到联合概率分布

$$P(X,Y)$$

  条件概率分布$P(X=x|Y=c_k)$有指数量级的参数，其估计实际是不可行的。事实上，假设$x^{(j)}$可取值有$S_j$个，$j=1,2,\ldots,n$，$Y$可取值有$K$个，那么参数个数为$K \prod_{j=1}^{n}S_j$。

  朴素贝叶斯算法对条件概率分布作了条件独立的假设，即用于分类的特征在类别确定的条件下都是条件独立的。

$$\begin{split} P(X=x|Y=c_k) &= P(X^{(1)}=x^{(1)},\ldots,X^{(n)}=x^{(n)}|Y=c_k)\\ &=\prod_{j=1}^{n}P(X^{(j)}=x^{(j)}|Y=c_k) \end{split} \tag{1}$$
这一假设使朴素贝叶斯算法变得简单，但是有时会牺牲一定的分类准确率。

  朴素贝叶斯法分类时，对给定的输入$x$，通过学习到的模型计算后验概率分布$P(Y=c_k|X=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)} \tag{2}$$
把1式代入2式，有

$$P(Y=c_k|X=x)=\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)}$$

所以输出为

$$y=f(x)=\mathop{\arg\max}_{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)} \tag{3}$$

由于3式中分母对所有$c_k$都相同，所以

$$y=f(x)=\mathop{\arg\max}_{c_k} P(Y=c_k)\prod_{j=1}^{n}P(X^{(j)}=x^{(j)}|Y=c_k) \tag{4}$$

参数估计

极大似然估计

  在朴素贝叶斯算法中，学习意味着估计$P(Y=c_k)$和$P(X^{(j)}=x^{(j)}|Y=c_k)$。可以使用极大似然估计法估计相应的概率。先验概率$P(Y=c_k)$的极大似然估计是

$$P(Y=c_k)=\frac{\sum_{i=1}^{N}I(y_i=c_k)}{N},\quad k=1,2,\ldots,K$$

设第$j$个特征$x^{(j)}$可能取值的集合为$\lbrace a_{j1},a_{},\ldots,a_{jS_j} \rbrace$，条件概率$P(X^{(j)}=a_{jl}|Y=c_k)$的极大似然估计是

$$P(X^{(j)}=a_{jl}|Y=c_k)= \frac{\sum_{i=1}^{N}I(x^{(j)}_{i}=a_{jl},y_i=c_k)}{\sum_{i=1}^{N}I(y_i=c_k)}$$
式中，$x_{i}^{(j)}$是第$i$个样本的第$j$个特征；$a_{jl}$是第$j$个特征可能取的第$l$个值；$I$为指示函数。

贝叶斯估计

  用极大似然估计可能出现所要估计的值为$0$的情况，这时会影响到后验概率的计算结果，使分类产生偏差。解决这一问题的方法是采用贝叶斯估计，条件概率的贝叶斯估计是

$$P_{\lambda}(X^{(j)}=a_{jl}|Y=c_k)=\frac{\sum_{i=1}^{N}I(x_{i}^{(j)}=a_{jl},y_i=c_k)+\lambda}{\sum_{i=1}^{N}I(y_i=c_k)+S_j\lambda}$$
式中$\lambda \geq 0$。等价于在随机变量各个取值的频数上赋予一个正数$\lambda>0$。当$\lambda=0$时就是极大似然估计。常取$\lambda=1$，这时称为拉普拉斯平滑(Laplace smoothing)。显然对于任何$l=1,2,\ldots,S_j,\quad k=1,2,\ldots,K$，有

$$P_{\lambda}(X^{(j)}=a_{jl}|Y=c_k)>0\\ \sum_{i=1}^{S_{jl}}P(X^{(j)}=a_{jl}|Y=c_k)=1$$
同样，先验概率的贝叶斯估计是
$$P_{\lambda}(Y=c_k)=\frac{\sum_{i=1}^{N}I(y_i=c_k)+\lambda}{N+K\lambda}$$

示例

  训练数据如下表，学习一个朴素贝叶斯分类器并确定$x=(2,S)^T$的类标记$y$，表中$X^{(1)}、X^{(2)}$为特征，取值的集合分别为$A_1={1,2,3},A_2={S,M,L}$，$Y$为类标记，$Y \in C=\lbrace 1,-1 \rbrace$

	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(Y=1)=\frac{9+1}{15+2}=\frac{10}{17},P(Y=-1)=\frac{6+1}{15+2}=\frac{7}{17}$$
$$P(X^{(1)}=1|Y=1)=\frac{2+1}{9+3}=\frac{3}{12},P(X^{(1)}=2|Y=1)=\frac{3+1}{9+3}=\frac{4}{12},P(X^{(1)}=3|Y=1)=\frac{4+1}{9+3}=\frac{5}{12}$$
$$P(X^{(1)}=1|Y=-1)=\frac{3+1}{6+3}=\frac{4}{9},P(X^{(1)}=2|Y=-1)=\frac{2+1}{6+3}=\frac{3}{9},P(X^{(1)}=3|Y=-1)=\frac{1+1}{6+3}=\frac{2}{9}$$
$$P(X^{(2)}=S|Y=1)=\frac{1+1}{9+3}=\frac{2}{12},P(X^{(2)}=M|Y=1)=\frac{4+1}{9+3}=\frac{5}{12},P(X^{(2)}=L|Y=1)=\frac{4+1}{9+3}=\frac{5}{12}$$
$$P(X^{(2)}=S|Y=-1)=\frac{3+1}{6+3}=\frac{4}{9},P(X^{(2)}=M|Y=-1)=\frac{2+1}{6+3}=\frac{3}{9},P(X^{(2)}=L|Y=-1)=\frac{1+1}{6+3}=\frac{2}{9}$$

对于给定的$x=(2,S)^T$计算，

$$P(Y=1)P(X^{(1)}=2|Y=1)P(X^{(2)}=S|Y=1)=\frac{10}{17} \times\frac{4}{12} \times \frac{2}{12}=0.0327$$
$$P(Y=-1)P(X^{(1)}=2|Y=-1)P(X^{(2)}=S|Y=-1)=\frac{7}{17} \times\frac{3}{9} \times \frac{4}{9}=0.0610$$
由于$P(Y=-1)P(X^{(1)}=2|Y=-1)P(X^{(2)}=S|Y=-1)$最大，所以$y=-1$

代码

# coding=utf-8
from numpy import *
import os


def createVocabList(dataSet):
    vocabSet = set()
    for document in dataSet:
        vocabSet = vocabSet | set(document)
    return list(vocabSet)


# 制作词袋模型
def setOfWords2Vec(vocabList, inputSet):
    returnVec = [0] * len(vocabList)
    for word in inputSet:
        if word in vocabList:
            returnVec[vocabList.index(word)] = 1
        else:
            print 'the word {} is not in my Vocabulary!'.format(word)
    return returnVec


def trainNB0(trainMatrix, trainCategory):
    numTrainDocs = len(trainMatrix)
    numWords = len(trainMatrix[0])
    # 包含敏感词的此条出现的概率
    pAbusive = sum(trainCategory) / float(numTrainDocs)
    # 构造单词出现次数列表 [1,1,...,1],这里使用了拉普拉斯平滑的思想
    p0Num = ones(numWords)
    p1Num = ones(numWords)
    p0Denom = 2.0
    p1Denom = 2.0
    for i in range(numTrainDocs):
        # 是否是侮辱性文件
        if trainCategory[i] == 1:
            # 对侮辱性文件的向量进行加和,[0,1,1,....] + [0,1,1,....]->[0,2,2,...]
            p1Num += trainMatrix[i]
            # 计算所有侮辱性文件中出现的单词总数
            p1Denom += sum(trainMatrix[i])
        else:
            p0Num += trainMatrix[i]
            p0Denom += sum(trainMatrix[i])
    # 在1类别下，每个单词出现的概率,[1,2,3,5]/90->[1/90,...]
    p1Vect = log(p1Num / p1Denom)
    #  在0类别下，每个单词出现的概率
    p0Vect = log(p0Num / p0Denom)
    return p0Vect, p1Vect, pAbusive


def classifyNB(vec2Classify, p0Vec, p1Vec, pClass1):
    # p1Vec是每个词在类别1里出现的概率，vec2Classify表示这个词在该文章中是否出现过，乘起来表示文章中每个词出现的概率
    # 这里进行了log变换，所以写成了加法
    p1 = sum(vec2Classify * p1Vec) + log(pClass1)
    p0 = sum(vec2Classify * p0Vec) + log(1 - pClass1)
    if p1 > p0:
        return 1
    else:
        return 0


def testParse(bigString):
    import re
    listOfTokens = re.split('r\W*', bigString)
    return [tok.lower() for tok in listOfTokens if len(tok) > 2]


def spamTest():
    docList = []
    classList = []
    fullText = []
    for i in range(1, 26):
        # 读取非垃圾邮件文本
        wordList = testParse(open(os.getcwd() + '/email/spam/{}.txt'.format(i)).read())
        docList.append(wordList)
        fullText.append(wordList)
        classList.append(1)
        # 读取垃圾邮件文本
        wordList = testParse(open(os.getcwd() + '/email/ham/{}.txt'.format(i)).read())
        docList.append(wordList)
        fullText.append(wordList)
        classList.append(0)
    # 建立词集
    vocabList = createVocabList(docList)
    # 训练集语料的下标
    trainingSet = range(50)
    testSet = []
    # 选10个语料作为测试集
    for i in range(10):
        randIndex = int(random.uniform(0, len(trainingSet)))
        testSet.append(trainingSet[randIndex])
        # 删除后训练集变成40个
        del (trainingSet[randIndex])
    trainMat = []
    trainClasses = []
    for docIndex in trainingSet:
        trainMat.append(setOfWords2Vec(vocabList, docList[docIndex]))
        trainClasses.append(classList[docIndex])
    # 计算各种概率
    p0V, p1V, pSpam = trainNB0(array(trainMat), array(trainClasses))
    errorCount = 0
    # 进行测试
    for docIndex in testSet:
        # 构建词向量
        wordVector = setOfWords2Vec(vocabList, docList[docIndex])
        # 进行分类
        if classifyNB(array(wordVector), p0V, p1V, pSpam) != classList[docIndex]:
            errorCount += 1
    print 'the error rate is {}'.format(float(errorCount) / len(testSet))


def main():
    spamTest()


if __name__ == '__main__':
    main()

代码详见朴素贝叶斯代码实现

参考

1. 机器学习实战，人民邮电出版社
2. 统计学习方法，清华大学出版社
3. AI Learning