最大熵模型简介

一 原理简介

最大熵原理是一种选择随机变量统计特性最符合客观情况的准则，也称为最大信息原理。在投资时常常讲不要把所有的鸡蛋放在一个篮子里，这样可以降低风险。在信息处理中，这个原理同样适用。在数学上，这个原理称为最大熵原理。

二 熵

可以考虑从$\log$函数的性质方面考虑熵的定义，下面是熵定义的公式。
$$H(X)=-\sum _{x}p(x)\log p(x)=-\sum _{i=1}^{n}p(x_i)\log p(x_i)$$

$H(X)$就被称为随机变量$X$的熵,它是表示随机变量不确定的度量，是对所有可能发生的事件产生的信息量的期望。从公式可得，随机变量的取值个数越多，状态数也就越多，信息熵就越大，混乱程度就越大。当随机分布为均匀分布时，熵最大，且$0\le H(X)\le \log n$。

• 简单案例
  现在考虑一个具有4种可能的状态的随机变量 { a , b , c , d } \begin{aligned}\left\{ a,b,c,d\right\}\end{aligned} {a,b,c,d}​，每个状态各自的概率为$\frac{1}{2}$。这种情形下的熵为：
  $$H(X)=-4\times \frac{1}{4}{\log }_2\frac{1}{4}=2bits$$
  假如各个状态的取值为$\begin{array}{r}\left(\frac{1}{2}\frac{1}{4}\frac{1}{8}\frac{1}{8}\right)\end{array}$，这时熵为:

$$H(X)=-\frac{1}{2}\times {\log }_2\frac{1}{2}-\frac{1}{4}\times {\log }_2\frac{1}{4}-\frac{1}{8}\times {\log }_2\frac{1}{8}-\frac{1}{8}\times {\log }_2\frac{1}{8}=1.75bits$$
我们可以看到，非均匀分布比均匀分布的熵要小。

• 证明 $0\le H(X)\le \log n$
  利用拉格朗日乘子法证明：由$p(1)+p(2)+\cdots +p(n)=1$得
  目标函数:$f(p(1),p(2),\cdots ,p(n))=-(p(1)\log p(1)+p(2)\log p(2)+\cdots +p(n)\log p(n))$
  约束条件$g(p(1),p(2),\cdots ,p(n),\lambda )=p(1)+p(2)+\cdots +p(n)-1=0$
  • 1 定义拉格朗日函数
    $L(p(1),p(2),\cdots ,p(n),\lambda )=-(p(1)\log p(1)+p(2)\log p(2)+\cdots +p(n)\log p(n))+\lambda (p(1)+p(2)+\cdots +p(n)-1)$
  • 2 $L(p(1),p(2),\cdots ,p(n),\lambda )$分别对$p(1),p(2),\cdots ,p(n),\lambda$求偏导数,令偏导数为0
    $$\begin{array}{rl}\lambda -\log (e\cdot p(1))& =0\\ \lambda -\log (e\cdot p(2))& =0\\ \cdots \cdots & \\ \lambda -\log (e\cdot p(n))& =0\\ p(1)+p(2)+\cdots +p(n)-1& =0\end{array}$$
  • 3 求出$p(1),p(2),\cdots ,p(n)$的值
    解方程容易求的：$p(1)=p(2)=\cdots =p(n)=\frac{1}{n}$,代入$f(p(1),p(2),\cdots ,p(n))$中得到目标函数的极值为$f(\frac{1}{n},\frac{1}{n},\cdots ,\frac{1}{n})=-\log \frac{1}{n}=\log n$
    易知$\log 1=0$,容易证得: $0\le H(X)\le \log n$。从熵的定义方面思考这一公式，熵是对信息不确定性的度量。当一件事情完全确定,概率为1时，这个时候熵自然最小；当事情等概率发生，这时候信息的不确定性最大，熵自然最大。

条件熵

条件熵$H(Y|X)$表示在已知随机变量$X$的条件下随机变量$Y$的不确定性。
$$\begin{array}{rl}H(Y|X)& =\sum _{x}p(x)H(Y|X=x)\\ & =-\sum _{x}p(x)\sum _{y}p(y|x)\log p(y|x)\\ & =-\sum _x\sum _{y}p(x,y)\log p(y|x)\\ & =-\sum _{x,y}p(x,y)\log p(y|x)\end{array}$$

