前言：

        本章的内容上的理解并不困难，但最大熵模型的代码实现较为困难，其主要原因是李航老师的《统计学习方法》中关于特征函数fi(x,y)的解释为了追求泛化而让人感到迷惑。

        这里感谢pku的大佬Dodo提供的思路和解释，本文章在Dodo实现的最大熵模型的基础上做了略微的改进，使得该模型可以对多标签（类别大于2）数据进行分类。

        Dodo的博客：统计学习方法|最大熵原理剖析及实现 | Dodo

一、实现思路

1.1 文章所用的简单数据集

为了便于后文的解释，使用的简单训练数据集如下：

1.2 最大熵模型的基本形式

由课本内容可知，上图为该最大熵模型的具体形式。其作用是当给定测试数据X时，该模型输出每种类别的概率。

例如输入：X=['sunny', 'hot', 'high', 'FALSE']，输出：[('no', 0.9094), ('yes', 0.0905)]。

1.3 最大化目标函数

我们的目标是通过更新每个特征函数前面的拉格朗日乘子wi，使得对数似然函数达到最大。

本文使用的方法是改进迭代尺度算法（IIS）。

由于最大熵模型的参数w是一个序列（列表），其长度等于特征函数的数量n，所以在IIS中每次迭代更新时我们也需要一个σ序列，让w序列中的每个wi = wi+σi。

所以现在我们需要在每次迭代中得到一个σ序列，而计算σ序列中的每个σi都需要三个值，

分别是M、fi(x,y)的经验期望、fi(x,y)的期望。

1.3.1 关于特征函数fi(x,y)的解释

实际上特征函数真正的形式并非是fi(X,y) = fi( X = ['sunny', 'hot', 'high', 'FALSE'], y = ['no'] )的样子。而是对于该训练数据(X, y)生成的特征函数如下：

f1(x1,y) = (['sunny'] , ['no'])

f2(x2,y) = (['hot'] , ['no'])

f3(x3,y) = (['high'] , ['no'])

f4(x4,y) = (['FALSE'] , ['no'])

同样书中的P_(x)和P_(x, y)也是针对单一特征而言，例如P_(x,y) = P_(['sunny'] , ['no'] )。

1.3.2 关于IIS迭代过程的细节

其实从IIS算法求σ序列的过程我们可以发现一个细节，那就是在计算σi的过程中，log分子上的经验期望实际上是一成不变的，这是因为fi(x,y)的经验期望公式中没有涉及到最大熵模型P(y | x)，而分母的fi(x,y)期望公式由于涉及到最大熵模型，所以每次迭代都需要重新计算fi(x,y)期望公式（因为上次迭代之后最大熵模型发生了更新）。

1.4 部分代码解释

1.4.1 fit函数

（1）首先我们给定训练数据和训练数据的标签。

（2）对训练数据和标签进行初始化，得到我们需要的各种变量。

（3）在指定迭代次数之内，按照IIS步骤进行：先求特征函数的期望(非经验)Epxy，再生成σ序列（利用M，fixy的经验期望，fixy的期望计算），然后利用σ序列对参数w进行更新。

 def fit(self, x_train, y_train):
        self.init_data(x_train, y_train)
        for i in range(self.max_iter):
            Epxy = self.cal_Epxy()

            sigmaList = [0] * self.fixy_num  # IIS算法所需要的σ更新变量表

            for j in range(self.fixy_num):
                sigmaList[j] = (1 / self.fixy_num) * np.log(self.Ep_xy[j] / Epxy[j])
            # 参数w的更新
            self.w = [self.w[i] + sigmaList[i] for i in range(self.fixy_num)]
            # 完成单次迭代

        # 训练结束

1.4.2 数据初始化

