造出一艘logistic模型 | 【logistic从生产到使用】（下） | 数说 · 算法

前几天飞扬博士更新了一篇算法文章，关于softmax regression的，它是logistic模型的扩展，因此要是能有些logistic regression的底子就看起来非常容易，因此在发softmax regression之前，重新复习一下logistic模型。

一句话介绍：

logistic regression，它用回归模型的形式来预测某种事物的可能性，并且使用优势（Odds）来考察“某事物发生的可能性大小”。

上篇介绍了logistic模型的原理，如果你只是想使用它，而不需要知道它的生产过程，即拟合方法及编程实现，那么上篇就足够了。如果你想知道它的上游生产，那么请继续。

本篇着重剖析logistic模型的内部生产流程、以及每一个流程的工作原理，暴力拆解。

上下两篇的大纲如下：

【上篇：使用篇】

1. Logistic回归模型的基本形式

2. logistic回归的意义

（1）优势

（2）优势比

（3）预测意义

3. 多分类变量的logistic回归

（1）无序多分类logistic回归

（2）有序多分类：比例优势模型

（3）有序多分类：偏比例优势模型

4．附：Logistic回归模型建模指南

【下篇：生产篇】

1. 模型的拟合

（1）回归模型的拟合流程

（2）logistic回归的拟合

2. 代码实现

（1）Python

（2）SAS

（3）Matlab

在微信后台回复【logistic】查看上下两篇

logistic回归：从生产到使用【下：生产篇】

1.模型的拟合

（1）回归模型的拟合流程

很多统计出身、尤其是经济统计出身的朋友，并不知道回归模型拟合的标准流程，只知道线性回归用最小二乘法。其实最小二乘问题、最小二乘法、极大似函数等，以及其他回归中用到的梯度下降算法、牛顿法等等，都是不同的东西，首先来看一下回归的一个标准拟合流程（点击查看大图）：

下面结合大家都用过的的一元线性模型为例，走一遍这个流程：

• 选择模型

首先我们要对具体分析的问题选择一个模型，对于连续变量我们用线性模型，对于定性变量我们可以用logistic模型（什么情况下用logistic模型，它可以解决什么问题？在我们的上篇重点介绍了），不同的模型，它的目标函数，以及后面选择的优化算法都会不同。这里我们以最常见的一元线性模型为例：

Y=a+bX

• 建立目标函数：

对于上面的例子，我们希望最终能得到这样的回归系数

：

使得由系数计算出来的

值与真实的y值之间差距绝对值最小，即

但其实，绝对值不便于后面的计算。在下一步“选择最优算法”的时候，要务实地考虑这个目标如何实现，绝对值不容易计算，我们用平方来代替：

这里，是的Q这个函数的值最小，就是我们的目标。Q就是目标函数。我们把目标变成一个求最小的问题，这个问题就是最小二乘问题。

对于logistic模型，我们的目标函数就不是最小二乘了，而是极大似然，其实它们之间不是对立的，最小二乘可以通过极大似然推导出来。这在后面会说。

• 选择最优算法

为了使得实现目标函数，即误差的平方最小，我们需要选择一个算法来实现。根据微积分，我们只需要把Q对a和b分别进行求导，另其导数为0，得出来的函数就是最小值（Q函数是二次函数，又是非负的）。

因此这里我们就选择求导为0的方法，也就是一般来说的最小二乘法。

并不是所有的情况都可以通过求导来获得，例如当目标函数求导为0之后，求不出解析解，那么只能找其他的方法了，一般来说可以通过迭代的方法。这个后面也会说。

• 对目标函数进行优化

这里的“优化”当然就是“求最小”，我们使用求导为0的方法。

• 拟合出最优的回归系数

求解上一步中的两个导数为零的函数，最终解得：

（2）logistic回归的拟合

弄清楚了回归模型的拟合流程，现在我们看一下logistic模型是如何“生产”出来的。与线性模型相比，logistic很多方法不一样。

• 选择模型——logistic模型

不多说，我们现在需要构建一个因变量y（0=否/1=是），m个自变量的logistic模型（点击查看大图）：

如果对这些基本形式不是很熟悉，可以查看上篇，有详细的解说。

• 建立目标函数——极大似然

Logistic的目标函数是极大似然函数，这是本【生产篇】的一个重头戏，我们要介绍极大似然的思想、logistic模型如何运用极大似然思想、以及极大似然与最小二乘之间的关系。

极大似然的思想如下：

假设我们有样本（X1，Y1），那么

当Y1=1时，这个样本发生的概率为

P{Y=1}，

当Y1=0时，这个样本发生的概率为

P{Y=0}=1-P{Y=1}

它可以直接用一个式子来表达：

各位可以想一下是不是，当Y1=1时，等价于第一个式子，Y1=0时，等价于第二个式子。

这时有朋友会问了：

“（X1，Y1）不是样本吗，既然是样本，那就说明它是一个已经发生过了的历史数据，它发生的概率应该是100%呀！”

是的，因为

这个式子是logistic回归模型得出的、样本（X1，Y1）发生的概率，而实际上这个样本是已经发生过了的“历史事实”，正因此，我们要使这个式子得到的值尽可能的大到100%，以使得模型的情况能最贴近现实，所以，

