统计学习方法|最大熵原理剖析及实现

最大熵的直观理解

为了引出最大熵，我们可能需要举一个所有博客都会举的例子：如果我手里拿着一个骰子，想问你扔下去后是每一面朝上的概率是多少？所有人都能抢答出来是1/6。那么为什么是1/6呢？为什么不是1/2或者1/3？

因为六个面的概率相同呀

emmm….我暂时承认你说的有点道理。可是如果这是老千手里的骰子呢？你还会答1/6吗？

可是你没说是老千手里的骰子呀

你让我哑口无言了，但为什么大家在猜之前不会去假设是不是老千的骰子这种情况呢？

因为你没说是老千手里的骰子呀

完蛋了，又绕回来了。不过我想想也是，如果要考虑老千，那么也可能要考虑骰子是不是破损了，桌面是不是有问题……完蛋了，没头了。

所以1/6最保险呀

如果我告诉你1朝上的概率是1/2呢？

那剩下的就是1/10呀

emmm…我承认在目前的对话中你是明智的，而我像个低龄的儿童。但是我要告诉你，上面几句对话，其实就是最大熵。
当我们在猜测概率时，不确定的部分我们都认为是等可能的，就好像骰子一样，我们知道有六个面，那么在概率的猜测上会倾向于1/6，也就是等可能，换句话讲，也就是倾向于均匀分布。为什么倾向均匀分布？你都告诉我了，因为这样最保险呀！当我们被告知1朝上的概率是1/2时，剩下的我们仍然会以最保险的形式去猜测是1/10。是啊，最大熵就是这样一个朴素的道理：

凡是我们知道的，就把它考虑进去，凡是不知道的，通通均匀分布！

从另一个角度看，均匀分布其实保证的是谁也不偏向谁。我说万一是老千的骰子呢？你说万一骰子有破损呢？骰子好端端地放在那咱俩地里咕噜瞎猜啥呢。咱们与其瞎猜是啥啥啥情况，干脆就均匀分布，谁也不偏向谁。
多朴素的道理啊！

朴素贝叶斯分类器的数学角度（配合《统计学习方法》食用更佳）
在最大熵章节中我们在公式上仍然参考《统计学习方法》，除此之外会补充一些必要知识点，建议阅读完毕后浏览一下代码知道详细的用法。

最大熵模型是一个分类模型，这样，我们照旧不管中间过程，统统都不管，我们最终可能会想要这样一个式子： P(Y|X)。是啊，P(Y|X)就像一个终极的大门一样，扔进入一个样本x，输出标签y为各种值的概率。我们现在手里有一堆训练数据，此外手握终极大式子P(Y|X)，是不是想到了朴素贝叶斯？我们对训练数据进行分析，获得先验概率，从而得到P(Y|X)模型进行预测。最大熵也有数据和相同的终极大式子，那么

最大熵和朴素贝叶斯有什么联系？

可以说在本质上是非常密切的，但具体的联系和区别，我希望读者正式了解了最大熵以后，我们一起放在文章末尾讨论。

文章的最初，我们需要给出最大熵公式，如果对下面这个式子不太熟悉的，可以查看《统计学习方法|决策树原理剖析及实现》一节 ：

H§表现为事件P的熵变，也称为事件P的不确定性。对于最大熵模型而言，我们一直记得那个终极式子P(Y|X)，如果放到骰子例子中，P(Y=1|X)表示扔了一个骰子，1朝上的概率，P(Y=2|X)表示2朝上的概率。我们之前怎么说的？我们只知道骰子是六面，至于其他情况我们一概不知，在不知道的情况下，我们就不能瞎做推测，换言之，我们就是要让这个P(Y|X)的不确定性程度最高。

不确定性程度最高？

