1.最大熵原理
最大熵原理认为,学习概率模型时,在所有可能的概率分布中,熵最大模型为最好的模型。因此在满足约束条件的情况下,最大熵原理选择的使熵最大的模型。
熵的公式为:
H
(
P
)
=
−
∑
x
P
(
x
)
l
o
g
P
(
x
)
H(P)=-\sum_x P(x)logP(x)
H(P)=−∑xP(x)logP(x)
2.最大熵模型
将上述最大熵原理运用到模型的选择上即为最大熵模型。
考虑现在有样本集
(
x
1
,
y
1
)
,
(
x
2
,
y
2
)
.
.
.
(
x
n
,
y
n
)
(x_1,y_1),(x_2,y_2)...(x_n,y_n)
(x1,y1),(x2,y2)...(xn,yn),如何来建立模型预测新的P(y|x)。
首先,原有的样本集即为我们的约束条件。
即对于上述的每一个样本,以
(
x
1
,
y
1
)
(x_1,y_1)
(x1,y1)为例,特征函数f(x,y)=1当前仅当
x
=
x
1
,
y
=
y
1
x=x_1,y=y_1
x=x1,y=y1,这为一个约束条件。
对于此f(x,y),我们可以求出其的关于经验分布
P
~
(
x
,
y
)
\widetilde{P}(x,y)
P
(x,y)的函数期望,即为:
E
p
~
(
f
)
=
∑
x
,
y
P
~
(
x
,
y
)
f
(
x
,
y
)
E_{\widetilde{p}}(f)=\sum_{x,y}\widetilde{P}(x,y)f(x,y)
Ep
(f)=x,y∑P
(x,y)f(x,y)
即求在所有的样本中,满足
x
=
x
1
,
y
=
y
1
x=x_1,y=y_1
x=x1,y=y1的样本所占的比例(因
x
!
=
x
1
∣
∣
y
!
=
y
1
x!=x1||y!=y1
x!=x1∣∣y!=y1时,f(x,y)=0)。同样地,我们要求的
P
(
y
∣
x
)
P(y|x)
P(y∣x)也应该上述约束条件,我们可以利用经验分布
P
~
(
x
)
\widetilde{P}(x)
P
(x)与P(y|x)的联合分布来对f(x,y)的期望值进行求解,
E
p
(
f
)
=
∑
x
,
y
P
~
(
x
)
P
(
y
∣
x
)
f
(
x
,
y
)
E_{{p}}(f)=\sum_{x,y}\widetilde{P}(x)P(y|x)f(x,y)
Ep(f)=∑x,yP
(x)P(y∣x)f(x,y).然后令
E
p
~
(
f
)
=
E
p
(
f
)
E_{\widetilde{p}}(f)=E_{{p}}(f)
Ep
(f)=Ep(f)。
即关于特征函数f(x,y)的P(y|x)应满足的约束条件如下所示:
∑
x
,
y
P
~
(
x
,
y
)
f
(
x
,
y
)
=
∑
x
,
y
P
~
(
x
)
P
(
y
∣
x
)
f
(
x
,
y
)
\sum_{x,y}\widetilde{P}(x,y)f(x,y) = \sum_{x,y}\widetilde{P}(x)P(y|x)f(x,y)
x,y∑P
(x,y)f(x,y)=x,y∑P
(x)P(y∣x)f(x,y)
而现在有n个样本,因此类似的有n个约束条件,因此会有n个上述类型的特征函数。
在要求的P(y|x)满足上述N个约束条件的前提下,然后利用最大熵原理,选择令
P(y|x)熵最大的模型,即要求的P为:
其中条件熵的计算公式可为
H
(
P
)
=
∑
x
P
~
(
x
)
∑
y
−
P
(
y
∣
x
)
l
o
g
P
(
y
∣
x
)
H(P)=\sum_{x}\widetilde{P}(x)\sum_{y}-P(y|x)logP(y|x)
H(P)=x∑P
(x)y∑−P(y∣x)logP(y∣x)
。
3.最大熵模型的学习
在上文介绍了最大熵原理以及对最大熵模型的定义之后,接下来涉及如何对上述最大熵模型进行求解。
在第2节对最大熵模型的定义中,我们可以将其看做一个带约束的最优化问题。因此可以利用经典的拉格朗日乘子法将其转换为无约束优化问题:
首先,我们列出拉格朗日函数:
L
(
P
,
w
)
=
−
H
(
P
)
+
w
0
(
1
−
∑
y
P
(
y
∣
x
)
)
+
∑
i
=
1
n
w
i
(
E
P
~
(
f
i
)
−
E
P
(
f
i
)
)
L(P,w)=-H(P)+w_0(1-\sum_y{P(y|x)})+\sum_{i=1}^{n}{w_i(E_{\widetilde{P}}(f_i)-E_P(f_i))}
L(P,w)=−H(P)+w0(1−y∑P(y∣x))+i=1∑nwi(EP
(fi)−EP(fi))
然后利用拉格朗日的对偶性(具体可参考链接:link link link link进行拉格朗日相关知识的学习)
此处和下文中公式推导求解不再赘述,主要是参考李航书中的内容,具体也可参考这两篇博客link,link。
我们可通过
∂
L
(
P
,
w
)
∂
P
=
0
\frac{\partial L(P,w)}{\partial P}=0
∂P∂L(P,w)=0来求出P(y|x)。
最终求出来的结果如上图,即我们先求出来了对偶问题
m
i
n
P
∈
C
L
(
P
,
w
)
\underset{P\in C}{min}L(P,w)
P∈CminL(P,w)的结果,然后将求出来的
P
w
(
y
∣
x
)
P_w(y|x)
Pw(y∣x)代入到
Ψ
(
w
)
=
L
(
P
w
,
w
)
\Psi(w)=L(P_w,w)
Ψ(w)=L(Pw,w)中,现在转换为了求w的问题,即求无约束优化问题——对对偶问题的极大化
m
a
x
w
Ψ
(
w
)
\underset{w}{max}\Psi(w)
wmaxΨ(w)。
之后的参数w的求解可利用常用的无约束优化问题求解方法——牛顿法或者《统计学习方法》中介绍的改进的迭代尺度法,此处不再赘述。
4.最大熵模型与逻辑回归的关系
同时,在统计学习方法的6.2.4节,证明了对对对偶问题的极大化
m
a
x
w
Ψ
(
w
)
\underset{w}{max}\Psi(w)
wmaxΨ(w)等价于对条件概率分布P(y|x)的极大对数似然估计。
而其实逻辑回归中P(w|x)参数的求解也是对条件概率分布P(w|x)的极大对数似然估计,二者唯一的区别是P(w|x)的形式不同,但我们可以通过对最大熵模型加入一些条件推出逻辑回归模型。见:link.
5.总结
本次首先介绍最大熵原理,然后由最大熵原理+n个样本带来的约束条件推出最大熵模型,即为一个有约束的最优化问题,然后利用拉格朗日乘子法求解上述的有约束的最优化问题,并利用其对偶性将其转换为了先求出P(y|x)在求w的无约束的最优化问题,同时我们也证明了求出P(y|x)之后求w的最优化过程等价于直接对条件概率分布P(y|x)的极大对数似然估计,然后引出了与逻辑回归的关系,且在对最大熵模型进行一些初始条件限制之后,即与逻辑回归等价。
本次并未详细介绍公式的推导,主要介绍了最大熵模型的大体流程,具体推导可参考《统计学习方法》和上述引用的一些链接,最后也可能会贴上一个利用改进的迭代尺度法实现的最大熵模型代码。
6 代码实现
import numpy as np
import time
from collections import defaultdict
'''
minist 数据集
50000 训练集(实际使用1000)
10000 测试集
训练结果:
the acc is 0.962400 :
time span 2157.652851 s:
'''
def load_data(fileName):#加载数据
'''
:param fileName:minist训练集/测试集文件
:return: 对应的特征值X和标签值Y
'''
fr = open(fileName, 'r')
dataX = [];dataY = []
for line in fr.readlines():
curline = line.strip().split(',')
if(int(curline[0])==0):# 二分类问题 转换为区分0和非0数字
dataY.append(0) # 标签
else:
dataY.append(1)
dataX.append([int(int(num)>128) for num in curline[1:]]) # 特征
# print(dataX)
return dataX, dataY
class maxEnt:
'''
最大熵类
'''
def __init__(self,trainX,trainY,testX,testY):
self.trainX =trainX
self.trainY =trainY
self.testX =testX
self.testY =testY
self.featureNum = len(trainX[0])
self.N = len(trainX)#总的训练集数目
self.n =0 #不同特征下的(xi,y)的总数目 784*2*2 二分类
self.M =10000
self.fixy = self.calc_fixy()# #计算不同特征下出现的(x,y)对的次数
self.w = [0]*self.n #每个(xi,y)对的出现 对应一个f 总共有n个(xi,y)对 对应n个f 也对应n个w
self.xy2idDict,self.id2xyDict = self.createSearchDict() #创建xy与id之间的对应关系
self.Ep_xy = self.calcEp_xy()# 特征函数f(x,y)关于经验分布P-(x,y)的期望
def maxEntropyTrain(self,iter=500):
'''
进行最大熵模型的训练
:param iter: 迭代次数
:return:
'''
for i in range(iter):
iterStart = time.time()
Epxy = self.calcEpxy()#计算P83页的期望
sig =[0]*self.n
for j in range(self.n):
sig[j] = (1/self.M)*np.log(self.Ep_xy[j]/Epxy[j])
self.w =[self.w[q] + sig[q] for q in range(self.n)]# 没有使用array 不能直接相加减
iterEnd = time.time()
print('iter:%d:%d,time:%d'%(i,iter,iterEnd-iterStart))
return self.w
def calcEpxy(self):
'''
计算特征函数f(x,y)关于模型P(y|x)与经验分布P(x)的期望值
:return:
'''
Epxy = [0.0]*self.n
for i in range(self.N):
Pwxy = [0]*2
##计算P(y = 0 } X)
#注:程序中X表示是一个样本的全部特征,x表示单个特征,这里是全部特征的一个样本
Pwxy[0] = self.calcPwy_x(self.trainX[i],0)#计算P(y|x)的概率
##计算P(y = 1 } X)
Pwxy[1] = self.calcPwy_x(self.trainX[i],1)
for j in range(self.featureNum):
for y in range(2):#对每个xi点判断所有情况
if (self.trainX[i][j],y) in self.fixy[j]:#如果真的存在这种情况
id = self.xy2idDict[j][(self.trainX[i][j],y)]
Epxy[id] += (1/self.N)*Pwxy[y]#然后计算对应情况下的关于模型P(y|x)与经验分布P(x)的期望值
return Epxy
def calcPwy_x(self,X,y):
'''
#计算条件概率 依据最大熵模型公式计算
:param X:训练样本X
:param y: 不同的标签
:return: P(y|X)
'''
Z =0
numerator = 0
for i in range(self.featureNum):
if (X[i],y) in self.xy2idDict[i]:#xi y 出现过
index = self.xy2idDict[i][(X[i],y)]
numerator += self.w[index]
if (X[i],1-y) in self.xy2idDict[i]:#xi 1-y 出现过
index = self.xy2idDict[i][(X[i], 1-y)]
Z += self.w[index]
numerator = np.exp(numerator)
Z = numerator+np.exp(Z)# 所有情况的
return numerator/Z
def calcEp_xy(self):# 特征函数f(x,y)关于经验分布P-(x,y)的期望
Ep_xy = [0]*self.n
for i in range(self.featureNum):
for (x,y) in self.fixy[i]:
id = self.xy2idDict[i][(x,y)]
Ep_xy[id]=self.fixy[i][(x,y)]/self.N#利用(xi,y)出现次数除以总的训练个数即为期望 fixy 已是对所有x,y的遍历
return Ep_xy
def createSearchDict(self):
# 对 self.fixy中的(x,y)建立一个标号 实现feature、(x,y)对到n的映射
xy2idDict = [{} for i in range(self.featureNum)]
id2xyDict = {}
index = 0
for i in range(self.featureNum):
for (x,y) in self.fixy[i]:
xy2idDict[i][(x,y)]=index
id2xyDict[index] = (x,y)
index+=1
return xy2idDict,id2xyDict
def calc_fixy(self):
'''
#计算不同特征下出现的(x,y)对的次数
:return:字典列表
'''
fxyDict = [defaultdict(int) for i in range(self.featureNum)]#建立featureNum个字典
for i in range(len(self.trainX)):
for j in range(self.featureNum):
fxyDict[j][(self.trainX[i][j],self.trainY[i])]+=1#计算j特征下的(xi,y)出现的频率
for i in fxyDict:#对整个大字典进行去重、计算每个特征下共有多少个特征对
self.n+=len(i)
return fxyDict
def predict(self,X):
# Pwy_x = self.calcPwy_x(X,y)
# if(Pwy_x>=0.5):
# return 1
# else:
# return 0
result =[0]*2
for i in range(2):
result[i] = self.calcPwy_x(X,i)
print(result[0],result[1],self.w[0],self.w[1])
return 0 if result[0]>result[1] else 1
def test(self):
errorCnt = 0
for i in range(len(self.testX)):
if self.predict(self.testX[i])!=self.testY[i]:
errorCnt+=1
return 1-errorCnt/len(self.testX)
if __name__=="__main__":
start =time.time()
print("start read trainData")
trainX, trainY = load_data("./mnist_train/mnist_train.csv")
print("start read TestData")
testX, testY = load_data("./mnist_test/mnist_test.csv")
maxEnt = maxEnt(trainX[:1000],trainY[:1000],testX,testY)
print("start to train")
maxEnt.maxEntropyTrain(500)
print("start to test")
acc = maxEnt.test()
print("the acc is %f :"%acc)
print("time span %f s:"%(time.time()-start))
参考链接:
https://www.pkudodo.com/2018/12/05/1-7/
https://www.jianshu.com/p/e7c13002440d