就是最大似然函数，因为我们不止有一个样本，且样本之间是相互独立的，所以n个样本的发生概率就是这n个样本概率的乘积：

至此，我们建立了logistic的目标函数。

还没完，这里还有人问，

“为什么logistic的目标函数不能是最小二乘？而是最大似然？”

线性回归中，因变量Y是连续的，因此我们用拟合出来的

与真实之间的Y的差别平方作为目标函数，目标是使误差平方最小。而logistic模型，因变量Y是分类函数，比如0、1模型中我们计算的缺是Y的发生概率P{Y=0}、P{Y=1}。因此适合用最大似然。

实际上，最小二乘和极大似然并不对立。最小二乘是可以用极大似然推导出来的。

下面给出推导过程，不敢兴趣的可以直接跳过，知道两者相关就好了：

现在有回归模型，模型希望通过参数θ和若干自变量X拟合出因变量

，与真实之间的Y之间有误差：

e是误差项，服从正态分布（回归模型的经典假设）：

因此有：

根据极大似然思想，我们希望拟合出来的参数θ，可以使Y这个“历史事实”发生的概率最高，这个目标函数为：

为了简化计算，我们取对数：

之所以简化，因为对于这两个单调函数来说，l(θ)=lnL(θ)达到最大时，L(θ)也达到最大。而让l(θ)最大，也就是让

取最小，这样就转化成了最小二乘的问题，回顾一下看是不是？

• 选择最优算法——梯度下降：

对于logistic的目标函数，我们可不可以再用求导为0的算法？答案是不行，不信？不信走两步——

首先需要取对数化简L(θ)（点击查看大图）：

这里，

，下表i表示第i条数据。

用L(θ)对θ求导：

这时我们发现，求导之后令其等于0，是找不到关于θ的解析解的，你可以试一下。因此求导为0的这个算法失效。只能用计算机迭代搜索来找到最优的解了。

梯度下降（上升）算法，就是从一个初始状态出发，一点一点的搜索到全局的最小（最大）值，我们以梯度下降，搜索最小值为例：

这个思想可不可以用一个数学表达式表示呢？可以，假设我们的目标函数为J(θ)

先看一下迭代的思想，再具体说一下每个项都代表什么：

迭代就是这么进行的，这里θ会不断进行更新，直到达到局部最小值点。那么后面的更新项

是怎么来的呢？它什么可以使参数达到局部最小值？

Answer1：α是学习率，代表了每一次迭代的更新程度。α过大过小都不好，过小则参数收敛速度很慢，要很久才能达到最小值点；过大则很容易在极小值点出徘徊，步子太大了，老是走不到关键点上。

Answer2：

这个设计很妙，它起什么作用呢？

想象一下，当J(θ) 接近最小值的时候，它的曲线会变的非常平滑——斜率，也即J(θ) 导数的值，将会接近于0，那么更新将趋近于停滞。当达到最小值的时候，θ就会停止更新了。

现在梯度下降算法基本搞明白了，但是，这里我们是要最大化似然函数啊，应该求的是最大值啊。不错，logistic模型中我们应该使用梯度上升算法，和梯度下降算法的原理是一样的，比如，求J(θ) 的最大值，其实也就是求-J(θ) 的最小值，加个负号，就可以用梯度下降算法了。

因此这里，我们的目标函数实际是：

注意加了负号，因此我们可以使用梯度下降算法，迭代的数学表达式为：

梯度下降法在具体实践上，分为“批量梯度下降”和“随机梯度下降”：

批量梯度下降，是进行迭代时，使用所有的样本，正如上面的式子中有一个sigma求和函数，考虑的是所有样本。这个方法处理100个样本的数据还可以，如果有数十亿的样本数据，这个方法复杂度就太高了。

随机梯度下降，是迭代时，只用一个样本来更新，使用这个样本迭代完之后，再用下一个样本再更新迭代。那么上面的式子中就没有sigma求和函数了，不明白不要紧，下面的python代码中，会有一个图，来解释两者的区别。

2.代码实现

下面我们以一个名字为testSet.txt的数据为例，它的样子如下图：

第一行、第二行是自变量X1，X2，第三行作为因变量Y。该数据来源于《Machine Learining》这本书中。

（1）Python

首先看一下Python如何实现梯度下降的一轮迭代的：

对于批量梯度下降：

对于随机梯度下降：

以上就是批量梯度下降和随机梯度下降中，每一轮迭代的思想，以及Python实现。下面要写出具体的代码：

（2）SAS

直接把testSet.txt文件导入SAS，自变量命名Y，因变量命名X1，X2，X3...

proc logistic desc data=a; model Y= X1 X2; run;

注意SAS默认把数值较高的作为对照值，即Y=0、1时，SAS默认1值作为不发生的情况，0值作为发生的情况，而事实相反，因此应该把顺序换一下，使用descending关键句。

（3）Matlab

不多说，Python代码出来了，在Matlab中稍修改一下就可以，代码如下图。

只是，数说君发现Matlab和Python的计算结果差的蛮大的。排查了一下觉得应该不是程序问题，因为迭代的时，初始值一样，第一轮迭代计算没有差别，第二轮迭代计算，出现了小的差别，第三轮第四轮....后面的差别越来越大。有找到原因的同学可以加数说君微信（在微信公众号“数说工作室”）私聊。