class MaxEntropy:
    def __init__(self, max_iter):
        self.train_data = None  # 训练数据
        self.train_label = None  # 训练标签
        self.label_set = None  # 标签y的set
        self.feature_num = None  # 特征数目，即训练数据x有多少列
        self.data_len = None  # 数据集长度
        self.fixy = None  # 一个list，针对数据的每一个特征都有一个字典，里面记录了该特征下具体的特征值和标签(xi,y)出现的次数
        self.fixy_num = 0  # 特征函数fixy的数量
        self.w = None  # 特征函数前面的拉格朗日乘子
        self.x_y2idx = None  # 双向索引
        self.idx2x_y = None  # 双向索引
        self.Ep_xy = None  # EP_xy的期望值EP_xy=sum_xy( P_(x,y) * f(x,y) )
        self.max_iter = max_iter  # 最大迭代次数

    def init_data(self, x_train, y_train):
        self.train_data = x_train
        self.train_label = y_train
        self.data_len = len(y_train)
        self.feature_num = x_train.shape[1]
        self.label_set = set(y_train)

        # fixyDict是一个list里面为每列特征都创建一个字典
        fixyDict = [defaultdict(int) for i in range(self.feature_num)]
        for i in range(self.data_len):
            for j in range(self.feature_num):
                fixyDict[j][(self.train_data[i][j], self.train_label[i])] += 1
        for i in fixyDict:
            self.fixy_num += len(i)
        self.fixy = fixyDict
        self.w = [0] * self.fixy_num

        # 双向索引绑定
        # 为了通过(x,y)能找到对应是第几列的特征，对每个特征创造一个字典
        self.x_y2idx = [{} for i in range(self.feature_num)]
        self.idx2x_y = {}
        index = 0
        for i in range(self.feature_num):  # i就是第几列特征的列索引
            for (x, y) in self.fixy[i]:  # 从大字典中取出第i列特征对应的特征函数fixy,直接对字典遍历得到的(x,y)是key值
                # 双向索引，为每一列特征的每一个特征函数赋予一个index
                self.x_y2idx[i][(x, y)] = index  # 通过第i列的特征函数(x,y)可以找到一个id
                self.idx2x_y[index] = (x, y)  # 通过id也能找到具体的特征函数(x, y)
                index += 1
        # 计算经验期望值
        self.Ep_xy = self.cal_Ep_xy()

1.4.3 最大熵模型函数

    def cal_Pwy_x(self, X, y):  # 此处传入的X是一个完整的训练数据(包含所有的特征(x1,x2,x3...,y)
        numerator = 0  # 分子
        Z = 0  # 分母 也就是规范化因子

        for i in range(self.feature_num):  # 求分子
            if (X[i], y) in self.fixy[i]:  # 如果存在(xi, y)的特征函数
                idx = self.x_y2idx[i][(X[i], y)]  # 该idx不光对照了每个特征函数的期望，还对照了特征函数前面的拉格朗日乘子wi
                numerator += self.w[idx]
        numerator = np.exp(numerator)

        for label in self.label_set:  # 求分母
            zi = 0
            for i in range(self.feature_num):
                if (X[i], label) in self.fixy[i]:
                    idx = self.x_y2idx[i][(X[i], label)]
                    zi += self.w[idx]
            Z += np.exp(zi)

        return numerator / Z

1.4.4 计算特征函数fi(x,y)的经验期望/非经验期望

    # 计算每个特征函数的经验期望
    def cal_Ep_xy(self):
        Ep_xy = [0] * self.fixy_num  # 一共有多少特征函数就有多少个期望
        for feature in range(self.feature_num):
            # 课本上fixy的期望公式太泛化，实际上fixy特征函数的期望应该是：[fi出现的次数/总数据量]（其实是P_(x,y)）* fi
            for (x, y) in self.fixy[feature]:
                idx = self.x_y2idx[feature][(x, y)]
                # 课本上的特征函数期望遍历的x,y,但是实际上只有（x,y）满足了fi，这个期望值才存在
                Ep_xy[idx] = self.fixy[feature][(x, y)] / self.data_len
        return Ep_xy

    # 计算每个特征函数的期望(非经验期望)
    def cal_Epxy(self):
        Epxy = [0] * self.fixy_num

        for i in range(self.data_len):
            Pwxy = {}
            for label in self.label_set:
                Pwxy[label] = self.cal_Pwy_x(self.train_data[i], label)

            for feature in range(self.feature_num):
                for label in self.label_set:
                    if (self.train_data[i][feature], label) in self.fixy[feature]:
                        idx = self.x_y2idx[feature][(self.train_data[i][feature], label)]
                        Epxy[idx] += Pwxy[label] * (1 / self.data_len)  # P_(x)是取1/N
        return Epxy