可能有些读者这么一绕没回过神来，如果我告诉你点数1朝上的概率是1/2，你是不是对这个概率分布开始有了一点认识，那么会将剩下的面认为是1/10（注：1朝上概率是1/2，那么1背面的概率一定会降低，但我们没有考虑这件事情，事实上我们不应该添加任何的主观假设，因为实际上有些骰子1的背面是3，也有些是4，任何主观假设都有可能造成偏见产生错误的推断）。你推断剩余面是1/10的概率是因为我告诉了你信息，你对这个骰子开始有了一些确定的认识，也就是信息的不确定性程度——熵减少了。但事实上我没有告诉你这一信息，你对P(Y|X)的分布是很不确定的，那么我们能不能用使用下面这个式子，让H§最大化来表示我们对P(Y|X)的极大不确定呢？事实上当我们令下式中H§最大时，P(Y|X)正好等于1/6。也就是说，如果我们要让熵最大，P(Y|X)必须为1/6。换句话说，只有我们把骰子的每一面概率当做1/6时，才能充分表示我们对骰子每一面的概率分布的极度不确定性，骰子的概率这件事情的混乱度最高，同时也传递了一个信息：我们对这件事情没有添加任何主观的假设，对任何未知的东西没有偏见和主观假象，这是最保险、最安全的做法。
上式可以理解成x=托马斯全旋法扔骰子，y为各类结果的概率。那么我们并不是说x只有托马斯全旋法这一种，很可能还有王思聪吃热狗法扔骰子、边骂作者啥也不懂边扔骰子法等等等等….我们要使得在所有情况下熵都是最大，所以引出下式：

加了一个P_(x)是什么意思？就是说我在求每种x对应的y时，要考虑到x出现的概率进行相乘。换句话说，也可以认为是当前P(y|x)的权值。因为我们要保证在所有情况下总的熵是最大的，所以每个部分单独的熵乘以权值后再相加，就是整个事件的熵。上式这个模型，其实就是最大熵模型，我们要让P(y|x)的熵最大，等价于让H§最大，我们好像获得了一个新的终极大式子。目前来看，只要能够让H§为最大值，那就能获得模型啦。

emmm….不过里面的各个子式好像也不是那么好求的。所以我们暂时先忘记上面的式子，好好回想一下咱们手里有哪些信息，好好整理一下。

我认为我手里应该会有一个训练集，包含所有的样本以及样本对应的标签（训练集都没有的话啥模型也训练不了啊，那可太惨了）。

那么我们首先可以得到训练集中x，y的联合概率分布以及x的边缘分布，数据都在我们手上，这两个很简单就能求出来对吧。

在上式中v表示频数，第一个式子的分子也就表示(x, y)组合出现的次数，处于总数N即为(x, y)组合出现的概率，也就是P(X, Y)。第二个式子同理，求得P(X=x)的概率。那么现在P(X, Y)和(X)我们知道了，有什么用呢？emmm….暂时好像还没什么头绪，我们看下式，书中说这是特征函数f(x, y)关于经验分布p(x, y)的期望值，f(x, y)在x和y满足某一事实时结果为1，反之输出0。比如说训练集中某一样本记录是托马斯回旋法扔骰子，结果是5朝上。那么f(x=托马斯回旋法扔骰子|Y=5)=1。其实更简单的理解，可以认为如果(x, y)出现在训练集中，那么它就是1，因为一旦出现在训练集中，说明(x, y)对已经符合了某一样本的事实了。P上面的短横表示这个模型是基于训练集得出来的，只能被称为经验分布。

于此同时我们可以直接将P_(x, y)拆开来，拆成下式：

那么我们将P_(x, y)转换成了P_(x)*P(y|x)，发现什么了没有？里面有P_(y|x)，这个和我们要得到的最终大式子P(Y|X)很像，那么如果我们只是将其替换进去变成下面这样会有什么结果呢？

我们发现，如果我们的预测模型P(y|x)是正确的话，那P(y|x)应该等于p_(y, x)，也就是说

这两个应该相等！或者说：

那这个有啥用？emmm….我们将我们目前知道的信息整理出来了上式，老实说，目前来看没啥用……

咦？等等，我们换种思维想一下，我们知道如果做出来的模型足够正确的画，那应该上面的两个式子相等，那….这算不算一种约束？约束我一定要得到正确的模型？还记得我们要最大化的那个熵的式子吗？这个代表了我们要求得的模型。

那加上我们刚刚整理出来的约束：
是不是看起来有点像那么回事了，第三个式子也是一个约束，概率和一定为1嘛。在机器学习中不止一次见到了上图的式子形式，但通常是求min而非上图的max，那么我们再转换一下，求最大取个负号就是求最小了嘛。
好了，现在是不是知道该怎么做了？转换成拉格朗日乘子法直接求最小值H§呀。
但是我们并不能直接求到H§的最优值，在使用拉格朗日乘子法后引入了变量w，对于w，我们需要先求得L(P,W)的最大值，此时将P看成了一个常数。之后再对P操作求得-H§的最小值。也就是说：

