第5章 CHAPTER 5 逻辑回归与最大熵模型 逻辑回归的核心思想是把原线性回归的取值范围通过Logistic函数映射到一个概率空间,从而将一个回归模型转换成了一个分类模型。本章将从最简单的二元逻辑回归出发,逐渐拓展到多元逻辑回归。非常重要的一点是,多元逻辑回归同时也是人工神经网络中的Softmax模型,或者从信息论的角度来解释,它又被称为最大熵模型。在推导最大熵模型时,还会用到凸优化(包含拉格朗日乘数法)的内容。 5.1逻辑回归 下面通过一个例子来引出逻辑回归的基本思想。为研究与急性心肌梗死急诊治疗预后情况有关的因素,现收集了200个急性心肌梗死的病例如表51所示。其中,x1用于指示救治前是否休克,x1=1表示救治前已休克,x1=0表示救治前未休克; x2用于指示救治前是否心衰,x2=1表示救治前已发生心衰,x2=0表示救治前未发生心衰; x3用于指示12小时内有无治疗措施,x3=1表示没有,x3=0表示有。最后P给出了患者的最终结局,当P=0时,表示患者生存; 当P=1时,表示患者死亡。 表51急性心肌梗死的病例数据 P=0 P=1 x1 x2 x3 N x1 x2 x3 N 0 0 0 35 0 0 0 4 0 0 1 34 0 0 1 10 0 1 0 17 0 1 0 4 0 1 1 19 0 1 1 15 1 0 0 17 1 0 0 6 1 0 1 6 1 0 1 9 1 1 0 6 1 1 0 6 1 1 1 6 1 1 1 6 如果要建立回归模型,进而预测不同情况下患者生存的概率,此时P的取值应该为0~1。考虑用线性回归来建模: P=w0+w1x1+w2x2+w3x3 显然将自变量代入上述回归方程,并不能保证概率P一定位于0~1。最直接的想法是使用某个函数将上述线性回归模型的预测值映射至0~1,逻辑回归模型就是使用Logistic函数完成这件事的。Logistic函数的定义如下: P=ez1+ez=11+e-z 或 z=lnP1-P 图51Logistic函数图形 其函数图像如图51所示。 用上面的Logistic函数定义式将多元线性回归中的因变量替换得到 lnP1-P=w0+w1x1+…+wnxn 或者 P=ew0+w1x1+…+wnxn1+ew0+w1x1+…+wnxn 或者 P=11+e-(w0+w1x1+…+wnxn) 可以简记为 P(z)=11+e-z 其中,z=w0+w1x1+…+wnxn,而且在当前所讨论的例子中n=3。 上式中的w0,w1,…,wn正是需要求解的参数,通常采用极大似然估计法对参数进行求解。对于每个观测到的样本,使用变量y标记样本的最终状态(0为生存,1为死亡)。采用上述模型可以简单得到每个样本生存和死亡的概率分别如下: P(y=1)=ez1+ez,P(y=0)=1-P(y=1)=11+ez 进而计算似然概率有 L=exp(w0)1+exp(w0)3511+exp(w0)4× exp(w0+w3)1+exp(w0+w3)3411+exp(w0+w3)10× … exp(w0+w1+w2+w3)1+exp(w0+w1+w2+w3)611+exp(w0+w1+w2+w3)6 =∏ki=1exp(w0+w1x1i+w2x2i+w3x3i)1+exp(w0+w1x1i+w2x2i+w3x3i)ni011+exp(w0+w1x1i+w2x2i+w3x3i)ni1 寻找最适宜的w^0,w^1,w^2,w^3使得L达到最大,最终得到估计模型为 lnP1-P=-2.0858+1.1098x1+0.7028x2+0.9751x3 或写成 P=11+e-(-2.0858+1.1098x1+0.7028x2+0.9751x3) 得到上面的公式后,如果再有一组观察样本,将其代入公式,就可以算得患者生存与否的概率。例如,现在有一名患者A没有休克,病发5小时后送医院,而且已出现了症状,即x1=0,x2=1,x3=0,则可据此计算其生存的概率为 PA=11+e-(-2.0858+0.7028)=0.200 同理,若另有一名患者B已经出现休克,病发18小时后送医院,出现了症状,即x1=1,x2=1,x3=1,则可据此计算其生存的概率为 PB=11+e-(-2.0858+1.1098+0.7028+0.9751)=0.669 前面在对参数进行估计时,我们其实假设了y=1的概率为 P(y=1|x; w)=hw(x)=g(wTx)=ewTx1+ewTx 于是还有 P(y=0|x; w)=1-hw(x)=1-g(wTx)=11+ewTx 由此,便可以利用逻辑回归从特征学习中得出一个非0即1的分类模型。当要判别一个新来的特征属于哪个类时,只需求hw(x)即可,若hw(x)大于0.5就可被归为y=1的类; 反之就被归为y=0类。 以上通过一个例子向读者演示了如何从原始的线性回归演化出逻辑回归。而且不难发现,逻辑回归可以用作机器学习中的分类器。当我们得到一个事件发生与否的概率时,自然就已经得出结论,其到底应该属于“发生”的那一类别,还是属于“不发生”的那一类别。 机器学习最终是希望机器自己学到一个可以用于问题解决的模型。而这个模型本质上是由一组参数(或权值)定义的,也就是前面讨论的w0,w1,…,wn。在得到测试数据时,将这组参数与测试数据线性加和得到 z=w0+w1x1+…+wnxn 这里x1,x2,…,xn是每个样本的n个特征。之后再按照Logistic函数的形式求出 P(z)=11+e-z 在给定特征向量x=(x1,x2,…,xn)时,条件概率P(y=1|x)为根据观测量某事件y发生的概率。那么逻辑回归模型可以表示为 P(y=1|x)=π(x)=11+e-z 相对应地,在给定条件x时,事件y不发生的概率为 P(y=0|x)=1-π(x)=11+ez 而且还可以得到事件发生与不发生的概率之比为 odds=P(y=1|x)P(y=0|x)=ez 这个比值称为事件的发生比。 概率论的知识告诉我们进行参数估计时可以采用最大似然法。假设有m个观测样本,观测值分别为y1,y2,…,yn,设pi=P(yi=1|xi)为给定条件下得到yi=1的概率。同样地,yi=0的概率为1-pi=P(yi=0|xi),所以得到一个观测值的概率为P(yi)=pyii(1-pi)1-yi。 各个观测样本之间相互独立,那么它们的联合分布为各边缘分布的乘积。得到似然函数为 L(w)=∏mi=1[π(xi)]yi[1-π(xi)]1-yi 我们的目标是求出使这一似然函数值最大的参数估计,于是对函数取对数得到 lnL(w)=∑mi=1{yiln[π(xi)]+(1-yi)ln[1-π(xi)]} =∑mi=1ln[1-π(xi)]+∑mi=1yilnπ(xi)1-π(xi) =∑mi=1ln[1-π(xi)]+∑mi=1yizi =∑mi=1-ln[1+ezi]+∑mi=1yizi 其中,zi=w·xi=w0+w1xi1+…+wnxin,根据多元函数求极值的方法,为了求出使得lnL(w)最大的向量w=(w0,w1,…,wn),对上述的似然函数求偏导后得到 lnL(w)wk=∑mi=1yixik-∑mi=111+eziezi·xik =∑mi=1xik[yi-π(xi)] 现在,我们所要做的就是通过上面已经得到的结论来求解使得似然函数最大化的参数向量。在实际中有很多方法可供选择,其中比较常用的有梯度下降法、牛顿法和拟牛顿法等。 5.2牛顿法解逻辑回归 现代计算中涉及大量的工程计算问题,这些计算问题往往很少采用我们通常在求解计算题甚至是考试时所采用的方法,因为计算机最擅长的无非就是“重复执行大量的简单任务”,所以数值计算方面的迭代法在计算机时代便有了很大的作用。一个典型的例子就是利用牛顿迭代法近似求解方程的方法。牛顿迭代又称为牛顿拉夫逊(NewtonRaphson)方法。 有时候某些方程的求根公式可能很复杂(甚至有些方程可能没有求根公式),导致求解困难,这时便可利用牛顿法进行迭代求解。 图52牛顿迭代法求方程的根 如图52所示,假设要求解方程f(x)=0的根,首先随便找一个初始值x0,如果x0不是解,做一个经过(x0,f(x0))这个点的切线,与x轴的交点为x1。同样的道理,如果x1不是解,做一个经过(x1,f(x1))这个点的切线,与x轴的交点为x2。以此类推。以这样的方式得到的xi会无限趋近于f(x)=0的解。 判断xi是否是f(x)=0的解有两种方法: 一是直接计算f(xi)的值并判断其是否为0,二是判断前后两个解xi和xi-1是否无限接近。经过(xi,f(xi))这个点的切线方程为(注意这也是一元函数的一阶泰勒展开式) f(x)=f(xi)+f′(xi)(x-xi) 其中,f′(x)为f(x)的导数。令切线方程等于 0,即可求出 xi+1=xi-f(xi)f′(xi) 于是得到了一个迭代公式,而且它必然在f(x*)=0处收敛,其中x*就是方程的根,由此便可对方程进行迭代求根。 接下来讨论牛顿法在最优化问题中的应用。假设当前任务是优化一个目标函数f,也就是求该函数的极大值或极小值问题,可以转化为求解函数f的导数f′=0的问题,这样就可以把优化问题看成方程f′=0的求解问题。剩下的问题就和前面提到的牛顿迭代法求解很相似了。 这次为了求解方程f′=0的根,把原函数f(x)做泰勒展开,展开到二阶形式(注意之前是一阶) f(x+Δx)=f(x)+f′(x)Δx+12f″(x)Δx2 当且仅当Δx无限趋近于0时,(可以舍去后面的无穷小项)使得等式成立。此时,上式等价于 f′(x)+12f″(x)Δx=0 注意因为Δx无限趋近于0,前面的的常数12将不再起作用,可以将其一并忽略,即 f′(x)+f″(x)Δx=0 求解 Δx=-f′(x)f″(x) 得出迭代公式 xn+1=xn-f′(x)f″(x),n=0,1,… 最优化问题除了用牛顿法来解之外,还可以用梯度下降法来解。但是通常来说,牛顿法可以利用到曲线本身的信息,比梯度下降法更容易收敛,即迭代更少次数。 再次联系到Hessian矩阵与多元函数极值,对于一个多维向量x,以及在点x0的邻域内有连续二阶偏导数的多元函数f(x),可以写出该函数在点x0处的(二阶)泰勒展开式 f(x)=f(x0)+(x-x0)Tf(x0)+12(x-x0)THf(x0)(x-x0)+o(‖x-x0‖2) 其中,o(‖x-x0‖2)是高阶无穷小表示的皮亚诺余项。而Hf(x0)是一个Hessian矩阵。依据之前的思路,忽略掉无穷小项,写出迭代公式即为 xn+1=xn-f(xn)Hf(xn),n≥0 由此可知,高维情况依然可以用牛顿迭代求解,但是问题是Hessian矩阵引入的复杂性,使得牛顿迭代求解的难度大大增加。所以人们又提出了所谓的拟牛顿法(QuasiNewton method),不再直接计算Hessian矩阵,而是每一步的时候使用梯度向量更新Hessian矩阵的近似。 现在尝试用牛顿法来求解逻辑回归。我们已经求出了逻辑回归的似然函数的一阶偏导数 lnL(w)wk=∑mi=1xik[yi-π(xi)] 由于lnL(w)是一个多元函数,变量是w=(w0,w1,…,wn),所以根据多元函数求极值问题的规则,易知极值点处的导数一定均为零,所以一共需要列出n+1个方程,联立解出所有的参数。下面列出方程组如下: lnL(w)w0=∑mi=1xi0[yi-π(xi)]=0 lnL(w)w1=∑mi=1xi1[yi-π(xi)]=0 lnL(w)w2=∑mi=1xi2[yi-π(xi)]=0  lnL(w)wn=∑mi=1xin[yi-π(xi)]=0 当然,在具体解方程组之前需要用Hessian矩阵来判断极值的存在性。求Hessian矩阵就得先求二阶偏导,即 2lnL(w)wkwj=∑mi=1xik[yi-π(xi)]wj =∑mi=1xik-ezi1+eziwj =-∑mi=1xik(1+e-zi)-1wj =-∑mi=1xik-(1+e-zi)-2·e-zi·-ziwj =-∑mi=1xik[(1+e-zi)-2·e-zi·xij] =-∑mi=1xik·11+ezi·ezi1+ezi·xij =∑mi=1xik·π(xi)·[π(xi)-1]·xij 显然可以用Hessian矩阵来表示以上多元函数的二阶偏导数,于是有 X=x10x11…x1n x20x21…x2n  xm0xm1…xmn A=π(x1)·[π(x1)-1]0…0 0π(x2)·[π(x2)-1]…0  00…π(xm)·[π(xm)-1] 所以得到Hessian矩阵H=XTAX,可以看出矩阵A是负定的。线性代数的知识告诉我们,如果A是负定的,那么Hessian矩阵H也是负定的。也就是说,多元函数存在局部极大值,这刚好与我们要求的最大似然估计相吻合。于是确信可以用牛顿迭代法来继续求解最优化问题。 对于多元函数求解零点,同样可以用牛顿迭代法,对于当前讨论的逻辑回归,可以得到如下迭代式 Xnew=Xold-UH=Xold-H-1U 其中H是Hessian矩阵,U的表达式如下: U=x10x11…x1n x20x21…x2n  xm0xm1…xmny1-π(x1) y2-π(x2)  ym-π(xm) 由于Hessian矩阵H是对称负定的,将矩阵A提取一个负号出来,得到 A′=π(x1)·[1-π(x1)]0…0 0π(x2)·[1-π(x2)]…0  00…π(xm)·[1-π(xm)] 然后Hessian矩阵H变为H′=XTA′X,这样H′就是对称正定的了。那么现在牛顿迭代公式变为 Xnew=Xold+(H′)-1U 现在我们需要考虑如何快速地算得(H′)-1U,即解方程组(H′)X=U。通常的做法是直接用高斯消元法求解,但是这种方法的效率一般比较低。而当前我们可以利用的一个有利条件是H′是对称正定的,所以可以用Cholesky矩阵分解法来解。 到这里,牛顿迭代法求解逻辑回归的原理已经介绍完了,但正如前面所提过的,在这个过程中因为要对Hessian矩阵求逆,计算量还是很大。于是研究人员又提出了拟牛顿法,它是针对牛顿法的弱点进行了改进,更具实践应用价值。 5.3应用实例: 二分类问题 本节将基于Python的机器学习包scikitlearn,给出一个使用逻辑回归模型处理二分类问题的实例。在此基础之上,会使用matplotlib和seaborn这两个常用的包对数据进行一定程度的可视化,方便读者了解数据理解模型。 5.3.1数据初探 一组来自世界银行的数据统计了30个国家的两项指标,第三产业增加值占GDP的比重和年龄大于或等于65岁的人口(也就是老龄人口)占总人口的比重。使用Pandas包的read_csv()方法,以DataFrame的格式读取数据。简单地打印前几行数据如下: import pandas as pd # 读取数据 data = pd.read_csv("countries_data.csv") # 以下为打印内容 countriesServices_of_GDPages65_of_totallabel 0Belgium76.7 181 1 France 78.9 18 1 2 Denmark 76.2 18 1 3 Spain 73.9 18 1 4 Japan 72.6 25 1 数据集包含四列: 第一列是国家名,第二列是“第三产业增加值占GDP比重”,第三列是“老龄人口占总人口的比重”,最后一列是数据的标注。简单根据数据集提供的两个特征绘制散点图,可以发现两个类别的数据之间区分度还是很明显的。此处引入matplotlib包,并使用seaborn包提供的scatterplot()函数,对数据集进行可视化。代码如下: import matplotlib.pyplot as plt import seaborn as sns # 绘制散点图 sns.scatterplot(data['Services_of_GDP'],data['ages65_of_total'],hue=data['label'],style=data['label'],data=data) # 打印绘制的散点图 plt.show() 上述代码中scatterplot()函数接收了5个参数,前两个参数分别指定了散点图的x轴和y轴所使用的数据; 第三和第四个参数分别指定了用以区分点的颜色和样式所依据的数据,此处使用label字段区分点的颜色和样式; 最后一个参数data指定绘图依赖的DataFrame,即data。 观察图53不难发现,数据集只包含两种标注0和1。右上方橙色的数据点第三产业增加值占比和老龄化人口占比都较高,应该是发达国家; 左下方蓝色的数据点,第三产业增加值占比和老龄化人口占比相对较低,应该是发展中国家。事实上,右边数据点代表的国家都是OECD(经济合作与发展组织)成员国,而左边数据点代表的国家则是发展中国家(包括一些东盟国家、南亚国家和拉美国家)。 5.3.2建模 在了解数据的大致分布后,使用scikitlearn下linear_model提供的LogisticRegression类直接根据数据集拟合一个逻辑回归模型。 图53彩图 图53数据集的简单可视化 from sklearn.linear_model import LogisticRegression # 区分特征和标注 X = data[['Services_of_GDP','Services_of_GDP']] Y = data['label'] # 训练模型 clf = LogisticRegression(penalty='l2',fit_intercept=True,solver='liblinear',max_iter=100).fit(X, Y) # 使用训练得到的模型预测训练数据集 prediction = clf.predict(X) 在上面的示例代码中,首先声明了一个LogisticRegression类并接收了4个参数。从左到右,第一个参数penalty指定模型的正则化方式,并提供l1正则、l2正则、elasticnet正则和none共4个选项,示例代码中使用了l2正则。第二个参数fit_intercept,配置决策函数中是否需要截距项,示例中设置为True表示需要截距项。第三个参数solver,用于指定模型的优化算法,有5个选项: newtoncg、lbfgs、liblinear、sag和saga。其中liblinear是基于C++的LIBLINEAR包实现的坐标下降法(coordinate descent),对于小数据集效果较好,上述示例中便使用了此优化算法。另外几种优化算法中,newtoncg是牛顿法的一种,lbfgs是拟牛顿法的一种,sag是随机平均梯度下降法,saga是sag的一种变体算法。其中newtoncg、lbfgs和sag算法,只能配合l1正则或无正则项; liblinear和saga均可配合l1和l2正则项,saga算法还能处理elasticnet正则项。对于小数据集的场景,liblinear、lbfgs和newtoncg都是不错的选择,而sag和saga则更适合大数据集的场景。第四个参数max_iter设置了优化算法的最大迭代次数,采用默认值100。更多模型参数,请查阅scikitlearn的官方文档。紧接着,使用了fit()方法,将训练数据X和标注数据Y传递给了方法,得到训练好的模型。 使用如下代码,将模型在训练集上的预测结果绘制在图54中。 # 将预测结果添加为data的一列 data['prediction'] = prediction # 使用颜色表示预测结果,样式表示真实标注 sns.scatterplot(data['Services_of_GDP'],data['ages65_of_total'],hue=data['prediction'],style=data['label'],data=data) plt.show() 图54彩图 图54模型分类结果 图54中使用颜色标识了模型的预测结果,而数据点的样式标识了数据的真实标注。圆点代表真实标注为“发展中国家”,十字代表真实标注为“发达国家”; 蓝色代表模型判断为“发展中国家”,橙色代表模型判断为“发达国家”。显然,一些发展中国家(圆点)被识别为了发达国家(橙色),即模型并未完全拟合训练集。若使用score()函数,输入待评估数据X和真实标注Y,则可以得到模型在训练数据集上的平均准确率。 # 使用score方法获取平均准确率 clf.score(X,Y) # 以下为打印结果 0.8666666666666667 为了让模型更好地拟合训练数据集,可以尝试增加模型的表征能力,即模型复杂度。具体到调参上,可以尝试增大模型的迭代次数; 通过参数C降低正则化项的权重; 还可以尝试使用收敛速度更快的优化算法。这个例子中,将参数C增大到2.5或者换用newtoncg和lbfgs算法均可在其他参数不变的情况下使模型的平均准确率达到1。感兴趣的读者可自行尝试。 此外scikitlearn还提供了另一种内置交叉验证的逻辑回归模型LogisticRegressionCV,其中CV代表交叉验证(crossvalidation)。这个方法在训练模型时会使用交叉验证根据score()方法的评估结果寻找最优的C和l1_ratio参数。使用同样的参数,这次换用LogisticRegressionCV训练模型并打印平均准确率。 from sklearn.linear_model import LogisticRegressionCV clf2 = LogisticRegressionCV(penalty='l2',fit_intercept=True,solver='liblinear',max_iter=100).fit(X, Y) clf2.score(X,Y) # 打印结果 1.0 换用LogisticRegressionCV建模后,几乎无须调参便完美拟合了训练数据集。实际应用中使用LogisticRegressionCV方法建模,可以在调参上节省不少的工作。此方法的参数列表和LogisticRegression方法大同小异,感兴趣的读者可自行查阅scikitlearn官方文档,了解更多细节。 5.4多元逻辑回归 在本章前面的介绍中,逻辑回归通常只用来处理二分类问题,所以又称这种逻辑回归为Binary(或Binomial)Logistic Regression。本节以此为基础,讨论如何处理多分类问题,这也就是通常所说的Multinomial Logistic Regression。MLR还有很多可能听起来更熟悉的名字,例如,在神经网络中,称其为Softmax; 从最大熵原理出发,又称其为最大熵模型。 人们之所以会提出逻辑回归的方法,其实是从线性回归中得到的启发。通常,线性回归可以表示为 P=w·x 其中,w是参数向量; x是样本观察值向量。对于一个二元分类器问题而言,我们希望P表示一个概率,也就是由x所表征的样本属于两个分类之一的概率。既然是概率,所以也就要求P介于0和1之间。但上述方程左侧的w·x可以取任何实数值。 为此,定义 odds=P1-P 为概率P对它的补1-P的比率,或正例对负例的比率。 然后对上式两边同时取对数,得到 logit 或 logodds如下: y=logit(P)=logP1-P=w·x 取反函数可得 P=ey1+ey 此时,概率P就介于0到1之间了。 在其他一些文献中,也有人从另外一个角度对逻辑回归进行了解读。既然我们最终是想求条件概率P(y|x),其中y是目标的类别,x是特征向量(或观察值向量)。那么根据贝叶斯定理,则有 P(y=1|x)=P(x|y=1)P(y=1)P(x|y=1)P(y=1)+P(x|y=0)P(y=0) =ey1+ey=11+e-y 其中 y=logP(x|y=1)P(y=1)P(x|y=0)P(y=0)=logP(y=1|x)P(y=0|x)=logP1-P 这与本书前面的解读是一致的。 逻辑回归遵循一个(单次的)二项分布,所以可以写出下列PMF: P(y|x)=Pw(x),y=1 1-Pw(x),y=0 其中,w是参数向量,且 Pw(x)=11+e-w·x 上述PMF也可写为如下形式: P(y|x)=[Pw(x)]y[1-Pw(x)]1-y 然后就可以使用N次观测样本的最大对数似然来对参数进行估计: lnL(w)=ln∏Ni=1[Pw(xi)]yi[1-Pw(xi)]1-yi =∑Ni=1{yilogPw(xi)+(1-yi)log[1-Pw(xi)]} 根据多元函数求极值的方法,为了求出使得lnL(w)最大的向量w=(w0,w1,…,wn),对上述的似然函数求偏导后得到 lnL(w)wk=wk∑Ni=1{yilogPw(xi)+(1-yi)log[1-Pw(xi)]} =∑Ni=1xik[yi-Pw(xi)] 因此,所要做的就是通过上面已经得到的结论来求解使得似然函数最大化的参数向量。我们已经在本章前面给出了基于牛顿迭代法求解的步骤。 现在考虑当目标类别数超过2时的情况,这时就需要使用多元逻辑回归(Multinomial Logistic Regression)。总的来说,  二元逻辑回归是一个针对二项分布样本的概率模型;  多元逻辑回归则将同样的原则扩展到多项分布的样本上。 在这个从二分类向多分类演进的过程中,参考之前的做法无疑会带给我们许多启示。由于参数向量w和特征向量x的线性组合w·x的取值范围是任意实数,为了能够引入概率值P,在二分类中我们的做法是将两种分类情况的概率之比(取对数后)作为与w·x相对应的取值,因为只有两个类别,所以类别1比上类别2也就等于类别1比上类别2的补,即 logP(y=1|x)P(y=0|x)=logP1-P=w·x 事实上,我们知道,如果P(y=1|x)大于P(y=0|x),自然就表示特征向量x所标准的样本属于y=1这一类别的概率更高。所以odds这个比值是相当有意义的。同理,如果现在的目标类别超过2个,即y=1,y=2,…,y=j,…,y=J,那么我们就应该从这些目标类别中找一个类别来作为基线,然后计算其他类别相对于这个基线的logodds。并让这个logodds作为线性函数的自变量。通常,我们选择最后一个类别作为基线。 在多元逻辑回归模型中,假设每一个响应(response)的logodds遵从一个线性模型,即 yj=logPjPJ=wj·x 其中,j=1,2,…,J-1。这个模型与之前的二元逻辑回归模型本质上是一致的。只是响应Y的分布由二项分布变成了多项分布。而且我们不再只有一个方程,而是J-1个。例如,当J=3时,我们有两个方程,分别比较j=1对j=3和j=2对j=3。你可能会疑惑是不是缺失了一个比较。缺失的那个比较是在类别1和类别2之间进行的,而它很容易从另外两个比较的结果中获得,因为log(P1/P2)=log(P1/P3)-log(P2/P3)。 多元逻辑回归也可以写成原始概率Pj(而不是logodds)的形式。从二项分布的逻辑回归中出发,有 y=logP(y=1)P(y=0)P=ey1+ey 类似地,在多项式分布的情况中,我们还可以得到 yj=logPjPJ=wj·xPj=eyj∑Jk=1eyk 其中,j=1,2,…,J。 下面我们来证明它。这里需要用到的一个事实是∑Jj=1Pj=1。首先证明 PJ=1∑Jj=1eyj 如下: PJ=PJ·∑Jj=1Pj∑Jj=1Pj=PJ(P1+P2+…+PJ)∑Jj=1Pj =P1P1+P2+…+PJPJ+P2P1+P2+…+PJPJ+…+PJP1+P2+…+PJPJ =P1ey1+ey2+…+eyJ+P2ey1+ey2+…+eyJ+…+PJey1+ey2+…+eyJ =P1+P2+…+PJey1+ey2+…+eyJ=∑Jj=1Pj∑Jj=1eyj=1∑Jj=1eyj 又因为Pj=PJ·eyj,而且yJ=0,综上便证明了之前的结论。 此外,如果从贝叶斯定理的角度出发,可得 P(y=j|x)=P(x|y=j)P(y=j)∑Jk=1P(x|y=k)P(y=k)=yj∑Jk=1eyk 其中,yk=log[P(x|y=k)P(y=k)]。 前面的推导就告诉我们,在多元逻辑回归中,响应变量y遵循一个多项分布。可以写出下列PMF(假设有J个类别): P(y|x)=ey1∑Jk=1eyk=ew1·x∑Jk=1ewk·x,y=1 ey2∑Jk=1eyk=ew2·x∑Jk=1ewk·x,y=2  eyj∑Jk=1eyk=ewj·x∑Jk=1ewk·x,y=j  eyJ∑Jk=1eyk=ewJ·x∑Jk=1ewk·x,y=J 二元逻辑回归和多元逻辑回归本质上是统一的,二元逻辑回归是多元逻辑回归的特殊情况。我们将这两类机器学习模型统称为逻辑回归。更进一步地,逻辑回归和“最大熵模型”本质上也是一致的、等同的。或者说最大熵模型是多元逻辑回归的另外一个称谓,我们将在5.5节讨论最大熵模型和多元逻辑回归的等同性。 5.5最大熵模型 我们在5.4节中导出的多元逻辑回归之一般形式为 P(y=j|x)=ewj·x∑Jk=1ewk·x 而本节将要介绍的最大熵模型的一般形式(其中的f为特征函数)为 Pw(y|x)=1Zw(x)exp∑Jk=1wkfk(x,y) 其中, Zw(x)=∑yexp∑Jk=1wkfk(x,y) 此处,x∈RJ,为输入,y={1,2,…,K}为输出,w=RJ,为权值向量,fk(x,y)为任意实值特征函数,其中k=1,2,…,J。 可见,多元逻辑回归和最大熵模型在形式上是统一的。事实上,尽管采用的方法不同,二者最终是殊途同归、万法归宗的。因此,无论是多元逻辑回归,还是最大熵模型,又或者是Softmax,它们本质上都是统一的。本节就将从最大熵原理这个角度来推导上述最大熵模型的一般形式。 5.5.1最大熵原理 本书前面已经讨论过熵的概念了。简单地说,假设离散随机变量X的概率分布是P(X),则其熵是 H(P)=-∑xP(x)logP(x) 而且熵满足下列不等式: 0≤H(P)≤log|X| 其中,|X|是X的取值个数,当且仅当X的分布是均匀分布时右边的等号成立。也就是说,当X服从均匀分布时,熵最大。 直观地,最大熵原理认为要选择的概率模型首先必须满足已有的事实,即约束条件。在没有更多信息的情况下,那些不确定的部分都是“等可能的”。最大熵原理通过熵的最大化来表示等可能性。“等可能性”不容易操作,而熵则是一个可以优化的数值指标。 吴军博士在其所著的《数学之美》一书中曾经谈到: “有一次,我去AT&T实验室作关最大熵模型的报告,随身带了一个骰子。我问听众‘每个面朝上的概率分别是多少’,所有人都说是等概率,即各种点数的概率均为1/6。这种猜测当然是对的。我问听众为什么,得到的回答是一致的: 对这个‘一无所知’的骰子,假定它每一面朝上概率均等是最安全的做法(你不应该主观假设它像韦小宝的骰子一样灌了铅)。从投资的角度看,就是风险最小的做法。从信息论的角度讲,就是保留了最大的不确定性,也就是说让熵达到最大。接着我又告诉听众,我的这个骰子被我特殊处理过,已知四点朝上的概率是1/3,在这种情况下,每个面朝上的概率是多少?这次,大部分人认为除去四点的概率是1/3,其余的均是2/15,也就是说,已知的条件(四点概率为1/3)必须满足,而对于其余各点的概率因为仍然无从知道,因此只好认为它们均等。注意,在猜测这两种不同情况下的概率分布时,大家都没有添加任何主观的假设,诸如四点的反面一定是三点等(事实上,有的骰子四点的反面不是三点而是一点)。这种基于直觉的猜测之所以准确,是因为它恰好符合了最大熵原理。” 通过上面这个关于骰子的例子,我们对最大熵原理应该已经有了一个基本的认识,即“对已有知识进行建模,对未知内容不做任何假设(model all that is known and assume nothing about that which is unknown)”。 5.5.2约束条件 最大熵原理是统计学习理论中的一般原理,将它应用到分类任务上就会得到最大熵模型。假设分类模型是一个条件概率分布P(Y|X),X∈InputRn表示输入(特征向量),Y∈Output表示输出(分类标签),Input和Output分别是输入和输出的集合。这个模型表示的是对于给定的输入X,输出为Y的概率是P(Y|X)。 给定一个训练数据集 T={(x1,y1),(x2,y2),…,(xN,yN)} 我们现在的目标是利用最大熵原理来选择最好的分类模型。首先考虑模型应该满足的条件。给定训练数据集,便可以据此确定联合分布P(X,Y)的经验分布P~(X,Y),以及边缘分布P(X)的经验分布。此处,则有 P~(X=x,Y=y)=v(X=x,Y=y)N P~(X=x)=v(X=x)N 其中,v(X=x,Y=y) 表示训练数据集中样本(x,y)出现的频率(也就是计数); v(X=x)表示训练数据集中输入x出现的频率(也就是计数),N是训练数据集的大小。 举个例子,在英汉翻译中,take有多种解释,例如下文中存在7种。 (t1) 抓住: The mother takes her child by the hand. 母亲抓住孩子的手。 (t2) 拿走: Take the book home. 把书拿回家。 (t3) 乘坐: I take a bus to work. 我乘坐公交车上班。 (t4) 量: Take your temperature. 量一量你的体温。 (t5) 装: The suitcase wouldnt take another thing. 这个手提箱不能装下别的东西了。 (t6) 花费: I takes a lot of money to buy a house. 买一座房子要花很多钱。 (t7) 理解: How do you take this passage? 你怎么理解这段话? 在没有任何限制的条件下,最大熵原理认为翻译成任何一种解释都是等概率的,即 P(t1|x)=P(t2|x)=…=P(t7|x)=17 实际中总有或多的限制条件,例如t1、t2比较常见,假设满足 P(t1|x)+P(t2|x)=25 同样根据最大熵原理,可以得出 P(t1|x)=P(t2|x)=15 P(t3|x)=P(t4|x)=P(t5|x)=P(t6|x)=P(t7|x)=325 通常可以用特征函数f(x,y)来描述输入x和输出y之间的某一个事实。一般来说,特征函数可以是任意实值函数,下面我们采用一种最简单的二值函数来定义我们的特征函数 f(x,y)=1,x,y满足某一事实 0,其他 它表示当x和y满足某一事实时,函数取值为1,否则取值为0。 在实际的统计模型中,我们通过引入特征(以特征函数的形式)提高准确率。例如take翻译为乘坐的概率小,但是当后面跟着交通工具的名词bus,概率就变得非常大。于是有 f(x,y)=1,y=乘坐,并且next(x)=bus 0,其他 特征函数f(x,y)关于经验分布P~(X,Y)的期望值,用EP~(f)表示: EP~(f)=∑x,yP~(x,y)f(x,y) 同理,EP(f)表示f(x,y)在模型上关于实际联合分布P(X,Y)的数学期望,类似地有 EP(f)=∑x,yP(x,y)f(x,y) 注意,P(x,y)是未知的,而建模的目标是生成P(y|x),因此我们希望将P(x,y)表示成P(y|x)的函数。因此,利用贝叶斯公式,有P(x,y)=P(x)P(y|x),但P(x)仍然是未知的。此时,只得利用P~(x)来近似。于是便可以将EP(f)重写成 EP(f)=∑x,yP~(x)P(y|x)f(x,y) 以上公式中的求和号∑x,y均是对∑(x,y)∈τ的简写,下同。 对于概率分布P(y|x),我们希望特征函数f的期望应该与从训练数据集中得到的特征期望值相一致,因此提出约束: EP(f)=EP~(f) 即 ∑x,yP~(x)P(y|x)f(x,y)=∑x,yP~(x,y)f(x,y) 把上式作为模型学习的约束条件。假如有J个特征函数fk(x,y),i=1,2,…,J,那么相应就有J个约束条件。 5.5.3模型推导 给定训练数据集,我们的目标是: 利用最大熵原理选择一个最好的分类模型,即对于任意给定的输入x∈Input,可以以概率P(y|x)输出y∈Output。要利用最大熵原理,而目标是获取一个条件分布,因此要采用相应的条件熵H(Y|X),或者记作H(P)。 H(Y|X)=-∑x,yP(x,y)logP(y|x) =-∑x,yP~(x)P(y|x)logP(y|x) 至此,就可以给出最大熵模型的完整描述了。对于给定的训练数据集以及特征函数fk(x,y),k=1,2,…,J,最大熵模型就是求解 maxP∈CH(P)=-∑x,yP~(x)P(y|x)logP(y|x) s.t.EP(fk)=EP~(fk),k=1,2,…,J ∑yP(y|x)=1 或者按照最优化问题的习惯,可以将上述求最大值的问题等价地转化为下面这个求最小值的问题: minP∈C-H(P)=∑x,yP~(x)P(y|x)logP(y|x) s.t.EP(fk)-EP~(fk)=0,k=1,2,…,J ∑yP(y|x)=1 其中的约束条件∑yP(y|x)=1是为了保证P(y|x)是一个合法的条件概率分布。注意,上面的式子其实还隐含一个不等式约束,即P(y|x)≥0,尽管无须显式地讨论它。现在便得到了一个带等式约束的最优化问题,求解这个带约束的最优化问题,所得之解即为最大熵模型学习的解。 这里需要使用拉格朗日乘数法,并将带约束的最优化之原始问题转换为无约束的最优化之对偶问题,并通过求解对偶问题来求解原始问题。首先,引入拉格朗日乘子w0,w1,…,wJ,并定义拉格朗日函数 L(P,w)=-H(P)+w01-∑yP(y|x)+∑Jk=1wk[EP(fk)-EP~(fk)] =∑x,yP~(x)P(y|x)logP(y|x)+w01-∑yP(y|x)+ ∑Jk=1wk∑x,yP~(x,y)fk(x,y)-∑x,yP~(x)P(y|x)fk(x,y) 为了找到该最优化问题的解,我们将诉诸KuhnTucker定理,该定理表明可以先将P看成是待求解的值,然后求解L(P,w)得到一个以w为参数形式的P*,然后把P*代回L(P,w)中,这时再求解w*。本书后面介绍SVM时,还会再遇到这个问题。 最优化的原始问题是 minP∈CmaxwL(P,w) 通过交换极大和极小的位置,可以得到如下这个对偶问题: maxwminP∈CL(P,w) 由于拉格朗日函数L(P,w)是关于P的凸函数,原始问题与对偶问题的解是等价的。这样便可以通过求解对偶问题来求解原始问题。 对偶问题内层的极小问题minP∈CL(P,w)是关于参数w的函数,将其记为 Ψ(w)= minP∈CL(P,w)=L(Pw,w) 同时将其解记为 Pw= argminP∈CL(P,w)=Pw(y|x) 接下来,根据费马定理,求L(P,w)对P(y|x)的偏导数 L(P,w)P(y|x)=∑x,yP~(x)[logP(y|x)+1]-∑yw0-∑x,yP~(x)∑Jk=1wkfk(x,y) =∑x,yP~(x)[logP(y|x)+1]-∑xP~(x)∑yw0-∑x,yP~(x)∑Jk=1wkfk(x,y) =∑x,yP~(x)logP(y|x)+1-w0-∑Jk=1wkfk(x,y) 注意上述推导中运用了下面这个事实: ∑xP~(x)=1 进一步地,令 L(P,w)P(y|x)=0 又因为P~(x)>0,于是有 logP(y|x)+1-w0-∑Jk=1wkfk(x,y)=0 进而有 P(y|x)=expw0-1+∑Jk=1wkfk(x,y)=exp(w0-1)·exp∑Jk=1wkfk(x,y) 又因为 ∑yP(y|x)=1 所以可得 ∑yP(y|x)=exp[w0-1]·∑yexp∑Jk=1wkfk(x,y)=1 即 exp(w0-1)=1∑yexp∑Jk=1wkfk(x,y) 将上面的式子代回前面P(y|x)的表达式,则得到 Pw=1Zw(x)exp∑Jk=1wkfk(x,y) 其中, Zw(x)=∑yexp∑Jk=1wkfk(x,y) Zw(x)称为规范化因子; fk(x,y)是特征函数; wk是特征的权值。由上述两式所表示的模型Pw=Pw(y|x)就是最大熵模型。这里,w是最大熵模型中的参数向量。注意,我们之前曾经提过,特征函数可以是任意实值函数,如果fk(x,y)=xk,那么这其实也就是5.5.2节中所说的多元逻辑回归模型,即 Pj=P(y=j|x)=ewj·x∑Jk=1ewk·x 5.5.4极大似然估计 下面,需要求解对偶问题中外部的极大化问题 maxwΨ(w) 将其解记为w*,即 w*= argmaxwΨ(w) 这就是说,可以应用最优化算法求对偶函数Ψ(w)的极大化,得到w*,用来表示P*∈C。这里,P*=Pw*=Pw*(y|x)是学习到的最优模型(最大熵模型)。于是,最大熵模型的学习算法现在就归结为对偶函数Ψ(w)的极大化问题上来。 前面已经给出了Ψ(w)的表达式: Ψ(w)=∑x,yP~(x)P(y|x)logP(y|x)+w01-∑yP(y|x)+ ∑Jk=1wk∑x,yP~(x,y)fk(x,y)-∑x,yP~(x)P(y|x)fk(x,y) 由于,其中 ∑yPw(y|x)=1 于是将Pw(y|x)代入Ψ(w),可得 Ψ(w)=∑x,yP~(x)Pw(y|x)logPw(y|x)+ ∑Jk=1wk∑x,yP~(x,y)fk(x,y)-∑x,yP~(x)Pw(y|x)fk(x,y) =∑Jk=1wk∑x,yP~(x,y)fk(x,y)+∑x,yP~(x)Pw(y|x)logPw(y|x)- ∑Jk=1wk∑x,yP~(x)Pw(y|x)fk(x,y) =∑x,yP~(x,y)∑Jk=1wkfk(x,y)+∑x,yP~(x)Pw(y|x)logPw(y|x)-∑Jk=1wkfk(x,y) =∑x,yP~(x,y)∑Jk=1wkfk(x,y)-∑x,yP~(x)Pw(y|x)logZw(x) =∑x,yP~(x,y)∑Jk=1wkfk(x,y)-∑xP~(x)logZw(x)∑yPw(y|x) =∑x,yP~(x,y)∑Jk=1wkfk(x,y)-∑ xP~(x)logZw(x) 注意其中倒数第4行至倒数第3行运用了下面这个推导: Pw(y|x)=1Zw(x)exp∑Jk=1wkfk(x,y)logPw(y|x)=∑Jk=1wkfk(x,y)-logZw(x) 下面来证明对偶函数的极大化等价于最大熵模型的极大似然估计。已知训练数据的经验概率分布P~(X,Y),条件概率分布P(Y|X)的对数似然函数表示为 LP~(Pw)=log∏x,yP(y|x)P~(x,y)=∑x,yP~(x,y)logP(y|x) 当条件概率分布P(y|x)是最大熵模型时,对数似然函数为 LP~(Pw)=∑x,yP~(x,y)logPw(y|x) =∑x,yP~(x,y)∑Jk=1wkfk(x,y)-logZw(x) =∑x,yP~(x,y)∑Jk=1wkfk(x,y)-∑x,yP~(x,y)logZw(x) 对比之后,不难发现 Ψ(w)=LP~(Pw) 既然对偶函数Ψ(w)等价于对数似然函数,也就证明了最大熵模型学习中的对偶函数极大化等价于最大熵模型的极大似然估计这一事实。由此,最大熵模型的学习问题就转换为具体求解“对数似然函数极大化的问题”或者“对偶函数极大化的问题”。 5.6应用实例: 多分类问题 本节将演示使用多元逻辑回归模型对手写数字进行识别。作为例子,这里要完成的任务是对0~9这10个手写数字进行分类。示例中使用到的数据集是scikitlearn内置的数据集。此数据集是创建于1998年的UCI手写数字数据集中测试集的复制,包含1797个样本,每个样本是一张8×8的手写数字图片。因此每个样本也可以看作有8×8 = 64维特征,每维特征是个0~16的整数。 5.6.1数据初探 首先加载数据,并可视化数据集中前50个样本,结果见图55。 from sklearn.datasets import load_digits # 加载内置数据集 digits = load_digits() # 可视化数据集的前几项 _, axes = plt.subplots(5, 10) axes = axes.reshape(1,-1) for ax, image in zip(axes[0], digits.images[:50]): ax.set_axis_off() ax.imshow(image, cmap=plt.cm.gray_r) plt.show() 图55手写数字数据集的样例 上面的代码中,首先使用load_digits()函数加载了scikitlearn内置的手写数字数据集。然后,使用matplotlib包中pyplot的subplots()方法生成一个5行10列的布局,返回的是每张子图的坐标。最后,将每个样本填入对应的坐标中,并使用imshow()方法将数据绘制为一张二维图片,便得到了图55的效果。显然,数据集的确是0~9等数字的手写图片,每个数字都有多张图片作为样本,数字本身则是标注。 5.6.2建模 显然这是个多分类问题,数据集一共有10个类别。仍然使用scikitlearn下linear_model模块中提供的LogisticRegression类建模,不过在参数的设置上需要格外注意。 from sklearn.linear_model import LogisticRegression from sklearn.model_selection import train_test_split from sklearn.preprocessing import StandardScaler X_raw = digits.images Y = digits.target # 将原始三维矩阵转化成二维矩阵 X = X_raw.reshape(X_raw.shape[0], -1) # 将原始数据划分成训练和测试集两个部分 X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size=0.2) # 归一化原始数据 scaler = StandardScaler() X_train = scaler.fit_transform(X_train) X_test = scaler.transform(X_test) # 训练模型 clf = LogisticRegression(penalty='l1',fit_intercept=True,solver='saga',tol=0.001,max_iter=1000) clf.fit(X_train, Y_train) 上述代码中,首先将原本三维的数据转化成了二维,方便后续处理。注意,这个操作和机器学习中的降维算法是截然不同的。此处只是将原本1797×8×8的三维数组变形为1797×64的格式,原本每个样本是64个特征,变化后还是64个特征,并不改变原始的特征取值。而后续章节将介绍的降维算法,属于机器学习中的流型学习(manifold learning),一般会改变数据的原始取值。 然后采用train_test_split()函数将数据分成训练集和测试集。函数的第一个参数是特征数据X,第二个参数是标注数据Y,第三个参数设定测试集的占比。此处划分比例被设置为0.2,即原始数据的80%将用于模型训练,而剩余20%的数据则是测试数据集。紧接着,使用StandardScaler类的fit_transform()方法,对训练集进行zscore标准化处理。fit_transform()方法等价于先调用fit()方法得到标准化模型scaler,再用这个实例的transform()方法对训练集进行标准化。然后,调用基于训练集得到的预处理模型scaler的transform()方法,对测试集的特征进行同样的标准化处理。标准化的目的是防止特征各维度值域的差异对模型造成影响,注意训练集和测试集的标准化计算必须一致。 最后,声明一个LogisticRegression类的实例,设置必要参数后,调用fit()方法基于标准化后的训练集和标注训练模型。这次建模采用的优化算法是saga,正则项为L1; 此外,模型的最大迭代次数提升到了1000次,而模型的停止迭代误差被提升到了0.001。为什么要这样调整模型的参数设置呢?这就不得不深入到优化算法各自的特点了。 之前二分类的示例中,采用的优化算法是在小数据集上效果不错的liblinear算法。虽然liblinear算法也可以处理多分类任务,但是其实现的坐标下降法并不能学习到一个真正的多分类模型。liblinear算法实际上使用的是“一对多”(onevsrest)框架处理多分类问题,即将每个类别和剩余所有类别看成一个二分类问题,然后训练多个二分类模型集成后处理多分类问题。第7章将详细介绍这种框架。 剩余的newtoncg、lbfgs、sag和saga均能学习到真正的多分类模型。我们知道,L1正则具有将参数稀疏化的能力,从而简化学习到的模型参数。观察图53的任一个数字样本,能发现虽然手写数字共64个特征,而实际上包含了有效信息的特征并不多。有效信息是图中黑色笔画的部分,白色部分并没有有效信息。因此,这个问题中采用L1正则是更好的选择,又因为剩余的算法中只有saga算法能处理L1正则,所以优化算法采用saga。剩余的max_iter和tol参数都是为了保证模型更快的收敛而设置的,感兴趣的读者可自行尝试。 最后,简单使用score()方法观察学习到的模型的平均精确度。 score = clf.score(X_test, Y_test) print("Test score with L1 penalty: %.4f" % score) # 以下为打印内容 Test score with L1 penalty: 0.9556 每个类别的详细分类效果,可使用metrics模块下的classification_report()方法,只需要提供预测结果和真实标注即可。 from sklearn import metrics predicted = clf.predict(X_test) print("Classification report for classifier %s:\n%s\n" % (clf, metrics.classification_report(Y_test, predicted))) # 打印结果如下: precisionrecallf1-scoresupport 01.001.001.0029 1 0.86 0.92 0.89 39 2 0.98 0.98 0.98 44 3 1.00 0.93 0.96 40 4 1.00 0.97 0.98 30 50.980.980.9842 6 1.00 0.97 0.99 36 7 0.97 1.00 0.98 30 8 0.94 0.89 0.92 38 9 0.86 0.94 0.90 32 accuracy 0.96 360 macro avg 0.96 0.96 0.96 360 weighted avg 0.96 0.96 0.96 360 上述打印内容中,第一列为类别名称,后面3列为precision、recall、f1值等指标,最后一列是测试集中该类别的数量。可以发现模型对数字0的识别效果最好,对数字1的识别效果最差。此外,metrics模块还提供了一个函数plot_confusion_matrix()能绘制模型的混淆矩阵,只需提供分类器、测试集和标注。 disp = metrics.plot_confusion_matrix(clf, X_test, Y_test) disp.figure_.suptitle("Confusion Matrix") #设置图片名称 plt.show() 运行上述代码,结果见图56。其中纵坐标为真实标注,横坐标为预测类别。显然,对角线上的数字代表分类正确的样本数,非对角线上的数字代表了误分类样本数。例如,真实标注为1的第二行,可以发现有两个样本被误识别为了数字“8”,还有一个样本被误识别为了数字“9”。此外,图56还使用不同的颜色代表样本的数量,能够非常直观地发现各类别误分类的具体情况。 图56模型的混淆矩阵