1.4.5 预测新数据

    def predict(self, x_test):
        result = []
        for label in self.label_set:
            result.append((label, self.cal_Pwy_x(x_test, label)))  # 结果里存放每一个预测的类别和对应的概率
        return result

1.5 完整代码

import numpy as np
from _collections import defaultdict

class MaxEntropy:
    def __init__(self, max_iter):
        self.train_data = None  # 训练数据
        self.train_label = None  # 训练标签
        self.label_set = None  # 标签y的set
        self.feature_num = None  # 特征数目，即训练数据x有多少列
        self.data_len = None  # 数据集长度
        self.fixy = None  # 一个list，针对数据的每一个特征都有一个字典，里面记录了该特征下具体的特征值和标签(xi,y)出现的次数
        self.fixy_num = 0  # 特征函数fixy的数量
        self.w = None  # 特征函数前面的拉格朗日乘子
        self.x_y2idx = None  # 双向索引
        self.idx2x_y = None  # 双向索引
        self.Ep_xy = None  # EP_xy的期望值EP_xy=sum_xy( P_(x,y) * f(x,y) )
        self.max_iter = max_iter  # 最大迭代次数

    def cal_Pwy_x(self, X, y):  # 此处传入的X是一个完整的训练数据(包含所有的特征(x1,x2,x3...,y)
        numerator = 0  # 分子
        Z = 0  # 分母 也就是规范化因子

        for i in range(self.feature_num):  # 求分子
            if (X[i], y) in self.fixy[i]:  # 如果存在(xi, y)的特征函数
                idx = self.x_y2idx[i][(X[i], y)]  # 该idx不光对照了每个特征函数的期望，还对照了特征函数前面的拉格朗日乘子wi
                numerator += self.w[idx]
        numerator = np.exp(numerator)

        for label in self.label_set:  # 求分母
            zi = 0
            for i in range(self.feature_num):
                if (X[i], label) in self.fixy[i]:
                    idx = self.x_y2idx[i][(X[i], label)]
                    zi += self.w[idx]
            Z += np.exp(zi)

        return numerator / Z

    # 计算每个特征函数的经验期望
    def cal_Ep_xy(self):
        Ep_xy = [0] * self.fixy_num  # 一共有多少特征函数就有多少个期望
        for feature in range(self.feature_num):
            # 课本上fixy的期望公式太泛化，实际上fixy特征函数的期望应该是：[fi出现的次数/总数据量]（其实是P_(x,y)）* fi
            for (x, y) in self.fixy[feature]:
                idx = self.x_y2idx[feature][(x, y)]
                # 课本上的特征函数期望遍历的x,y,但是实际上只有（x,y）满足了fi，这个期望值才存在
                Ep_xy[idx] = self.fixy[feature][(x, y)] / self.data_len
        return Ep_xy

    # 计算每个特征函数的期望(非经验期望)
    def cal_Epxy(self):
        Epxy = [0] * self.fixy_num

        for i in range(self.data_len):
            Pwxy = {}
            for label in self.label_set:
                Pwxy[label] = self.cal_Pwy_x(self.train_data[i], label)

            for feature in range(self.feature_num):
                for label in self.label_set:
                    if (self.train_data[i][feature], label) in self.fixy[feature]:
                        idx = self.x_y2idx[feature][(self.train_data[i][feature], label)]
                        Epxy[idx] += Pwxy[label] * (1 / self.data_len)  # P_(x)是取1/N
        return Epxy

    def init_data(self, x_train, y_train):
        self.train_data = x_train
        self.train_label = y_train
        self.data_len = len(y_train)
        self.feature_num = x_train.shape[1]
        self.label_set = set(y_train)

        # fixyDict是一个list里面为每列特征都创建一个字典
        fixyDict = [defaultdict(int) for i in range(self.feature_num)]
        for i in range(self.data_len):
            for j in range(self.feature_num):
                fixyDict[j][(self.train_data[i][j], self.train_label[i])] += 1
        for i in fixyDict:
            self.fixy_num += len(i)
        self.fixy = fixyDict
        self.w = [0] * self.fixy_num

        # 双向索引绑定
        # 为了通过(x,y)能找到对应是第几列的特征，对每个特征创造一个字典
        self.x_y2idx = [{} for i in range(self.feature_num)]
        self.idx2x_y = {}
        index = 0
        for i in range(self.feature_num):  # i就是第几列特征的列索引
            for (x, y) in self.fixy[i]:  # 从大字典中取出第i列特征对应的特征函数fixy,直接对字典遍历得到的(x,y)是key值
                # 双向索引，为每一列特征的每一个特征函数赋予一个index
                self.x_y2idx[i][(x, y)] = index  # 通过第i列的特征函数(x,y)可以找到一个id
                self.idx2x_y[index] = (x, y)  # 通过id也能找到具体的特征函数(x, y)
                index += 1
        # 计算经验期望值
        self.Ep_xy = self.cal_Ep_xy()

    def fit(self, x_train, y_train):
        self.init_data(x_train, y_train)
        for i in range(self.max_iter):
            Epxy = self.cal_Epxy()

            sigmaList = [0] * self.fixy_num  # IIS算法所需要的σ更新变量表

            for j in range(self.fixy_num):
                sigmaList[j] = (1 / self.fixy_num) * np.log(self.Ep_xy[j] / Epxy[j])
            # 参数w的更新
            self.w = [self.w[i] + sigmaList[i] for i in range(self.fixy_num)]
            # 完成单次迭代

        # 训练结束

    def predict(self, x_test):
        result = []
        for label in self.label_set:
            result.append((label, self.cal_Pwy_x(x_test, label)))  # 结果里存放每一个预测的类别和对应的概率
        return result