一定要注意min和max下方的参数是不一样的，因为L是考虑到约束函数从生成的目标桉树，在求解最优的参数w以此使得我们的模型能够符合约束时P只是一个定值，此时需要求得最优的w以此保证约束条件对模型的约束性。当求得w后P满足了约束性，此时再求最优的P以此使得-H§最小化，也就是H§最大化以此得到最优模型。

与此同时，我们可以将上式转换为下式，因为L是P的凸函数，因此两式等价。

所以我们先求内部的极小化问题，即min L(P， w)，因为L是此时是关于P的凸函数，w是一个定值，所以可以直接对P求导求得最小值。H§中的P是什么？P就是P(y|x)啊（前文在H§的式子中可以显然看到，事实上我们全文也一直在为P(y|x)进行准备）。因为L是P的凸函数，那么直接求导并且令求导结果为0就好啦，从而得到下式。

我们再整理一下：

上式没有什么区别，只是把分母变成了Zw，从而看起来简化不少。那么我们已经得到最优的P(y|x)了，将其代入L式中再求最优的w不就好了嘛，这个步骤也就是求个导，书上不详细展开了，因为是重复性工作，因此在这里我觉得也没有必要再讲求导式写出来，读者可以自己写一下。

因此综上，只要能够得到最优的Pw(y|x)，随后再求得最优的w，全部工作也就迎刃而解了。最主要的是Pw(y|x)问题。对于如何求解最优的Pw(y|x)，这里介绍改进的迭代尺度法(IIS)

改进的迭代尺度法(IIS)

我们已知目前的解决目标是：
所有的式子连乘取对数转换为对数似然函数为：

（在《统计学习方法》6.2.4节中说明了为何上文中关于w进行拉格朗日极大化等价于对数似然函数最大化，但由于只是一些简单的公式，这里不再展开。唯一难点在于书中给出的对数似然形式，可以查看该链接学习：最大熵模型中的对数似然函数的解释）

那么我们要干啥？当然得保证差值大于0呀。只有每次都比前一次的大才能保证不停地往上走。除此之外呢？我们还希望每次增加的量都尽可能大一点，这样才能以最快的速度收敛到极大值。利用不等式-lnα≥1-α我们可以得到下式：

A也可以称之为该变量的下届，就是对于任意的δ，它们的差值一定是大于等于A(δ|w)的。我们应该求一个δ使得尽量增大A(δ|w)，这样能保证我们每步都会走得尽可能跨度大，程序耗时也就越少。此外由于直接求A的最大值并不太好求，因此再将A松开一点，建立一个不太紧的方便求极大值的下界B，书中有详细公式推导，并不太难，因此不再讲解。给出最终的式子：
对B求导并令其为0，得到最优的δ，从而每次迭代过程中w=w+δ，以此提高L值，得到最优解。这就是算法的全过程。
补充（极其重要——针对mnist数据集）：

mnist是一个手写数字的数据集，里面有很多的样本，每个样本是28*28的图片，展开是一个784维的向量，代表一个数字，同时给了一个标签告诉你标签时多少。

好的，那么我的问题是：我们全文一直在说P(y|x)，可是你真正考虑过在mnist中里面的y和x都是什么吗？

有人说x是每一个样本，y是标签。看起来确实是这样的，那f(x, y)是什么？就是样本x和y是否成对出现过吗？那我如果出现一个训练集中没出现的样本，就变成0了吗？我们都知道手写每个人都不一样，不可能写出来的样本在训练集中一定有一个一模一样的，那么它就变成0了吗？

接下来我要说的事情，只针对Mnist数据集，作者进入机器学习时间尚短，不清楚书中是否是为了提高公式在所有情况下的泛化性而简略了写，因此对于其他方向的使用，保留一些意见。

事实上，对于Mnist而言所有的x都不是样本x！它表示的是样本的特征。以f(x, y)为例，或许应该写成f0(x, y)、f1(x, y), …., fn(x, y)，其中n为特征的数目。再详细下去以f0(x, y)为例，f0(x = 0, y = 0)表示是否存在这样一个事实：训练集所有样本的特征中有样本的第0个特征值为0，与此同时y为0。是不是有点理解了？