三 最大熵原理

最大熵原理是概率模型学习的一个准则，最大熵原理认为，学习概率模型时，在所有可能的概率模型（分布）中，熵最大的模型是最好的模型，通常用约束条件来确定概率模型的集合，所以，最大熵原理也可以表述为在满足约束条件的模型集合中选取最大的模型。
“最大熵”这个名词听起来很深奥，但是它的原理很简单,就是要保留全部的不确定性，将风险降到最小。让我们来看看一个实际案例。

假如现在有一个色子，让你猜测每个面向上的概率分别是多少?
最合理的猜测：各个面的概率均为1/6，从投资的角度来说,这是风险最小的做法。
假如现在告诉你，这个色子被做过手脚，已知一点向上的概率是1/3，在这种情况下，其他面向上的概率是多少?
除去一点的概率是1/3，其余的均是2/15。也就是在满足已知条件的情况下，将其余各点的概率均分。
这种基于直觉的猜测恰恰符合了最大熵的原理。

四 最大熵模型的简单实例

我们通过一个简单的例子来阐述下最大熵原理的应用。
我们以英汉翻译为例：对于英语中的“take”，其对应汉语的翻译有：
$(t_1)$“抓住”：The mother takes her child by the hand. 母亲抓住孩子的手。
$(t_2)$“拿走”：Take the book home.把书拿回家。
$(t_3)$“乘坐”：to take a bus to work.乘坐公共汽车上班。
$(t_4)$“量”：Take your temperature.量一量你的体温
$(t_5)$“装”：The suitcase wouldn’t take another thing.这个衣箱不能装别的东西了。
$(t_6)$“花费”：It takes a lot of money to buy a house.买一所房子要花一大笔钱。
$(t_7)$“理解、领会”：How do you take this package? 你怎么理解这段话？
假设对于所有的英文“take”,只有这七种翻译。则存在着如下限制：
$$p(t_1|x)+p(t_2|x)+\cdots +p(t_7|x)=1$$
$p(t_i|x)(1\le i\le 7)$表示在一个含有单词$take$的英文句子中，$take$翻译成$t_i$的概率。
在这个限制下，对每种翻译赋予均等一致的几率为:$p(t_1|x)=p(t_2|x)=\cdots =p(t_7|x)=1/7$
但是对于“take”,我们通过统计发现它的前两种翻译$(t_1)$和$(t_2)$是常见的，假设满足如下条件：$p(t_1|x)+p(t_2|x)=2/5$
那么在上面两个条件的约束下，分配给各个翻译的概率分布情况有许多，但是最一致的分布为：
$p(t_1|x)=p(t_2|x)=1/5$
$p(t_3|x)=p(t_4|x)=p(t_5|x)=p(t_6|x)=p(t_7|x)=3/25$
前面已经证明，最一致的分布具有最大的熵值。

