#朴素贝叶斯的特点
##朴素贝叶斯的用途
##朴素贝叶斯的适用范围
###数据类型
离散数据和连续数据都可用,前者是用频数的那一套方法,后者是用贝叶斯估计的那套方法
###朴素贝叶斯的优点
-
对缺失数据不敏感,在数据较少的情况下依然可以使用该方法
-
可以处理多个类别 的分类问题
-
适用于标称型数据
###朴素贝叶斯的局限性 -
监督学习,需要确定分类的目标(分类中)
-
对输入数据的形式比较敏感
-
由于用先验数据去预测分类,因此存在误差
-
运算量很大,在样本很多,变量很多的时候学习速度慢,而且朴素贝叶斯也是一种随性学习方法
#朴素贝叶斯分类器原理
本部分内容来源《数据挖掘导论》4.4节
贝叶斯分类器是基于贝叶斯定理构建出来的分类器,是一个统计分类器。对分类方法进行比较的有关研究表明,简单贝叶斯分类器在分类性能上与决策树和神经网络都是可比拟的。在处理大规模数据时,贝叶斯分类器已经表现出较高的准确性和运算性能。
简单贝叶斯分类器的假设是:指定类别中,样本各特征之间相互独立,即某一特征不同取值的概率与其他特征的取值没有任何关系。
假定:每个数据样本均是由一个n维特征向量
X
=
x
1
,
x
2
,
…
,
x
n
X={x_1,x_2,…,x_n}
X=x1,x2,…,xn来表示与属性
(
A
1
,
A
2
,
.
.
.
,
A
n
)
(A_1,A_2,...,A_n)
(A1,A2,...,An)对应的值,所有的样本共有m个类别,
C
1
,
C
2
,
.
.
.
,
C
m
C_1,C_2,...,C_m
C1,C2,...,Cm。对于任一个样本,我们要预测其属于哪一类样本。
我们的目标是计算
P
(
C
i
/
X
)
P(C_i/X)
P(Ci/X),式中
C
i
C_i
Ci表示某类样本,
P
(
C
i
/
X
)
P(C_i/X)
P(Ci/X)表示已知样本X的基础上样本为
C
i
C_i
Ci假设发生的概率,比如X={“红色”,“圆形”},
C
i
C_i
Ci:这是一个苹果,
P
(
C
i
/
X
)
P(C_i/X)
P(Ci/X)则表示在已知某物品是红色和圆形的前提下这是一个苹果的概率;在分类问题中,将
P
(
C
i
/
X
)
P(C_i/X)
P(Ci/X)最大值对应的
C
i
C_i
Ci作为样本X的最终类别。
如何求取
m
a
x
(
P
(
C
i
/
X
)
)
max(P(C_i/X))
max(P(Ci/X))?
根据贝叶斯定理,
P
(
C
i
/
X
)
=
P
(
X
/
C
i
)
∗
P
(
C
i
)
/
P
(
X
)
P(C_i/X)=P(X/C_i)*P(C_i)/P(X)
P(Ci/X)=P(X/Ci)∗P(Ci)/P(X),这也是此类分类器被称为贝叶斯分类器的原因。
上式中,对于给定样本X,P(X)是定值,故
m
a
x
(
P
(
C
i
/
X
)
)
max(P(C_i/X))
max(P(Ci/X))等价于求
m
a
x
(
P
(
X
/
C
i
)
∗
P
(
C
i
)
)
max(P(X/C_i)*P(C_i))
max(P(X/Ci)∗P(Ci)),对于给定训练样本集,
P
(
C
i
)
P(C_i)
P(Ci)可利用统计方法求得,而
P
(
X
/
C
i
)
P(X/C_i)
P(X/Ci)可在类别
C
i
C_i
Ci中求X的概率,但是在样本量较大、特征数目较多的时候,该方法求X会带来维数灾难,而且也不利于将这种算法拓展到分布式实现。这里通过引入贝叶斯假设:在同一类别中,不同特征之间相互独立。该假设的引入简化了计算
P
(
X
/
C
i
)
=
P
(
x
1
/
C
i
)
∗
P
(
x
2
/
C
i
)
∗
…
∗
P
(
x
n
/
C
i
)
P(X/C_i)=P(x_1/C_i)*P(x_2/C_i)*…*P(x_n/C_i)
P(X/Ci)=P(x1/Ci)∗P(x2/Ci)∗…∗P(xn/Ci)。即一个样本在某一类样本出现的概率就变成了计算该样本各特征的值在该类样本中出现概率。
##计算步骤
###计算类别
C
i
C_i
Ci的先验概率
P
(
C
i
)
=
s
i
s
P(C_i)=\frac{s_i}{s}
P(Ci)=ssi
式中
s
i
s_i
si表示类别为
C
i
C_i
Ci的样本数目,s表示样本总体数目。
###计算类别
C
i
C_i
Ci中属性
A
k
A_k
Ak取值为
v
j
v_j
vj的概率
一、若属性
A
k
A_k
Ak是符号变量,则:
P
(
x
k
/
C
i
)
=
s
i
k
s
i
P(x_k/C_i)=\frac{s_{ik}}{s_i}
P(xk/Ci)=sisik
式中
s
i
k
s_{ik}
sik表示类别
C
i
C_i
Ci的样本中属性
A
k
A_k
Ak取值为
v
j
v_j
vj的样本数目,
s
i
s_i
si表示类别
C
i
C_i
Ci的样本数。
二、若属性
A
k
A_k
Ak是连续变量,在假设属性具有高斯分布的情况下,有:
P
(
x
k
/
C
i
)
=
g
(
x
k
,
μ
C
i
,
ρ
C
i
)
=
1
2
π
ρ
C
i
e
−
(
x
−
μ
C
i
)
2
2
ρ
C
i
2
P(x_k/C_i)=g(x_k,\mu_{C_i},\rho_{C_i})=\frac{1}{\sqrt{2\pi}\rho_{C_i}}e^{-\frac{(x-\mu_{C_i})^2}{2\rho_{C_i}^2}}
P(xk/Ci)=g(xk,μCi,ρCi)=2πρCi1e−2ρCi2(x−μCi)2
式中
μ
C
i
\mu_{C_i}
μCi和
ρ
C
i
\rho_{C_i}
ρCi分别表示类别
C
i
C_i
Ci的样本中属性
A
k
A_k
Ak的均值和方差。
###计算
P
(
X
/
C
i
)
P(X/C_i)
P(X/Ci)
P
(
X
/
C
i
)
=
∏
k
=
1
k
P
(
x
k
/
C
i
)
P(X/C_i)=\prod_{k=1}^kP(x_k/C_i)
P(X/Ci)=∏k=1kP(xk/Ci)
对每一类样本都计算
P
(
X
/
C
i
)
P(X/C_i)
P(X/Ci)
###确定样本类别
样本X归属类别
C
I
C_I
CI当且仅当:
P
(
X
/
C
I
)
>
P
(
X
/
C
j
)
P(X/C_I)>P(X/C_j)
P(X/CI)>P(X/Cj)
式中,
1
≤
I
≤
m
,
1
≤
j
≤
m
1\leq I\leq m,1\leq j \leq m
1≤I≤m,1≤j≤m
贝叶斯算法流程图
从理论上讲与其它分类器相比,贝叶斯分类器具有最小的错误率。但实际上由于其所依据的类别独立性假设和缺乏某些数据的准确概率分布,从而使得贝叶斯分类器预测准确率受到影响。
##算法实现
算法来源网络,没找到具体链接,可以参考该文。
数据存在本地,文件名为"pima-indians-diabetes.data.csv"。
数据大致如下,共768行9列,最后一列为标签列:
6 148 72 35 0 33.6 0.627 50 1
1 85 66 29 0 26.6 0.351 31 0
8 183 64 0 0 23.3 0.672 32 1
1 89 66 23 94 28.1 0.167 21 0
0 137 40 35 168 43.1 2.288 33 1
程序代码如下:
# -*- coding: UTF-8 -*-
import csv
import random
import math
import numpy as np
'''
1.处理数据
首先加载数据文件。CSV格式的数据没有标题行和任何引号。我们可以使用csv模块中的open函数打开文件,使用reader函数读取行数据。
我们也需要将以字符串类型加载进来属性转换为我们可以使用的数字。下面是用来加载匹马印第安人数据集(Pima indians dataset)的loadCsv()函数。
'''
'''
1.处理数据
首先加载数据文件。CSV格式的数据没有标题行和任何引号。我们可以使用csv模块中的open函数打开文件,使用reader函数读取行数据。
我们也需要将以字符串类型加载进来属性转换为我们可以使用的数字。下面是用来加载匹马印第安人数据集(Pima indians dataset)的loadCsv()函数。
'''
def loadCsv(filename):
lines = csv.reader(open(filename, encoding='iso-8859-1'))
dataset = list(lines)
for i in range(len(dataset)):
dataset[i] = [float(x) for x in dataset[i]]#将字符串类型转换成浮点型数据
return dataset
##我们可以通过加载皮马印第安人数据集,然后打印出数据样本的个数,以此测试这个函数。
#filename = './pima-indians-diabetes.data.csv'
#dataset = loadCsv(filename)
#print(('Loaded data file {0} with {1} rows').format(filename, len(dataset)))
#打印结果如下
#'Loaded data file ./pima-indians-diabetes.data.csv with 768 rows'
'''
下一步,我们将数据分为用于朴素贝叶斯预测的训练数据集,以及用来评估模型精度的测试数据集。我们需要将数据集随机分为包含67%的训练集合和包含33%的测试集(这是在此数据集上测试算法的通常比率)。
下面是splitDataset()函数,它以给定的划分比例将数据集进行划分。
dataset为数据集,splitRatio为分裂的比例即训练集合在总数据中的比例
'''
def splitDataset(dataset, splitRatio):
trainSize = int(len(dataset) * splitRatio)
trainSet = []
copy = list(dataset)
while len(trainSet) < trainSize:
index = random.randrange(len(copy))
trainSet.append(copy.pop(index))#一下子得到两个数据集
return [trainSet, copy]
##我们可以定义一个具有5个样例的数据集来进行测试,首先它分为训练数据集和测试数据集,然后打印出来,看看每个数据样本最终落在哪个数据集。
#dataset = [[1], [2], [3], [4], [5]]
#splitRatio = 0.67
#train, test = splitDataset(dataset, splitRatio)
#print(('Split {0} rows into train with {1} and test with {2}').format(len(dataset), train, test))
#打印结果如下:
#'Split 5 rows into train with [[5], [2], [3]] and test with [[1], [4]]'
'''
提取数据特征
朴素贝叶斯模型包含训练数据集中数据的特征,然后使用这个数据特征来做预测。
所收集的训练数据的特征,包含相对于每个类的每个属性的均值和标准差。举例来说,如果如果有2个类和7个数值属性,然后我们需要每一个属性(7)和类(2)的组合的均值和标准差,也就是14个属性特征。
也就是求取特征的概率分布
在对特定的属性归属于每个类的概率做计算、预测时,将用到这些特征。
我们将数据特征的获取划分为以下的子任务:
1 按类别划分数据
2 计算均值
3 计算标准差
4 提取数据集特征
5 按类别提取属性特征
'''
'''
按类别划分数据
首先将训练数据集中的样本按照类别进行划分,然后计算出每个类的统计数据。我们可以创建一个类别到属于此类别的样本列表的的映射,并将整个数据集中的样本分类到相应的列表。
下面的SeparateByClass()函数可以完成这个任务:
'''
def separateByClass(dataset):
separated = {}
for i in range(len(dataset)):
vector = dataset[i]
if (vector[-1] not in separated):
separated[vector[-1]] = []#求到不同的类
separated[vector[-1]].append(vector)#将某类样本添加到该类的列表中
return separated
##可以看出,函数假设样本中最后一个属性(-1)为类别值,返回一个类别值到数据样本列表的映射。
##我们可以用一些样本数据测试如下:
#dataset = [[1,20,1], [2,21,0], [3,22,1]]
#separated = separateByClass(dataset)#得到一字典,以类别标签为key,以样本为元素组成的列表作为值
#print(('Separated instances: {0}').format(separated))
#测试结果如下:
#Separated instances: {1: [[1, 20, 1], [3, 22, 1]], 0: [[2, 21, 0]]}
'''
计算均值、方差
我们需要计算在每个类中每个属性的均值。均值是数据的中点或者集中趋势,在计算概率时,我们用它作为高斯分布的中值。
我们也需要计算每个类中每个属性的标准差。标准差描述了数据散布的偏差,在计算概率时,我们用它来刻画高斯分布中,每个属性所期望的散布。
标准差是方差的平方根。方差是每个属性值与均值的离差平方的平均数。注意我们使用N-1的方法(译者注:参见无偏估计),也就是在在计算方差时,属性值的个数减1。
'''
def mean(numbers):
return sum(numbers)/float(len(numbers))
def stdev(numbers):
avg = mean(numbers)
variance = sum([pow(x-avg,2) for x in numbers])/float(len(numbers)-1)
return math.sqrt(variance)
##通过计算从1到5这5个数的均值来测试函数。
#numbers = [1,2,3,4,5]
#print(('Summary of {0}: mean={1}, stdev={2}').format(numbers, mean(numbers), stdev(numbers)))
#测试结果如下:
#Summary of [1, 2, 3, 4, 5]: mean=3.0, stdev=1.5811388300841898
#其实没有必要自己来定义这些基础的函数,虽然math里没有但numpy里就有这些函数
'''
提取数据集的特征
现在我们可以提取数据集特征。对于一个给定的样本列表(对应于某个类),我们可以计算每个属性的均值和标准差。
zip函数将数据样本按照属性分组为一个个列表,然后可以对每个属性计算均值和标准差。
'''
def summarize(dataset):#几个属性就得到几个均值、方差对,形式为[(2.0, 1.0), ... ,(21.0, 1.0)]
summaries = [(mean(attribute), stdev(attribute)) for attribute in zip(*dataset)]
del summaries[-1]
return summaries
##我们可以使用一些测试数据来测试这个summarize()函数,测试数据对于第一个和第二个数据属性的均值和标准差显示出显著的不同。
#dataset = [[1,20,0], [2,21,1], [3,22,0]]
#summary = summarize(dataset)
#print(('Attribute summaries: {0}').format(summary))
#测试结果如下:
#Attribute summaries: [(2.0, 1.0), (21.0, 1.0)]
'''
按类别提取属性特征
合并代码,我们首先将训练数据集按照类别进行划分,然后计算每个属性的摘要。
'''
def summarizeByClass(dataset):#不同类别及各属性的均值方差,形式为{1: [(2.0, 1.4142135623730951), ... ,(21.0, 1.4142135623730951)], 0: [(3.0, 1.4142135623730951), ... ,(21.5, 0.7071067811865476)]}
separated = separateByClass(dataset)
summaries = {}
for classValue, instances in separated.items():
summaries[classValue] = summarize(instances)
return summaries
##使用小的测试数据集来测试summarizeByClass()函数。
#dataset = [[1,20,1], [2,21,0], [3,22,1], [4,22,0]]
#summary = summarizeByClass(dataset)
#print(('Summary by class value: {0}').format(summary))
#测试结果如下:
#Summary by class value: {1: [(2.0, 1.4142135623730951), (21.0, 1.4142135623730951)], 0: [(3.0, 1.4142135623730951), (21.5, 0.7071067811865476)]}
'''
预测
我们现在可以使用从训练数据中得到的摘要来做预测。做预测涉及到对于给定的数据样本,计算其归属于每个类的概率,然后选择具有最大概率的类作为预测结果。
我们可以将这部分划分成以下任务:
1 计算高斯概率密度函数
2 计算对应类的概率
3 单一预测
4 评估精度
'''
'''
计算高斯概率密度函数
给定来自训练数据中已知属性的均值和标准差,我们可以使用高斯函数来评估一个给定的属性值的概率。
已知每个属性和类值的属性特征,在给定类值的条件下,可以得到给定属性值的条件概率。
关于高斯概率密度函数,可以查看参考文献。总之,我们要把已知的细节融入到高斯函数(属性值,均值,标准差),并得到属性值归属于某个类的似然(译者注:即可能性)。
在calculateProbability()函数中,我们首先计算指数部分,然后计算等式的主干。这样可以将其很好地组织成2行。
'''
import math
def calculateProbability(x, mean, stdev):
exponent = math.exp(-(math.pow(x-mean,2)/(2*math.pow(stdev,2))))
return (1 / (math.sqrt(2*math.pi) * stdev)) * exponent
##使用一些简单的数据测试如下:
#x = 71.5
#mean = 73
#stdev = 6.2
#probability = calculateProbability(x, mean, stdev)
#print(('Probability of belonging to this class: {0}').format(probability))
#测试结果如下:
#Probability of belonging to this class: 0.06248965759370005
'''
计算所属类的概率
既然我们可以计算一个属性属于某个类的概率,那么合并一个数据样本中所有属性的概率,最后便得到整个数据样本属于某个类的概率。
使用乘法合并概率,在下面的calculClassProbilities()函数中,给定一个数据样本,它所属每个类别的概率,可以通过将其属性概率相乘得到。结果是一个类值到概率的映射。
'''
def calculateClassProbabilities(summaries, inputVector):
probabilities = {}
for classValue, classSummaries in summaries.items():#逐个类别内部计算某样本属于某一类的概率
probabilities[classValue] = 1#初始化样本属于某类的概率为1
for i in range(len(classSummaries)):#逐个属性计算某样本属于某一类的概率
mean, stdev = classSummaries[i]
x = inputVector[i]
probabilities[classValue] *= calculateProbability(x, mean, stdev)#连乘,会直接改变probabilities[classValue]的值
return probabilities
##测试calculateClassProbabilities()函数。
#summaries = {0:[(1, 0.5)], 1:[(20, 5.0)]}
#inputVector = [1.1, '?']
#probabilities = calculateClassProbabilities(summaries, inputVector)
#print(('Probabilities for each class: {0}').format(probabilities))
#测试结果如下:
#Probabilities for each class: {0: 0.7820853879509118, 1: 6.298736258150442e-05}
'''
单一预测
既然可以计算一个数据样本属于每个类的概率,那么我们可以找到最大的概率值,并返回关联的类。
下面的predict()函数可以完成以上任务。
'''
def predict(summaries, inputVector):
probabilities = calculateClassProbabilities(summaries, inputVector)
bestLabel, bestProb = None, -1
for classValue, probability in probabilities.items():
if bestLabel is None or probability > bestProb:
bestProb = probability
bestLabel = classValue
return bestLabel
##测试predict()函数如下:
#summaries = {'A':[(1, 0.5)], 'B':[(20, 5.0)]}
#inputVector = [1.1, '?']
#result = predict(summaries, inputVector)
#print(('Prediction: {0}').format(result))
#测试结果如下:
#Prediction: A
'''
多重预测
最后,通过对测试数据集中每个数据样本的预测,我们可以评估模型精度。getPredictions()函数可以实现这个功能,并返回每个测试样本的预测列表。
'''
def getPredictions(summaries, testSet):
predictions = []
for i in range(len(testSet)):
result = predict(summaries, testSet[i])
predictions.append(result)
return predictions
##测试getPredictions()函数如下。
#summaries = {'A':[(1, 0.5)], 'B':[(20, 5.0)]}
#testSet = [[1.1, '?'], [19.1, '?']]
#predictions = getPredictions(summaries, testSet)
#print(('Predictions: {0}').format(predictions))
##测试结果如下:
##Predictions: ['A', 'B']
'''
计算精度
预测值和测试数据集中的类别值进行比较,可以计算得到一个介于0%~100%精确率作为分类的精确度。getAccuracy()函数可以计算出这个精确率。
'''
def getAccuracy(testSet, predictions):
correct = 0
for x in range(len(testSet)):
if testSet[x][-1] == predictions[x]:
correct += 1
return (correct/float(len(testSet))) * 100.0
##我们可以使用如下简单的代码来测试getAccuracy()函数。
#testSet = [[1,1,1,'a'], [2,2,2,'a'], [3,3,3,'b']]
#predictions = ['a', 'a', 'a']
#accuracy = getAccuracy(testSet, predictions)
#print(('Accuracy: {0}').format(accuracy))
#测试结果如下:
#Accuracy: 66.66666666666666
#以上就完成了相应的功能模块,下面用这些模块来完成一个预测
def main():
filename = './pima-indians-diabetes.data.csv'
splitRatio = 0.67
dataset = loadCsv(filename)
trainingSet, testSet = splitDataset(dataset, splitRatio)
print(('Split {0} rows into train={1} and test={2} rows').format(len(dataset),len(trainingSet),len(testSet)))
# prepare model
summaries = summarizeByClass(trainingSet)
# test model
predictions = getPredictions(summaries, testSet)
accuracy = getAccuracy(testSet, predictions)
print(('Accuracy: {0}%').format(accuracy))
main()
#组合后测试结果如下:
#Split 768 rows into train=514 and test=254 rows
#Accuracy: 75.59055118110236%
也可以参考下面的说明,差不多。
#朴素贝叶斯原理
来源:
朴素贝叶斯法是基于贝叶斯定理与特征条件独立假设的分类方法。训练的时候,学习输入输出的联合概率分布;分类的时候,利用贝叶斯定理计算后验概率最大的输出。
##朴素贝叶斯法的学习与分类
基本方法
设输入空间
χ
⊆
R
n
\chi \subseteq R^n
χ⊆Rn,为n维向量的集合,输出空间为类标记集合
y
=
{
c
1
…
…
c
k
}
y=\{c_1……c_k\}
y={c1……ck}。输入特征向量
x
x
x和输出类标记
y
y
y分属于这两个集合。
X
X
X是输入空间上的随机变量,
Y
Y
Y是输出空间上的随机变量。
P
(
X
,
Y
)
P(X,Y)
P(X,Y)是
X
X
X和
Y
Y
Y的联合概率分布,训练数据集
T
=
{
(
x
1
,
y
1
)
,
(
x
2
,
y
2
)
,
…
,
(
x
N
,
y
N
)
}
T=\{(x_1,y_1),(x_2,y_2),\dots,(x_N,y_N)\}
T={(x1,y1),(x2,y2),…,(xN,yN)}
由P(X,Y)独立同分布产生。
朴素贝叶斯法通过T学习联合概率分布P(X,Y)。具体来讲,学习以下先验概率:
P
(
Y
c
k
)
,
k
=
1
,
2
,
…
,
K
P(Y_{c_k}),k=1,2,\dots,K
P(Yck),k=1,2,…,K
以及条件概率分布:
P
(
X
=
x
∣
Y
=
c
k
)
=
P
(
X
(
1
)
=
x
(
1
)
,
…
,
X
(
n
)
=
x
(
n
)
∣
Y
=
c
k
)
,
k
=
1
,
2
,
…
,
K
P(X=x|Y=c_k)=P(X^{(1)}=x^{(1)},\dots,X^{(n)}=x^{(n)}|Y=c_k),k=1,2,\dots,K
P(X=x∣Y=ck)=P(X(1)=x(1),…,X(n)=x(n)∣Y=ck),k=1,2,…,K
于是根据联合概率密度分布函数:
P
(
X
=
x
,
Y
=
y
)
=
P
(
X
=
x
)
P
(
Y
=
y
∣
X
=
x
)
=
P
(
Y
=
y
)
P
(
X
=
x
∣
Y
=
y
)
P(X=x,Y=y)=P(X=x)P(Y=y|X=x)=P(Y=y)P(X=x|Y=y)
P(X=x,Y=y)=P(X=x)P(Y=y∣X=x)=P(Y=y)P(X=x∣Y=y)
从而有
P
(
Y
=
y
∣
X
=
x
)
=
P
(
Y
=
y
)
P
(
X
=
x
∣
Y
=
y
)
P
(
X
=
x
)
P(Y=y|X=x)=\frac{P(Y=y)P(X=x|Y=y)}{P(X=x)}
P(Y=y∣X=x)=P(X=x)P(Y=y)P(X=x∣Y=y)
上面的式子就是我们的求解目标,也就是在
X
=
x
X=x
X=x的前提下不同
y
y
y的取值概率,在分类里,我们可以概率最大的y当作样本
x
x
x的类别。
维数灾难
条件概率分布
P
(
X
=
x
∣
Y
=
c
k
)
P(X=x|Y=c_k)
P(X=x∣Y=ck)的参数数量是指数级的,也就是X和Y的组合很多,假设
x
j
x_j
xj可能取值
S
j
S_j
Sj个,
Y
Y
Y可能取值有
K
K
K个,那么参数的个数是
K
∏
i
=
1
n
S
j
K\prod_{i=1}^n S_j
K∏i=1nSj 。特别地,取
x
j
=
S
x_j=S
xj=S,那么参数个数为
K
S
n
K S^n
KSn,当维数
n
n
n很大的时候,就会发生维数灾难。
假如一个特征的取值有
1
0
2
10^2
102种可能,而在10维度空间中,就
1
0
20
10^{20}
1020种可能。计算方式用Python描述如下:
dimensionality = 10
print(1 / (0.01 ** dimensionality))
也可以如下可视化:
# -*- coding:utf-8 -*-
# Filename: dimensionality.py
# Author:hankcs
# Date: 2015/2/6 14:40
from matplotlib import pyplot as plt
import numpy as np
max_dimensionality = 10
ax = plt.axes(xlim=(0, max_dimensionality), ylim=(0, 1 / (0.01 ** max_dimensionality)))
x = np.linspace(0, max_dimensionality, 1000)
y = 1 / (0.01 ** x)
plt.plot(x, y, lw=2)
plt.show()
可视化图像:
这种指数级的复杂度增长被称为维数灾难。
看完这个图大概就能理解为什么条件概率分布按原本的定义无法计算了。
为了计算它,朴素贝叶斯法对它做了条件独立性的假设:
P
(
X
(
1
)
=
x
(
1
)
,
…
,
X
(
n
)
=
x
(
n
)
∣
Y
=
c
k
)
=
∏
j
=
1
n
P
(
X
(
j
)
=
x
(
j
)
∣
Y
=
c
k
)
P(X^{(1)}=x^{(1)},\dots,X^{(n)}=x^{(n)}|Y=c_k)=\prod_{j=1}^n P(X^{(j)}=x^{(j)}|Y=c_k)
P(X(1)=x(1),…,X(n)=x(n)∣Y=ck)=j=1∏nP(X(j)=x(j)∣Y=ck)
也就是说样本某一维特征的取值与其他特征的取值没有任务关系,不会相互影响。这一假设简化了计算,也牺牲了一定的分类准确率。
基于此假设,以及贝叶斯定理,后验概率为:
P
(
Y
=
c
k
∣
X
=
x
)
=
P
(
X
=
x
∣
Y
=
c
k
)
P
(
Y
=
c
k
)
P
(
X
=
x
)
=
P
(
X
=
x
∣
Y
=
c
k
)
P
(
Y
=
c
k
)
∑
k
P
(
X
=
x
∣
Y
=
c
k
)
P
(
Y
=
c
k
)
P(Y=c_k|X=x)=\frac{P(X=x|Y=c_k)P(Y=c_k)}{P(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=ck∣X=x)=P(X=x)P(X=x∣Y=ck)P(Y=ck)=∑kP(X=x∣Y=ck)P(Y=ck)P(X=x∣Y=ck)P(Y=ck)
将独立性假设代入上式,得到
P
(
Y
=
c
k
∣
X
=
x
)
=
P
(
Y
=
c
k
)
∏
i
=
1
n
P
(
X
(
i
)
=
x
(
i
)
∣
Y
=
c
k
)
P
(
X
=
x
)
P(Y=c_k|X=x)=\frac{P(Y=c_k)\prod_{i=1}^n P(X^{(i)}=x^{(i)}|Y=c_k)}{P(X=x)}
P(Y=ck∣X=x)=P(X=x)P(Y=ck)∏i=1nP(X(i)=x(i)∣Y=ck)
朴素贝叶斯分类器可以表示为:
y
=
f
(
x
)
=
a
r
g
m
a
x
c
k
P
(
Y
=
c
k
)
∏
i
=
1
n
P
(
X
(
i
)
=
x
(
i
)
∣
Y
=
c
k
)
P
(
X
=
x
)
y=f(x)=arg~ max_{c_k} \frac{P(Y=c_k)\prod_{i=1}^n P(X^{(i)}=x^{(i)}|Y=c_k)}{P(X=x)}
y=f(x)=arg maxckP(X=x)P(Y=ck)∏i=1nP(X(i)=x(i)∣Y=ck)
事实上,对于给定的样本,等式右边的分母值是固定的,跟
c
k
c_k
ck一点关系都没有,所以分母可以去掉,得到:
y
=
a
r
g
m
a
x
c
k
P
(
Y
=
c
k
)
∏
i
=
1
n
P
(
X
(
i
)
=
x
(
i
)
∣
Y
=
c
k
)
y=arg~ max_{c_k} P(Y=c_k)\prod_{i=1}^n P(X^{(i)}=x^{(i)}|Y=c_k)
y=arg maxckP(Y=ck)i=1∏nP(X(i)=x(i)∣Y=ck)
后验概率最大化的含义
参原文
##朴素贝叶斯法的参数估计
###极大似然估计
朴素贝叶斯法要学习的东西就是
P
(
Y
=
c
k
)
P(Y=c_k)
P(Y=ck)和
P
(
X
=
x
∣
Y
=
c
k
)
P(X=x|Y=c_k)
P(X=x∣Y=ck),这两个概率的估计用极大似然估计法(简单讲,就是用样本猜测模型参数,或者说使得似然函数最大的参数)进行:
P
(
Y
=
c
k
)
=
∑
i
=
1
N
I
(
y
(
i
)
=
c
k
)
N
,
k
=
1
,
…
,
K
P(Y=c_k)=\frac{\sum_{i=1}^N I(y^{(i)}=c_k)}{N},k=1,\dots,K
P(Y=ck)=N∑i=1NI(y(i)=ck),k=1,…,K
函数I用来判断样本是否属于类别
c
k
c_k
ck
也就是用样本中
c
k
c_k
ck的出现次数除以样本容量。
P
(
X
(
i
)
=
x
(
i
j
)
∣
Y
=
c
k
)
=
∑
i
=
1
N
I
(
X
(
i
)
=
x
(
i
j
)
,
y
(
i
)
=
c
k
)
∑
i
=
1
N
I
(
y
(
i
)
=
c
k
)
P(X^{(i)}=x^{(ij)}|Y=c_k)=\frac{\sum_{i=1}^N I(X^{(i)}=x^{(ij)},y^{(i)}=c_k)}{\sum_{i=1}^N I(y^{(i)}=c_k)}
P(X(i)=x(ij)∣Y=ck)=∑i=1NI(y(i)=ck)∑i=1NI(X(i)=x(ij),y(i)=ck)
x
(
i
j
)
x^{(ij)}
x(ij)表示第
i
i
i个样本的第j种取值。
分子是样本中变量组合的出现次数,分母是上面说过的样本中
c
k
c_k
ck的出现次数。
###学习与分类算法
于是就有朴素贝叶斯算法,先从训练数据中计算先验概率和条件概率,然后对于给定的实例计算最大的条件概率,输出该条件对应的类别。形式化的描述如下:
贝叶斯估计
最大似然估计有个隐患,假设训练数据中没有出现某种参数和类别的组合怎么办?此时估计的概率值为0,但是这不代表真实数据中就没有这样的组合。解决办法是采用贝叶斯估计
1、条件概率的贝叶斯估计:
P
(
X
(
i
)
=
x
(
i
j
)
∣
Y
=
c
k
)
=
∑
i
=
1
N
I
(
X
(
i
)
=
x
(
i
j
)
,
y
(
i
)
=
c
k
)
+
λ
∑
i
=
1
N
I
(
y
(
i
)
=
c
k
)
+
S
j
λ
P(X^{(i)}=x^{(ij)}|Y=c_k)=\frac{\sum_{i=1}^N I(X^{(i)}=x^{(ij)},y^{(i)}=c_k)+\lambda}{\sum_{i=1}^N I(y^{(i)}=c_k)+S_j \lambda}
P(X(i)=x(ij)∣Y=ck)=∑i=1NI(y(i)=ck)+Sjλ∑i=1NI(X(i)=x(ij),y(i)=ck)+λ
其中
λ
≥
0
\lambda\geq 0
λ≥0,
S
j
S_j
Sj表示
x
(
i
)
x^{(i)}
x(i)可能的取值种数。分子和分母分别比最大似然估计多了一点东西,其意义是在随机变量每个取值的频数上加一个常量
λ
\lambda
λ。当此常量取0时,就是最大似然估计,当此常量取1时,称为拉普拉斯平滑。
2、先验概率的贝叶斯估计:
P ( Y = c k ) = ∑ i = 1 N I ( y ( i ) = c k ) + λ N + K λ , k = 1 , … , K P(Y=c_k)=\frac{\sum_{i=1}^N I(y^{(i)}=c_k)+\lambda}{N+K\lambda},k=1,\dots,K P(Y=ck)=N+Kλ∑i=1NI(y(i)=ck)+λ,k=1,…,K
#贝叶斯多分类
贝叶斯天然就是一个多分类器
#贝叶斯回归
贝叶斯线性回归:思路还是线性回归的思路
y
(
x
,
w
)
=
w
T
ϕ
(
x
)
y(x,w)=w^T\phi (x)
y(x,w)=wTϕ(x)
然后用贝叶斯估计去求解参数
w
w
w,而传统的方法是用最小二乘法。
而贝叶斯回归里有一类特殊的也是出现最多,当先验分布为
p
(
w
)
=
N
(
w
∣
0
,
α
−
1
I
)
p(w)=N(w|0,α^{−1}I)
p(w)=N(w∣0,α−1I)时,这就是贝叶斯岭回归。
https://blog.csdn.net/daunxx/article/details/51725086
https://blog.csdn.net/u010016927/article/details/75000152
https://blog.csdn.net/u010016927/article/details/75000207
http://lib.csdn.net/article/machinelearning/37615
#朴素贝叶斯模型的检验(及自变量的选择)
#朴素贝叶斯模型的评估方法
#朴素贝叶斯的正则化
#并行朴素贝叶斯
http://www.aboutyun.com/thread-20839-1-1.html
#朴素贝叶斯的python实现
其他参考资料:
李航《统计学习方法》第四章——用Python实现朴素贝叶斯分类器(MNIST数据集):里面有多种算法的python实现。
机器学习实战笔记(Python实现)-03-朴素贝叶斯
机器学习之朴素贝叶斯(NB)分类算法与Python实现
使用Spark进行数据挖掘之实现朴素贝叶斯算法
来自《机器学习》第7章
7.5 贝叶斯网
贝叶斯网是一种经典的概率图模型。
为了简化讨论,本节假定所有属性均为离数型。对于连续属性,条件概率表可推广为条件概率密度函数。
贝叶斯网( Bayesian network)亦称“信念网"( belief network),它借助有向无环图( Directed Acyclic Graph,简称DAG)来刻画属性之间的依赖关系,并使用条件概率表( Conditional Probability Table,简称CPT来描述属性的联合概率分布。
具体来说,一个贝叶斯网
B
B
B由结构
G
G
G和参数
Θ
\Theta
Θ两部分构成,即
B
=
<
G
,
Θ
>
B=<G,\Theta>
B=<G,Θ>。网络结构
G
G
G是一个有向无环图,其每个结点对应于一个属性,若两个属性有直接依赖关系,则它们由一条边连接起来;参数
Θ
\Theta
Θ定量描述这种依赖关系,假设属性
x
i
x_i
xi,在
G
G
G中的父结点集为
π
i
\pi_i
πi,则
Θ
\Theta
Θ包含了每个属性的条件概率表
θ
x
i
∣
π
i
=
P
B
(
x
i
∣
π
i
)
\theta_{x_i|\pi_i}= P_B(x_i|\pi_i)
θxi∣πi=PB(xi∣πi)
作为一个例子,图7.2给出了西瓜问题的一种贝叶斯网结构和属性“根蒂”的条件概率表,从图中网络结构可看出,“色泽”直接依赖于“好瓜”和“甜度”,而“根蒂”则直接依赖于“甜度”;进一步从条件概率表能得到”根蒂”对“甜度”量化依赖关系,如
P
(
根蒂
=
硬挺
∣
甜度
=
高
)
=
0.1
P(根蒂=硬挺|甜度=高)=0.1
P(根蒂=硬挺∣甜度=高)=0.1等。
- 找出有向图中的所有V型结构,在V型结构的两个父结点之间加上一条无向边;
- 将所有有向边改为无向边。
同父、顺序和V型结构的发现以及有向分离的提出推动了因果发现方面的研究.
由此产生的无向图称为“道德图”( moral graph),也有译为“端正图”;令父结点相连的过程称为“道德化”(moralization).
“道德化”的含义:孩子的父母应建立牢靠的关系,否则是不道的。
基于道德图能直观、迅速地找到变量间的条件独立性.假定道德图中有变量
x
,
y
x,y
x,y和变量集合
z
⃗
{
z
i
}
\vec{z}\{z_i \}
z{zi},若变量
x
x
x和
y
y
y能在图上被
z
z
z分开,即从道德图中将变量集合
z
z
z去除后,
z
z
z和
y
y
y分属两个连通分支,则称变量
x
x
x和
y
y
y被
z
z
z有向分离,
x
⊥
y
|
z
x⊥y|z
x⊥y|z成立.例如,图7.2所对应的道德图如图7.4所示,从图中能容易地找出所有的条件独立关系:
x
3
⊥
x
4
∣
x
1
,
x
4
⊥
x
5
∣
x
2
,
x
3
⊥
x
2
∣
x
1
,
x
3
⊥
x
5
∣
x
1
,
x
3
⊥
x
5
∣
x
2
x_3⊥x_4|x_1,x_4⊥x_5|x_2,x_3⊥x_2|x_1,x_3⊥x_5|x_1,x_3⊥x_5|x_2
x3⊥x4∣x1,x4⊥x5∣x2,x3⊥x2∣x1,x3⊥x5∣x1,x3⊥x5∣x2等。
7.5.2 学习
若网络结构已知,即属性间的依赖关系已知,则贝叶斯网的学习过程相对简单,只需通过对训练样本“计数”,估计出每个结点的条件概率表即可,但在现实应用中我们往往并不知晓网络结构,于是,贝叶斯网学习的首要任务就是根据训练数据集来找出结构最“恰当”的贝叶斯网.“评分搜索”是求解这一
问题的常用办法.具体来说,我们先定义一个评分函数( score function),以此来评估贝叶斯网与训练数据的契合程度,然后基于这个评分函数来寻找结构最优的贝叶斯网.显然,评分函数引入了关于我们希望获得什么样的贝叶斯网的归纳偏好。
常用评分函数通常基于信息论准则,此类准则将学习问题看作一个数据压缩任务,学习的目标是找到一个能以最短编码长度描述训练数据的模型,此时编码的长度包括了描述模型自身所需的字节长度和使用该模型描述数据所需的字节长度.对贝叶斯网学习而言,模型就是一个贝叶斯网,同时,每个贝叶斯网描述了一个在训练数据上的概率分布,自有一套编码机制能使那些经常出现的样本有更短的编码.于是,我们应选择那个综合编码长度(包括描述网络和编码数据)最短的贝叶斯网,这就是“最小描述长度"( Minimal Description Length,简称MDL)准则。
这里我们把类别也看作一个属性,即
x
i
x_i
xi,是一个包括示例和类别的向量。
给定训练集
D
=
{
x
1
,
x
2
,
⋯
,
x
m
}
D=\{x_1,x_2,\cdots,x_m\}
D={x1,x2,⋯,xm},贝叶斯网
B
=
<
G
,
Θ
>
B=<G,\Theta>
B=<G,Θ>在D上的评分函数可写为
s
(
B
∣
D
)
=
f
(
θ
)
∣
B
∣
-
L
L
(
B
∣
D
)
−
−
−
(
7.28
)
s(B|D)=f(\theta)|B|-LL(B|D)---(7.28)
s(B∣D)=f(θ)∣B∣-LL(B∣D)−−−(7.28)
其中,|B|是贝叶斯网的参数个数;
f
(
θ
)
f(\theta)
f(θ)表示描述每个参数
θ
\theta
θ所需的字节数;而
L
L
(
B
∣
D
)
=
∑
i
=
1
m
l
o
g
P
B
(
x
i
)
−
−
−
(
7.29
)
LL(B|D)=∑_{i=1}^m logP_B(x_i)---(7.29)
LL(B∣D)=i=1∑mlogPB(xi)−−−(7.29)
是贝叶斯网
B
B
B的对数似然.显然,式(7.28)的第一项是计算编码贝叶斯网
B
B
B所需的字节数,第二项是计算
B
B
B所对应的概率分布
P
B
P_B
PB需多少字节来描述
D
D
D。于是,学习任务就转化为一个优化任务;即寻找一个贝叶斯网
B
B
B使评分函数
s
(
B
∣
D
)
s(B|D)
s(B∣D)最小。
若
f
(
θ
)
=
1
f(\theta)=1
f(θ)=1,即每个参数用1字节描述,则得到AIC( Akaike Information Criterion)j评分函数
A
I
C
(
B
∣
D
)
=
∣
B
∣
−
L
L
(
B
∣
D
)
−
−
−
(
7.30
)
AIC(B|D)=|B|-LL(B|D)---(7.30)
AIC(B∣D)=∣B∣−LL(B∣D)−−−(7.30)
若
f
(
θ
)
=
1
2
l
o
g
m
f(\theta)=\frac{1}{2} log~m
f(θ)=21log m,即每个参数用
1
2
l
o
g
m
\frac{1}{2} log~m
21log m字节描述,则得到BIC( Bayesian Information Criterion)评分函数
B
I
C
(
B
∣
D
)
=
l
o
g
m
2
∣
B
∣
−
L
L
(
B
∣
D
)
−
−
−
(
7.31
)
BIC(B|D)=\frac{log~m}{2} |B|-LL(B|D)---(7.31)
BIC(B∣D)=2log m∣B∣−LL(B∣D)−−−(7.31)
显然,若
f
(
θ
)
=
0
f(\theta)=0
f(θ)=0,即不计算对网络进行编码的长度,则评分函数退化为负对数似然,相应的,学习任务退化为极大似然估计.
不难发现,若贝叶斯网
B
=
<
G
,
Θ
>
B=<G,\Theta>
B=<G,Θ>的网络结构
G
G
G固定,则评分函数
s
(
B
∣
D
)
s(B|D)
s(B∣D)的第一项为常数(如果网络结构不同,第一项就不固定,说明对不同的网络,参与其中的参数个数不同,而不是所有的属性都参与其中).此时,最小化
s
(
B
∣
D
)
s(B|D)
s(B∣D)等价于对参数
Θ
\Theta
Θ的极大似然估计.由式(7.29)和(7-26)可知,参数
θ
x
i
∣
π
i
\theta_{x_i|\pi_i}
θxi∣πi,能直接在训练数据
D
D
D上通过经验估计获得,即
θ
x
i
∣
π
i
=
P
^
D
(
x
i
∣
π
i
)
−
−
−
(
7.32
)
\theta_{x_i|\pi_i}=\hat{P}_D(x_i|\pi_i)---(7.32)
θxi∣πi=P^D(xi∣πi)−−−(7.32)
其中
P
^
D
(
⋅
)
\hat{P}_D(\cdot)
P^D(⋅)是
D
D
D上的经验分布,即事件在训练数据上出现的频率.因此,为了最小化评分函数
s
(
B
∣
D
)
s(B|D)
s(B∣D),只需对网络结构进行搜索,而候选结构的最优参数可直接在训练集上计算得到。
不幸的是,从所有可能的网络结构空间搜索最优贝叶斯网结构是一个NP难问题,难以快速求解.有两种常用的策略能在有限时间内求得近似解:第一种是贪心法,例如从某个网络结构出发,每次调整一条边(增加、删除或调整方
向,直到评分函数值不再降低为止;第二种是通过给网络结构施加约束来削减搜索空间,例如将网络结构限定为树形结构等(半朴素叶分器可看作贝叶斯网的特例)。
7.5.3 推断
贝叶斯网训练好之后就能用来回答“查询”( query),即通过一些属性变量的观测值来推测其他属性变量的取值(类别也可着作一个属性).例如在西瓜问题中,若我们观测到西瓜色泽青绿、敲声浊响、根蒂缩,想知道它是否成熟、甜度如何,这样通过已知变量观测值来推测待查询变量的过程称为“推断”"( inference),已知变量观测值称为“证据”( evidence)。
最理想的是直接根据贝叶斯网定义的联合概率分布来精确计算后验概率,不幸的是,这样的“精确推断”已被证明是NP难的 ;换言之,当网络结点较多、连接稠密时,难以进行精确推断,此时需借助“近似推断”,通过降低精度要求,在有限时间内求得近似解,在现实应用中,贝叶斯网的近似推断常使用吉布斯采样 (Gibbs sampling)来完成,这是一种随机采样方法,我们来看看它是如何工作的。
令
Q
=
{
Q
1
,
Q
2
,
⋯
,
Q
n
}
Q=\{Q_1,Q_2,\cdots,Q_n\}
Q={Q1,Q2,⋯,Qn}表示特查询变量,
E
=
{
E
1
,
E
2
,
⋯
,
E
n
}
E=\{E_1,E_2,\cdots,E_n\}
E={E1,E2,⋯,En}为证据变量,已知其取值为
e
=
{
e
1
,
e
2
,
⋯
,
e
k
}
e=\{e_1,e_2,\cdots,e_k\}
e={e1,e2,⋯,ek}.目标是计算后验概率
P
(
Q
=
q
∣
E
=
e
)
P(Q=q|E=e)
P(Q=q∣E=e),其中
q
=
{
q
1
,
q
2
,
⋯
,
q
k
}
q=\{q_1,q_2,\cdots,q_k\}
q={q1,q2,⋯,qk}是待查询变量的一组取值.以西瓜问题为例,待查询变量为
Q
=
{
好瓜
,
甜度
}
Q=\{好瓜,甜度\}
Q={好瓜,甜度},证据变量为
E
=
(
色泽,敲声,根蒂
)
E=(色泽,敲声,根蒂)
E=(色泽,敲声,根蒂)且己知其取值为
e
=
{
青绿
,
浊响
,
蜷缩
}
e=\{青绿,浊响, 蜷缩\}
e={青绿,浊响,蜷缩},查询的目标值是
q
=
{
是
,
高
}
q=\{是,高\}
q={是,高},即这是好瓜且甜度高的概率有多大。
如图7.5所示,吉布斯采样算法先随机产生一个与证据
E
=
e
E=e
E=e一致的样本
q
0
q^0
q0作为初始点,然后每步从当前样本出发产生下一个样本。具体来说,在第
t
t
t次采样中,算法先假设
q
t
=
q
t
−
1
q^t=q^{t-1}
qt=qt−1,然后对非证据变量逐个进行采样改变其取值,采样概率根据贝叶斯网
B
B
B和其他变量的当前取值(即
Z
=
z
)
Z=z)
Z=z)计算获得.假定经过
T
T
T次采样得到的与
q
q
q一致的样本共有
n
q
n_q
nq个,则可近似估算出后验概率
P
(
Q
=
q
∣
E
=
e
)
≃
n
q
T
−
−
−
(
7.33
)
P(Q=q|E=e)\simeq \frac{n_q}{T}---(7.33)
P(Q=q∣E=e)≃Tnq−−−(7.33)
实质上,吉布斯采样是在贝叶斯网所有变量的联合状态空间与证据
E
=
e
E=e
E=e是一致的子空间中进行“随机漫步”( random walk).每一步仅依赖于前一步的状态,这是一个“马尔可夫链”( Markov chain)。在一定条件下,无论从什么初始状态开始,马尔可夫链第
t
t
t步的状态分布在
t
→
∞
t→\infty
t→∞时必收敛于一个平稳分布( stationary distribution;对于吉布斯采样来说,这个分布恰好是
P
(
Q
∣
E
=
e
)
P(Q|E=e)
P(Q∣E=e).因此,在
T
T
T很大时,吉布斯采样相当于根据
P
(
Q
∣
E
=
e
)
P(Q|E=e)
P(Q∣E=e)采样,从而保证了式(7.3)收于
P
(
Q
=
q
∣
E
=
e
)
P(Q=q|E=e)
P(Q=q∣E=e)。
—###############################################
**输入:**贝叶斯网
B
=
<
G
,
Θ
>
B=<G,\Theta>
B=<G,Θ>;
~~~~~~~~~~
采样次数
T
T
T;
~~~~~~~~~~
证据变量
E
E
E及其取值
e
e
e;
~~~~~~~~~~
待查询变量
Q
Q
Q及其取值
q
q
q;
过程:
1:
n
q
=
0
n_q=0
nq=0
2:
q
0
=
对
Q
随机赋初值
q^0=对Q随机赋初值
q0=对Q随机赋初值
3:for
t
=
1
,
2
,
⋯
,
T
t=1,2,\cdots,T
t=1,2,⋯,T do
4:
~~~~~
for
Q
i
∈
Q
Q_i \in Q
Qi∈Q do
5:
~~~~~
~~~~~
Z
=
E
∪
Q
Z=E \cup Q
Z=E∪Q \
{
Q
i
}
\{Q_i\}
{Qi}
6:
~~~~~
~~~~~
z
=
e
∪
q
t
−
1
z=e \cup q^{t-1}
z=e∪qt−1 \
{
q
i
t
−
1
}
\{q^{t-1}_i\}
{qit−1}
7:
~~~~~
~~~~~
根据
B
B
B计算分布
P
B
(
Q
i
∣
Z
=
z
)
P_B(Q_i|Z=z)
PB(Qi∣Z=z)
8:
~~~~~
~~~~~
q
i
t
=
根据
P
B
(
Q
i
∣
Z
=
z
)
采样所获
Q
i
取值
q^t_i=根据P_B(Q_i|Z=z)采样所获Q_i取值
qit=根据PB(Qi∣Z=z)采样所获Qi取值;
9:
~~~~~
~~~~~
q
t
=
将
q
t
−
1
中的
q
i
t
−
1
用
q
i
t
替换
q^t=将q^{t-1}中的q^{t-1}_i用q^{t}_i替换
qt=将qt−1中的qit−1用qit替换
10:
~~~~~
end for
11:
~~~~~
if
q
t
=
q
q^t=q
qt=q then
12:
~~~~~
~~~~~
n
q
=
n
q
+
1
n_q=n_q+1
nq=nq+1
13:
~~~~~
end if
14:end for
输出:
P
(
Q
=
q
∣
E
=
e
)
≃
n
q
T
P(Q=q|E=e) \simeq \frac{n_q}{T}
P(Q=q∣E=e)≃Tnq
—###############################################
7.6 EM算法
在前面的讨论中,我们一直假设训练样本所有属性变量的值都已被观测到。即训练样本是“完整”的,但在现实应用中往往会遇到“不完整”的训练样本,例如由于西瓜的根蒂已脱落,无法看出是“蜷缩”还是“硬挺”,则训练样本的“根蒂”属性变量值未知.在这种存在“未观测”变量的情形下,是否仍能对模型参数进行估计呢?
未观测变量的学名是“隐变量”( Intent variable)令
X
X
X表示已观测变量集,
Z
Z
Z表示隐变量集,
Θ
\Theta
Θ表示模型参数.若欲对
Θ
\Theta
Θ做极大似然估计,则应最大化对数似然(由于“似然”常基于指数簇函数来定义,因此对数似然及后续EM迭代过程中一般是使用自然对数
l
n
(
⋅
)
ln(\cdot)
ln(⋅))
L
L
(
Θ
∣
X
,
Z
)
=
l
n
P
(
X
,
Z
∣
Θ
)
−
−
−
(
7.34
)
LL(\Theta|X,Z)=lnP(X,Z|\Theta)---(7.34)
LL(Θ∣X,Z)=lnP(X,Z∣Θ)−−−(7.34)
然而由于
Z
Z
Z是隐变量,上式无法直接求解,此时我们可通过对
Z
Z
Z计算期望,来最大化已观测数据的对数“边际似然”( marginal likelihood)
L
L
(
Θ
∣
X
)
=
l
n
P
(
X
∣
Θ
)
=
l
n
∑
Z
P
(
X
,
Z
∣
Θ
)
−
−
−
(
7.35
)
LL(\Theta|X)=lnP(X|\Theta)=ln \sum_Z P(X,Z|\Theta)---(7.35)
LL(Θ∣X)=lnP(X∣Θ)=lnZ∑P(X,Z∣Θ)−−−(7.35)
EM(expectation-maximization)算法是常用的估计参数隐变量的利器,它是一种迭代式的方法,其基本想法是:若参数
Θ
\Theta
Θ已知,则可根据训练数据推断出最优隐变量
Z
Z
Z的值(E步);反之,若
Z
Z
Z的值已知,则可方便地对参数
Θ
\Theta
Θ做极大似然估计(M步)
于是,以初始值
Θ
0
\Theta^0
Θ0为起点,对式(7.35),可选代执行以下步骤直至收敛。
- 基于 Θ t \Theta^t Θt推断隐变量 Z Z Z的期望,记为 Z t Z^t Zt;
- 基于已观测变量 X X X和 Z t Z^t Zt对参数 Θ t \Theta^t Θt做极大似然估计,记为 Θ t + 1 \Theta^{t+1} Θt+1
这就是EM算法的原型.
进一步,若我们不是取
Z
Z
Z的期望,而是基于
Θ
t
\Theta^t
Θt计算隐变量
Z
Z
Z的概率分布
P
(
Z
∣
X
,
Θ
t
)
P(Z|X,\Theta^t)
P(Z∣X,Θt),则EM算法的两个步骤是
- E步( Expectation):以当前参数
Θ
t
\Theta^t
Θt(推断隐变量分布
P
(
Z
∣
X
,
Θ
t
)
P(Z|X,\Theta^t)
P(Z∣X,Θt),并计算对数似然
L
L
(
Θ
∣
X
,
Z
)
LL(\Theta|X,Z)
LL(Θ∣X,Z)关于
Z
Z
Z的期望
Q ( Θ , Θ t ) = E Z ∣ X , Θ t L L ( Θ ∣ X , Z ) − − − ( 7.36 ) Q(\Theta,\Theta^t)=E_{Z|X,\Theta^t}LL(\Theta|X,Z)---(7.36) Q(Θ,Θt)=EZ∣X,ΘtLL(Θ∣X,Z)−−−(7.36) - M步(Maximization):寻找参数最大化期望似然,即
Θ t + 1 = a r g m a x Θ Q ( Θ , Θ t ) − − − ( 7.37 ) \Theta^{t+1}=arg~max_{\Theta} Q(\Theta,\Theta^t)---(7.37) Θt+1=arg maxΘQ(Θ,Θt)−−−(7.37)
简要来说,EM算法使用两个步骤交替计算:第一步是期望(E)步,利用当前估计的参数值来计算对数似然的期望值;第二步是最大化(M)步,寻找能使
E步产生的似然期望最大化的参数值.然后,新得到的参数值重新被用于E步……直至收敛到局部最优解
事实上,隐变量估计问题也可通过梯度下降等优化算法求解,但由于求和的项数将随着隐变量的数目以指数级上升,会给梯度计算带来麻烦;而EM算法则可看作一种非梯度优化方法.