我们再以P(y|x)为例，例如P1(y=0|x=0)表示的是在训练集中所有样本中，第一个特征为0的样本其中y为0的概率。

我当时也是写程序的时候卡在这里总是不知道该怎么写，直到看了别人的实现以后才发现了这个没有提到的地方，关于具体实现可以查看我的代码，有详细注释。

最后还是需要补充一点，我并不太清楚在别的地方（例如NLP）中最大熵是如何应用的，也不清楚上面提到的这个小tip是否只是针对mnist而言做出的一个特例修改，因此保留对于这件事情的意见。

关于最大熵和朴素贝叶斯的联系

两者都使用了P(y|x)，区别在于朴素贝叶斯使用先验概率求解，最大熵利用最大化条件熵求解。

朴素贝叶斯倾向于从基本信息中提取新的信息，而最大熵将提取信息的程度进行了量化（就好像强迫自己获得概率一定是要均匀分布一样）。

# coding=utf-8
# Author:Dodo
# Date:2018-11-30
# Email:lvtengchao@pku.edu.cn
# Blog:www.pkudodo.com

'''
数据集：Mnist
训练集数量：60000(实际使用:20000)
测试集数量：10000
------------------------------
运行结果：
    正确率：96.9%
    运行时长：8.8h

备注：对于mnist而言，李航的统计学习方法中有一些关键细节没有阐述，
建议先阅读我的个人博客，其中有详细阐述。阅读结束后再看该程序。
Blog：www.pkudodo.com
'''

import time
import numpy as np
from collections import defaultdict

def loadData(fileName):
    '''
    加载Mnist数据集
    :param fileName:要加载的数据集路径
    :return: list形式的数据集及标记
    '''
    # 存放数据及标记的list
    dataList = []; labelList = []
    # 打开文件
    fr = open(fileName, 'r')
    # 将文件按行读取
    for line in fr.readlines():
        # 对每一行数据按切割福','进行切割，返回字段列表
        curLine = line.strip().split(',')
        #二分类，list中放置标签
        if int(curLine[0]) == 0:
            labelList.append(1)
        else:
            labelList.append(0)
        #二值化
        dataList.append([int(int(num) > 128) for num in curLine[1:]])

    #返回data和label
    return dataList, labelList