但上面的限制，都没有考虑上下文的环境，翻译效果不好，因此我们引入特征。
例如，英文“take”翻译为“乘坐”的概率很小，但是当“take”后面跟一个交通工具的名字“bus”时，它翻译成“乘坐”的概率就变得非常大。为了表示take跟有bus时，翻译成“乘坐”的事件，我们引入二值函数：
$$f(x,y)=\{\begin{array}{ll}1& \text{1 if y=“乘坐” and ˆnext(x)=“bus”}\\ 0& \end{array}$$
x表示上下文环境，这里可以看作是含有单词take的一个英文短语，而y代表输出，对应着“take”的中文翻译。^next (x) 看作是上下文环境x的一个函数，表示x中跟在单词take后的一个单词为“bus”。这样一个函数我们称作一个特征函数，或者简称一个特征。特征是为了使我们的模型具有更强的泛化能力。
引入诸如上述公式的特征，它们对概率分布模型加以限制，求在限制条件下具有一致分布的模型，该模型熵值最大。

五 拉格朗日对偶问题

对偶问题是利用拉格朗日对偶性将原始问题转换为对偶问题，通过解对偶问题得到原始问题的解。
考虑一般化问题:
$$\begin{array}{r}mi{n}_{x\in R^n}f(x)\\ s.t.c_i(x)\le 0i=1,2,\cdots k\\ h_j(x)=0j=1,2,\cdots t\end{array}$$
引入拉格朗日函数,已知约束条件为$k+t$个,所以:
$$\begin{array}{r}L(x,\alpha ,\beta )=f(x)+\sum _{i=1}^k{a}_i{c}_i(x)+\sum _{j=1}^t{\beta }_j{h}_j(x)\\ a_i,{\beta }_j是拉格朗日乘子,a_i\ge 0\end{array}$$
引入拉格朗日函数就是为了将对$x$有限制条件的$f(x)$的最优化问题转为对$x,\alpha ,\beta$没有限制条件的$L(x,\alpha ,\beta )$极值问题。我们再定义一个函数:
$$\begin{array}{rl}{\theta }_p(x)& =ma{x}_{\alpha ,\beta }L(x,\alpha ,\beta )\\ & =ma{x}_{\alpha ,\beta }(f(x)+\sum _{i=1}^k{a}_i{c}_i(x)+\sum _{j=1}^t{\beta }_j{h}_j(x))\end{array}$$

• x不满足原始问题的约束条件,即有$c(x)>0$或$h(x)\ne 0$。那么上式求解max的解为无穷大,因为这个时候可以取$a_i=\infty ,{\beta }_j{h}_j=\infty$

• x满足原始问题的约束条件,则上式取max时,$\alpha =0,h(x)=0$,这时发现上式的最大值条件不存在了，最大值就是定值$f(x)$
  所以在满足原始问题约束的条件下,${\theta }_p(x)=f(x)$。去掉原始问题中的$s.t.$条件，得到原始问题的等价问题:
  $$mi{n}_{x}f(x)=mi{n}_x{\theta }_p(x)=mi{n}_{x}ma{x}_{\alpha ,\beta }L(x,\alpha ,\beta )$$

• 对偶问题：若原始问题和对偶问题都有最优值，则对偶问题最优值$d^{\ast }$ <=原始问题最优值$p^{\ast }$
  $$d^{\ast }=ma{x}_{\alpha ,\beta }minL(x,\alpha ,\beta )\le mi{n}_{x}ma{x}_{\alpha ,\beta }L(x,\alpha ,\beta )=p^{\ast }$$
  直观理解：清华北大的差生的成绩也比一般学校的优生成绩好。或者说:从一群人中挑选10个人，先挑20个矮个子再从中选择10个高个子所得到人的身高是小于等于先挑20个高个子再从其中选择10个矮个子的人身高。

• 求解对偶问题的好处
  对于对偶问题来说，我们求解最小化部分时没有任何限制条件，而没有限制条件的最小化问题的解一定是在求得x的偏导数=0处，那我们就能得到一些等式，将这些等式代入拉格朗日函数中就可以简化计算。

六 模型框架形式化描述

假设对于训练数据有一个样本集合 { ( x 1 , y 1 ) , ( x 2 , y 2 ) , ⋯   , ( x n , y n ) } \begin{aligned}\left\{ (x_1,y_1),(x_2,y_2),\cdots,(x_n,y_n) \right\}\end{aligned} {(x1​,y1​),(x2​,y2​),⋯,(xn​,yn​)}​，每一个$x_i(1\le i\le n)$表示一个上下文，$y_i(1\le i\le n)$表示对应的结果。对于这个训练样本，我们得到$(x,y)$的经验分布，定义如下：
$$\begin{array}{r}\tilde{p}(x,y)=\frac{1}{N}\times numberoftimesthat(x,y)occursinthesample\end{array}$$
要对上面大小为N的训练样本集合建立统计模型，可利用的是样本集合的统计数据。模型中特征函数的引入，使模型依赖于上下文的信息。假设我们给出n个特征函数$f_i$，对每个特征条件进行限制：期望概率值等于经验值
p ( f i ) = p ~ ( f i ) i ∈ { 1 , 2 , ⋯   , n } \begin{aligned} p(f_i)= \widetilde p(f_i) i \in \left\{1,2,\cdots,n \right\} \end{aligned} p(fi​)=p ​(fi​)i∈{1,2,⋯,n}​
其中期望值和经验值分别为:
$$\begin{array}{rl}p(f)& =\sum _{x,y}\tilde{p}(x)p(y|x)f(x,y)\\ \tilde{p}(f)& =\tilde{p}(x,y)f(x,y)\end{array}$$
要求出最优的$P(y|x)$值，我们要得到一个最为一致分布的模型，条件熵作为衡量一致的标准。熵的最小值是0，这时模型没有任何不确定性；熵的最大值是$log|Y|$, 即在所有可能的$y$上的均匀分布。
求在限制条件下具有最大熵值的模型，$C$表示所有可能满足限制条件的概率分布模型的集合
m a x p ∈ C H ( P ) = − ∑ x , y p ~ ( x ) p ( y ∣ x ) log ⁡ p ( y ∣ x ) C = { p ∈ P ∣ p ( f i ) = p ~ ( f i ) f o r   i   i n   ( 1 , 2 , ⋯   , n ) } \begin{aligned} max_{p \in C} H(P) = - \sum_{x,y} \widetilde{p}(x)p(y|x)\log p(y|x)\\ C = \left \{p \in P|p(f_i)=\widetilde{p}(f_i) for\ i\ in\ (1,2,\cdots,n)\right\} \end{aligned} maxp∈C​H(P)=−x,y∑​p ​(x)p(y∣x)logp(y∣x)C={p∈P∣p(fi​)=p ​(fi​)for i in (1,2,⋯,n)}​
$H(P)$满足以下限制:
p ( y ∣ x ) ≥ 0   f o r   a l l   x ,   y ∑ y p ( y ∣ x ) = 1   f o r   a l l   x ∑ x , y p ~ ( x ) p ( y ∣ x ) f ( x , y ) = ∑ x , y p ~ ( y ∣ x ) f ( x , y )   f o r   i ∈ { 1 , 2 , ⋯   , n } \begin{aligned} p(y|x) &\geq 0\ for \ all \ x,\ y\\ \sum_y p(y|x) &= 1\ for\ all\ x \\ \sum_{x,y} \widetilde{p}(x)p(y|x)f(x,y) &= \sum_{x,y} \widetilde{p}(y|x)f(x,y) \ for\ i \in \left\{ 1,2,\cdots,n \right\} \end{aligned} p(y∣x)y∑​p(y∣x)x,y∑​p ​(x)p(y∣x)f(x,y)​≥0 for all x, y=1 for all x=x,y∑​p ​(y∣x)f(x,y) for i∈{1,2,⋯,n}​

按照最优化的习惯，将最大值问题改写为等价的求最小值问题。
$$\begin{array}{r}mi{n}_{p\in C}-H(P)=\sum _{x,y}\tilde{p}(x)p(y|x)\log p(y|x)\end{array}$$
将约束最优化的原始问题转化为无约束最优化的对偶问题，通过求解对偶问题求解原始问题。具体作法为：为每一个$f_i$引入一个参数$w_i$。因此定义拉格朗日函数$L(p,w)$
$$\begin{array}{rl}L(p,w)& =-H(p)+w_0(1-\sum _{y}p(y|x))+\sum _{i=1}^n{w}_i(p(f_i)-{\tilde{p}}_{f_i})\\ & =\sum _{x,y}\tilde{p}(x)p(y|x)\log p(y|x)+w_0(1-\sum _{y}p(y|x))+\sum _{i=1}^n{w}_i(\sum _{x,y}\tilde{p}(x_i)p(y|x)f_i(x,y)-\tilde{p}(x,y)f_i(x,y))\end{array}$$
最优化的原始问题是
$$\begin{array}{r}mi{n}_{p\in c}ma{x}_{w}L(p,w)\end{array}$$
对偶问题是：
$$\begin{array}{r}ma{x}_{w}mi{n}_{p\in c}L(p,w)\end{array}$$
首先，求解对偶问题内部的极小化问题$mi{n}_{p\in c}L(p,w)$，它是$w$的函数，将其记作:
$$\begin{array}{r}{\psi }_w=mi{n}_{p\in c}L(p,w)=L(p_w,w)\end{array}$$
${\psi }_w$称为对偶函数。同时，将其解记作
$$\begin{array}{r}{P}_w=argmi{n}_{p\in c}L(p,w)=P_w(y|x)\end{array}$$
具体地，求L(P,W)对P(y|x)的偏导数，再令其偏导数为0时，解得：
$$\begin{array}{r}P(y|x)=\frac{exp({\sum }_{i=1}^n{w}_i{f}_i(x,y))}{exp(1-w_0)}\end{array}$$
又因为 ${\sum }_{y}P(y|x)=1$
得到:
$$\begin{array}{r}{P}_w(y|x)=\frac{exp({\sum }_{i=1}^n{w}_i{f}_i(x,y))}{z_w}\\ z_w(x)=\sum _{y}exp(\sum _{i=1}^n{w}_i{f}_i(x,y))\end{array}$$
之后，求解对偶问题外部的极大化问题
$$ma{x}_w\psi (w)$$
将其记为$w^{\ast }$，即
$$\begin{array}{r}{w}^{\ast }=argma{x}_w\psi (w)\end{array}\text{\hspace{1em}(1)}$$
这就是说，可以应用最优化算法求对偶函数$\psi (w)$的极大化，得到$w^{\ast }$，用来表示$P^{\ast }=C$，这里，$P^{\ast }=P_w^{\ast }=P_w^{\ast }(y|x)$是学习到的最优模型（最大熵模型）。公式（1）没有显示的解析解,由于目标函数是一个 凸函数,所以可以借助多种优化方法来进行求解,并且能保证得到全局最优解。

七 GIS算法原理

最大熵模型参数训练的任务就是选取有效的特征$f_i$及其权重$w_i$。由于特征数量大，模型复杂，因此需要依赖数值解的方法，下面介绍GIS算法。GIS算法要求对训练样本集中每个实例的任意$(x,y)\in X\times Y$，特征函数之和为常数，即对每个实例的k个特征函数均满足${\sum }_{i=1}^k{f}_i(x,y)=C$（C为一常数）。如果该条件不能满足，则在训练集中取:$C=ma{x}_{x\in X,y\in Y}{\sum }_{i=1}^k{f}_i(x,y)$，并增加一个特征$f_l:f_l(x,y)=C-{\sum }_{i=1}^k{f}_i(x,y)$。其中$l=k+1$。与其他特征函数不一样,$f_l(x,y)$的取值范围为：0~C。
GIS的算法流程如下：(N表示样本总数，n表示特征总数)
（1）初始化:$w[1\cdots l]=0$；
（2）根据公式$\tilde{p}(f)=\tilde{p}(x,y)f(x,y)$计算每个特征函数$f_i$的训练样本期望值$\tilde{p}(f)$；
（3）执行如下循环，迭代计算特征函数的模型期望值：

• [1] 利用如下公式计算$\tilde{p}(x,y)$
  $$\begin{array}{r}\tilde{P}(y|x)=\frac{exp({\sum }_{i=1}^l{w}_i{f}_i(x,y))}{z(x)}\\ z(x)=\sum _{y}exp(\sum _{i=1}^l{w}_i{f}_i(x,y))\end{array}$$
• [2]若满足终止条件，则结束迭代；否则修正$w$
  $$\begin{array}{r}{w}^{n+1}=w^n+\frac{1}{C}\ln (\frac{\tilde{P}(f_i)}{P_{x(n)}(f_i)})n为循环迭代次数\end{array}$$

（4）算法结束，确定$w$，算出每个$\tilde{p}(y|x)$
迭代终止的条件可以为限定的迭代次数，也可以是对数似然$L(p)$的变化值小于某个阈值$ϵ$。
$$\begin{array}{r}|L_{n+1}-L_n|<ϵ\\ L(p)=\sum _{x,y}\tilde{p}(x,y)\log p(y|x)\end{array}$$

八 python代码实现

1. 定义MaxEnt类，主要属性如下：

import math
from collections import defaultdict

class MaxEnt:
    def __init__(self):
        self._samples = []  # 样本集, 元素是[y,x1,x2,...,xn]的元组
        self._Y = set([])  # 标签集合,相当于去重之后的y
        self._numXY = defaultdict(int)  # Key是(xi,yi)对，Value是count(xi,yi)
        self._N = 0  # 样本数量
        self._n = 0  # 特征对(xi,yi)总数量
        self._xyID = {}  # 对(x,y)对做的顺序编号(ID), Key是(xi,yi)对,Value是ID
        self._C = 0  # 样本最大的特征数量,用于求参数时的迭代，见IIS原理说明
        self._ep_ = []  # 样本分布的特征期望值
        self._ep = []  # 模型分布的特征期望值
        self._w = []  # 对应n个特征的权值
        self._lastw = []  # 上一轮迭代的权值
        self._EPS = 0.01  # 判断是否收敛的阈值

2. 加载训练数据
   训练数据的格式为“标注+特征”的格式，特征之间以空格分隔。
   

# 加载训练数据
def load_data(self, filename):
    for line in open(filename, "r", encoding='utf-8'):
        sample = line.strip().split()
        if len(sample) < 2:  # 至少：标签+一个特征
            continue
        y = sample[0]  # 标注
        X = sample[1:] # 特征
        self._samples.append(sample)  # label + features
        self._Y.add(y)  # label
        for x in set(X):  # set给X去重
            self._numXY[(x, y)] += 1

3. 初始化参数
   样本数量、特征对数、特征权重、上一轮迭代的权重初始化为0

# 初始参数
def _initparams(self):
    self._N = len(self._samples)  # 样本数量
    self._n = len(self._numXY)  # (xi,yi)总数量，没有做任何特征提取操作，直接操作特征
    self._C = max([len(sample) - 1
                    for sample in self._samples])  # 样本最大的特征数量
    self._w = [0.0] * self._n  # 初始权重设为零
    self._lastw = self._w[:]
    self._sample_ep()  # 计算样本期望

4. 判断是否收敛

# 判断是否收敛
def _convergence(self):
    for w, lw in zip(self._w, self._lastw):
        if math.fabs(w - lw) >= self._EPS:
            return False
    return True

5. 计算z(x)
   输入的$X$为每个标签对应的特征

def _zx(self, X):
    ZX = 0.0
    for y in self._Y:
        sum = 0.0
        for x in X:
            if (x, y) in self._numXY:  # 遍历键(xi,yi)
                sum += self._w[self._xyID[(x, y)]]
        ZX += math.exp(sum)
    return ZX

6. 初始化样本分布中的特征期望值
   每个特征的期望初始化为：$\frac{当前特征数}{总特征数}$为每个特征标记一个ID

def _sample_ep(self):
    self._ep_ = [0.0] * self._n

    for i, xy in enumerate(self._numXY):
        self._ep_[i] = self._numXY[xy] * 1.0 / self._N
        self._xyID[xy] = i

7. 计算p(y|x)

def _pyx(self, X):
    ZX = self._zx(X)
    results = []
    for y in self._Y:
        sum = 0.0
        for x in X:
            if (x, y) in self._numXY:  # 这个判断相当于指示函数的作用
                sum += self._w[self._xyID[(x, y)]]
        pyx = 1.0 / ZX * math.exp(sum)
        results.append((y, pyx))
    return results

8. 计算模型分布的特征期望

def _model_ep(self):
    self._ep = [0.0] * self._n

    for sample in self._samples:
        X = sample[1:]
        pyx = self._pyx(X)
        for y, p in pyx:
            for x in X:
                if (x, y) in self._numXY:
                    self._ep[self._xyID[(x, y)]] += p * 1.0 / self._N

9. 训练数据
   设置迭代次数maxiter=100，每轮迭代，保存上一轮权重。

def train(self, maxiter=1000):
    self._initparams()
    for i in range(0, maxiter):
        print("Iter:%d..." % i)
        self._lastw = self._w[:]  # 保存上一轮权值
        self._model_ep()
        # 更新每个特征的权值
        for i, w in enumerate(self._w):
            # 参考GIS迭代公式
            self._w[i] += 1.0 / self._C * math.log(
                self._ep_[i] / self._ep[i])
        print(self._w, "数量:", len(self._w))
        # 检查是否收敛
        if self._convergence():
            break

10. 进行预测

def predict(self, input):
    X = input.strip().split("\t")
    prob = self._pyx(X)
    return prob

11. 主函数调用

if __name__ == "__main__":
    maxent = MaxEnt()
    maxent.load_data(r"./data.txt")
    maxent.train()
    # 特征 play	outlook	temperature	humidity windy
    print(maxent.predict("sunny\thot\thigh\tFALSE"))
    print(maxent.predict("overcast\thot\thigh\tFALSE"))
    print(maxent.predict("sunny\tcool\thigh\tTRUE"))
    print(maxent.predict(""))

12. 实验结果
    

参考资料

统计学习方法
最大熵模型简介[例子+推导+GIS求解]