if __name__ == '__main__':
    # 第一列是标签(代表是否要外出游玩)
    # 第二至第五列分别为特征：天气状态、气温状态、湿度水平、是否有风
    data_set = [['no', 'sunny', 'hot', 'high', 'FALSE'],
                ['no', 'sunny', 'hot', 'high', 'TRUE'],
                ['yes', 'overcast', 'hot', 'high', 'FALSE'],
                ['yes', 'rainy', 'mild', 'high', 'FALSE'],
                ['yes', 'rainy', 'cool', 'normal', 'FALSE'],
                ['no', 'rainy', 'cool', 'normal', 'TRUE'],
                ['yes', 'overcast', 'cool', 'normal', 'TRUE'],
                ['no', 'sunny', 'mild', 'high', 'FALSE'],
                ['yes', 'sunny', 'cool', 'normal', 'FALSE'],
                ['yes', 'rainy', 'mild', 'normal', 'FALSE'],
                ['yes', 'sunny', 'mild', 'normal', 'TRUE'],
                ['yes', 'overcast', 'mild', 'high', 'TRUE'],
                ['yes', 'overcast', 'hot', 'normal', 'FALSE'],
                ['no', 'rainy', 'mild', 'high', 'TRUE']]
    data_set = np.array(data_set)
    x_train, y_train = data_set[:, 1:], data_set[:, 0]
    clf = MaxEntropy(max_iter=200)
    clf.fit(x_train, y_train)
    predict_1 = clf.predict(['sunny', 'hot', 'high', 'FALSE'])
    predict_2 = clf.predict(['sunny', 'hot', 'high', 'TRUE'])
    predict_3 = clf.predict(['overcast', 'hot', 'high', 'FALSE'])
    print(predict_1)
    print(predict_2)
    print(predict_3)

1.6 运行效果

由于数据比较匮乏，我们直接使用训练数据的前三个数据进行效果测试。

即使增加标签种类也可以预测。