class maxEnt:
    '''
    最大熵类
    '''
    def __init__(self, trainDataList, trainLabelList, testDataList, testLabelList):
        '''
        各参数初始化
        '''
        self.trainDataList = trainDataList          #训练数据集
        self.trainLabelList = trainLabelList        #训练标签集
        self.testDataList = testDataList            #测试数据集
        self.testLabelList = testLabelList          #测试标签集
        self.featureNum = len(trainDataList[0])     #特征数量

        self.N = len(trainDataList)                 #总训练集长度
        self.n = 0                                  #训练集中（xi，y）对数量
        self.M = 10000                              #
        self.fixy = self.calc_fixy()                #所有(x, y)对出现的次数
        self.w = [0] * self.n                       #Pw(y|x)中的w
        self.xy2idDict, self.id2xyDict = self.createSearchDict()        #(x, y)->id和id->(x, y)的搜索字典
        self.Ep_xy = self.calcEp_xy()               #Ep_xy期望值

    def calcEpxy(self):
        '''
        计算特征函数f(x, y)关于模型P(Y|X)与经验分布P_(X, Y)的期望值（P后带下划线“_”表示P上方的横线
        程序中部分下划线表示“|”，部分表示上方横线，请根据具体公式自行判断,）
        即“6.2.2 最大熵模型的定义”中第二个期望（83页最上方的期望）
        :return:
        '''
        #初始化期望存放列表，对于每一个xy对都有一个期望
        #这里的x是单个的特征，不是一个样本的全部特征。例如x={x1，x2，x3.....，xk}，实际上是（x1，y），（x2，y），。。。
        #但是在存放过程中需要将不同特诊的分开存放，李航的书可能是为了公式的泛化性高一点，所以没有对这部分提及
        #具体可以看我的博客，里面有详细介绍  www.pkudodo.com
        Epxy = [0] * self.n
        #对于每一个样本进行遍历
        for i in range(self.N):
            #初始化公式中的P(y|x)列表
            Pwxy = [0] * 2
            #计算P(y = 0 } X)
            #注：程序中X表示是一个样本的全部特征，x表示单个特征，这里是全部特征的一个样本
            Pwxy[0] = self.calcPwy_x(self.trainDataList[i], 0)
            #计算P(y = 1 } X)
            Pwxy[1] = self.calcPwy_x(self.trainDataList[i], 1)

            for feature in range(self.featureNum):
                for y in range(2):
                    if (self.trainDataList[i][feature], y) in self.fixy[feature]:
                        id = self.xy2idDict[feature][(self.trainDataList[i][feature], y)]
                        Epxy[id] += (1 / self.N) * Pwxy[y]
        return Epxy

    def calcEp_xy(self):
        '''
        计算特征函数f(x, y)关于经验分布P_(x, y)的期望值（下划线表示P上方的横线，
        同理Ep_xy中的“_”也表示p上方的横线）
        即“6.2.2 最大熵的定义”中第一个期望（82页最下方那个式子）
        :return: 计算得到的Ep_xy
        '''
        #初始化Ep_xy列表，长度为n
        Ep_xy = [0] * self.n

        #遍历每一个特征
        for feature in range(self.featureNum):
            #遍历每个特征中的(x, y)对
            for (x, y) in self.fixy[feature]:
                #获得其id
                id = self.xy2idDict[feature][(x, y)]
                #将计算得到的Ep_xy写入对应的位置中
                #fixy中存放所有对在训练集中出现过的次数，处于训练集总长度N就是概率了
                Ep_xy[id] = self.fixy[feature][(x, y)] / self.N

        #返回期望
        return Ep_xy

    def createSearchDict(self):
        '''
        创建查询字典
        xy2idDict：通过(x,y)对找到其id,所有出现过的xy对都有一个id
        id2xyDict：通过id找到对应的(x,y)对
        '''
        #设置xy搜多id字典
        #这里的x指的是单个的特征，而不是某个样本，因此将特征存入字典时也需要存入这是第几个特征
        #这一信息，这是为了后续的方便，否则会乱套。
        #比如说一个样本X = (0, 1, 1) label =(1)
        #生成的标签对有(0, 1), (1, 1), (1, 1)，三个(x，y)对并不能判断属于哪个特征的，后续就没法往下写
        #不可能通过(1, 1)就能找到对应的id，因为对于(1, 1),字典中有多重映射
        #所以在生成字典的时总共生成了特征数个字典，例如在mnist中样本有784维特征，所以生成784个字典，属于
        #不同特征的xy存入不同特征内的字典中，使其不会混淆
        xy2idDict = [{} for i in range(self.featureNum)]
        #初始化id到xy对的字典。因为id与(x，y)的指向是唯一的，所以可以使用一个字典
        id2xyDict = {}

        #设置缩影，其实就是最后的id
        index = 0
        #对特征进行遍历
        for feature in range(self.featureNum):
            #对出现过的每一个(x, y)对进行遍历
            #fixy：内部存放特征数目个字典，对于遍历的每一个特征，单独读取对应字典内的(x, y)对
            for (x, y) in self.fixy[feature]:
                #将该(x, y)对存入字典中，要注意存入时通过[feature]指定了存入哪个特征内部的字典
                #同时将index作为该对的id号
                xy2idDict[feature][(x, y)] = index
                #同时在id->xy字典中写入id号，val为(x, y)对
                id2xyDict[index] = (x, y)
                #id加一
                index += 1

        #返回创建的两个字典
        return xy2idDict, id2xyDict


    def calc_fixy(self):
        '''
        计算(x, y)在训练集中出现过的次数
        :return:
        '''
        #建立特征数目个字典，属于不同特征的(x, y)对存入不同的字典中，保证不被混淆
        fixyDict = [defaultdict(int) for i in range(self.featureNum)]
        #遍历训练集中所有样本
        for i in range(len(self.trainDataList)):
            #遍历样本中所有特征
            for j in range(self.featureNum):
                #将出现过的(x, y)对放入字典中并计数值加1
                fixyDict[j][(self.trainDataList[i][j], self.trainLabelList[i])] += 1
        #对整个大字典进行计数，判断去重后还有多少(x, y)对，写入n
        for i in fixyDict:
            self.n += len(i)
        #返回大字典
        return fixyDict


    def calcPwy_x(self, X, y):
        '''
        计算“6.23 最大熵模型的学习” 式6.22
        :param X: 要计算的样本X（一个包含全部特征的样本）
        :param y: 该样本的标签
        :return: 计算得到的Pw(Y|X)
        '''
        #分子
        numerator = 0
        #分母
        Z = 0
        #对每个特征进行遍历
        for i in range(self.featureNum):
            #如果该(xi,y)对在训练集中出现过
            if (X[i], y) in self.xy2idDict[i]:
                #在xy->id字典中指定当前特征i，以及(x, y)对：(X[i], y)，读取其id
                index = self.xy2idDict[i][(X[i], y)]
                #分子是wi和fi(x，y)的连乘再求和，最后指数
                #由于当(x, y)存在时fi(x，y)为1，因为xy对肯定存在，所以直接就是1
                #对于分子来说，就是n个wi累加，最后再指数就可以了
                #因为有n个w，所以通过id将w与xy绑定，前文的两个搜索字典中的id就是用在这里
                numerator += self.w[index]
            #同时计算其他一种标签y时候的分子，下面的z并不是全部的分母，再加上上式的分子以后
            #才是完整的分母，即z = z + numerator
            if (X[i], 1-y) in self.xy2idDict[i]:
                #原理与上式相同
                index = self.xy2idDict[i][(X[i], 1-y)]
                Z += self.w[index]
        #计算分子的指数
        numerator = np.exp(numerator)
        #计算分母的z
        Z = np.exp(Z) + numerator
        #返回Pw(y|x)
        return numerator / Z


    def maxEntropyTrain(self, iter = 500):
        #设置迭代次数寻找最优解
        for i in range(iter):
            #单次迭代起始时间点
            iterStart = time.time()

            #计算“6.2.3 最大熵模型的学习”中的第二个期望（83页最上方哪个）
            Epxy = self.calcEpxy()

            #使用的是IIS，所以设置sigma列表
            sigmaList = [0] * self.n
            #对于所有的n进行一次遍历
            for j in range(self.n):
                #依据“6.3.1 改进的迭代尺度法” 式6.34计算
                sigmaList[j] = (1 / self.M) * np.log(self.Ep_xy[j] / Epxy[j])

            #按照算法6.1步骤二中的（b）更新w
            self.w = [self.w[i] + sigmaList[i] for i in range(self.n)]

            #单次迭代结束
            iterEnd = time.time()
            #打印运行时长信息
            print('iter:%d:%d, time:%d'%(i, iter, iterStart - iterEnd))

    def predict(self, X):
        '''
        预测标签
        :param X:要预测的样本
        :return: 预测值
        '''
        #因为y只有0和1，所有建立两个长度的概率列表
        result = [0] * 2
        #循环计算两个概率
        for i in range(2):
            #计算样本x的标签为i的概率
            result[i] = self.calcPwy_x(X, i)
        #返回标签
        #max(result)：找到result中最大的那个概率值
        #result.index(max(result))：通过最大的那个概率值再找到其索引，索引是0就返回0，1就返回1
        return result.index(max(result))

    def test(self):
        '''
        对测试集进行测试
        :return:
        '''
        #错误值计数
        errorCnt = 0
        #对测试集中所有样本进行遍历
        for i in range(len(self.testDataList)):
            #预测该样本对应的标签
            result = self.predict(self.testDataList[i])
            #如果错误，计数值加1
            if result != self.testLabelList[i]:   errorCnt += 1
        #返回准确率
        return 1 - errorCnt / len(self.testDataList)

if __name__ == '__main__':
    start = time.time()

    # 获取训练集及标签
    print('start read transSet')
    trainData, trainLabel = loadData('../Mnist/mnist_train.csv')

    # 获取测试集及标签
    print('start read testSet')
    testData, testLabel = loadData('../Mnist/mnist_test.csv')

    #初始化最大熵类
    maxEnt = maxEnt(trainData[:20000], trainLabel[:20000], testData, testLabel)

    #开始训练
    print('start to train')
    maxEnt.maxEntropyTrain()

    #开始测试
    print('start to test')
    accuracy = maxEnt.test()
    print('the accuracy is:', accuracy)

    # 打印时间
    print('time span:', time.time() - start)