[机器学习]逻辑回归公式推导及其梯度下降法的Python实现

一般来说，二项逻辑斯谛回归模型是一个二分类判别模型，由条件概率分布$P(Y|X)$表示，随机变量$X$为实数，$Y$取值0或者1。我们通过比较$P(Y=1|x)$和$P(Y=0|x)$值大小来判断给定x的类别为1还是0。

从线性模型推导

我们先说广义的线性回归：$y=wx+b$，这里 $y$为回归的目标，x为输入的数据，$w$和$b$分别为需要学习的参数和偏置，我们简化为$y=wx$ ，此时$w$为参数向量，包括了偏置$b$。如果此时我们把线性回归的目标$y$设为一个事件的对数几率(log odds)，即事件发生的概率$p$与不发生的概率$1-p$的比值再取对数，也就是：

$$\log(\frac{p}{1-p})=wx$$
那么这就是逻辑斯谛回归模型，此时 $p$可以为类别1事件发生的概率$P(Y=1|x)$， 则 $1-p$为类别0事件发生的概率 $P(Y=0|x)$，我们用 $p$表示上述式子：
$$p=\frac{exp(wx)}{1+exp(wx)}= \sigma(wx)$$
就得到了我们平时接触最多的logistic函数（ $\sigma$）的形式，那么换一种角度理解，这个模型就是一个线性模型 $wx$经过一个非线性变换 $\sigma$得到的，更通俗的说，是一个神经网络，这个网络除了输入层只有一个输出神经元，并且使用了sigmoid激活函数 $\sigma$，最后输出的p表示类别为1的概率。

极大似然法估计与交叉墒

学习逻辑斯谛回归模型时，对于给定的数据集$T=\{(x_1,y_1),(x_2,y_2),\dots,(x_N,y_N)\}$，我们通过极大似然估计法来估计出一套参数$\hat{w}$，使模型对于此数据集中样本值的发生概率最大，也就是使模型能尽可能拟合现有数据集中现有的样本。

由于数据集中每个样本$(x_i,y_i)$之间相互独立，用$p(x_i)$表示模型判别$x_i$为类别1的概率，则$1-p(x_i)$表示模型判别$x_i$为类别0的概率，那么整个数据集$T$所有样本发生的概率为每个样本对应实际类别概率的连乘：

$$\prod _1 ^N [p(x_i)]^{y_i}[1-p(x_i)]^{1-y_i}$$
这就是我们需要最大化的似然函数，由于连乘的性质，我们一般采用对数似然函数将乘法变成加法，于是有：

$$L(w)=\sum _1 ^N [y_i \log p(x_i) + (1-y_i) \log (1-p(x_i))]$$
对 $L(w)$求极大值，得到对w的估计值 $\hat w$。

从另外一个角度来看，带负号的$L(w)$其实非常像交叉墒：交叉墒在一定程度上是衡量着两个分布之间的差异，这里具体是指$p(x_i)$和$y_i$。当模型学到的$p(x_i)$和真实目标分布$y_i$越相近，模型的效果越好，交叉墒越小。理想情况下，当两个分布完全一致时，此时相对墒为0，交叉墒达到最低值，即真实目标分布的墒 $y_i \log (y_i)$，详情见如何通俗的解释交叉熵与相对熵?。当带负号的$L(w)$达到极小值时，也就是$L(w)$达到极大值，此时的模型参数是理想的，也是我们最终需要达到的目标，这与我们前面用似然函数估计分析的结果一致，可以看出这两者实际上是相通的。

梯度下降法最优化

我们把求$L(w)$极大值问题转化为求$-L(w)$的极小值的最优化问题。最优化的方法通常是梯度下降法和牛顿法，我们具体给出梯度下降法的做法。我们把概率$p(x_i)$表示为非线性函数$\sigma(wx_i)$，因此目标函数重新写成如下：

$$L= \sum_i^N -y_i \log \sigma(wx_i)-(1-y_i) \log( 1-\sigma(wx_i))$$
对 $wx_i$进行求导，有：

$$\begin{align} \frac{\partial L}{\partial (wx_i)} & = \frac{\partial L}{\partial (\sigma (wx_i))} \frac{\partial (\sigma (wx_i))}{\partial (wx_i)} \\ & = [-y_i \frac {1}{\sigma (wx_i)}- (1-y_i) \frac {-1} {1-\sigma(wx_i)}] \sigma(wx_i)(1-\sigma(wx_i)) \\ & = -y_i(1-\sigma(wx_i))+(1-y_i)\sigma(wx_i) \\ & = -y_i+ y_i \sigma(wx_i)+ \sigma(wx_i)-y_i \sigma(wx_i) \\ & = \sigma(wx_i)-y_i \end{align}$$
其中利用了sigmoid函数的导数 $\sigma'=\sigma(1-\sigma)$。从直观上解释：对于单个样本 $x_i$, 线性函数 $wx_i$的梯度为其非线性激活后的值减去真实标签值，也就是输出值 $p(x_i)$减去 $y_i$，更具体来说，当 $y_i$为1时，梯度为 $p(x_i)-1$，当 $y_i$为0时，梯度为 $p(x_i)$。

如果用批量梯度下降法，那么每个epoch，参数$w$的梯度为：

$$\begin{align} \frac{\partial L}{\partial (w))} & = \sum_i^N \frac{\partial L} {\partial (wx_i)} \frac{\partial (wx_i)} {\partial w} \\ & = \sum_i^N (\sigma(wx_i)-y_i) \times x_i \end{align}$$

Python实现（批量梯度下降法）

我们使用iris数据集中的前100行数据，实现一个简单的逻辑回归二分类器，使用批量梯度下降法：

from sklearn import datasets
import pandas as pd
import numpy as np
from matplotlib import pyplot as plt
from numpy.linalg import inv

iris = datasets.load_iris()
X = iris.data[:100, :]
y = iris.target[:100].reshape((100, -1))


def logit(x):
    return 1. / (1 + np.exp(-x))


m, n = X.shape
alpha = 0.0065 // 步长
w = np.random.random((n, 1)) // 参数矩阵
maxCycles = 30
J = pd.Series(np.arange(maxCycles, dtype=float))

for i in range(maxCycles):
    h = logit(np.dot(X, w)) // 输出估计值h
    J[i] = -(1 / 100.) * np.sum(y * np.log(h) + (1 - y) * np.log(1 - h)) // 记录目标函数的值
    error = h - y //计算wx的梯度，输出值减去真实值
    grad = np.dot(X.T, error) //计算w的梯度
    w -= alpha * grad // 更新参数w，使用负梯度最小化J
print w
J.plot()
plt.show()

牛顿法与梯度下降法比较

牛顿法是直接利用函数有极小点的必要条件，即一阶导数 f’导数为零，所以每次迭代后的点 $x_{k+1}$ 的特点就是 $f'(x_{k+1})=0$ ，又由于一阶导数可以通过二阶导数与步长的乘积来计算（二阶泰勒展开）：$f'(x_{k+1})=f'(x_k)+f''(x_k)(x_{k+1}-x_k)$ ，因此就有 $f'(x_k)+f''(x_k)(x_{k+1}-x_k)=0$ 即每次更新的步长等于 $\Delta x = x_{k+1}-x_{k}=-f'(x_k)/f''(x_k)$ ，这个 $\Delta x$ 在梯度下降法中的体现就是 $\alpha$ ，只不过梯度下降法每次更新都是以一个非常小的固定步长 $\alpha$ （一般情况下）在更新，而牛顿法一步到位直接找到一个相对最优的点，利用了二阶导数的信息，因此迭代次数更少，收敛快。

举例子，比如一个简单的强凸函数 $f=x^2$ , 无论起始点是哪，牛顿法总能一次迭代就找到最小值点0, 而梯度下降法只能根据步长和起始点得位置慢慢逼近0，但是实际上函数 $f$ 可能要复杂得多，函数可能是整体上非凸，但是局部是凸函数（比如最优点附近），直接用牛顿法迭代多次也不能保证收敛到最优点，因此需要先用梯度下降找到一个相对好的解后再用牛顿法可能效果比较好（根据曾文俊的回答），从这里我们也可以看出牛顿法对函数的性质以及初始点的位置选择比较挑剔，另外一个是二阶导数矩阵的逆矩阵计算量比较大（当 $x_k$ 为参数矩阵时，此时 $f''(x_k)$ 为一个二阶偏导数矩阵，即Hesse 矩阵），通常使用拟牛顿法。

个人思考：前面的 $f=x^2$ 最大阶数只有2，我们利用二阶导数能完美解决它的最优化问题；但是当 $f$ 的阶数 n 大于2时，此时用泰勒展开式展开到2阶（牛顿法只展开到2阶），后面还有一个高阶无穷小 $R_n$ ，因此现在展开的等式两边都是约等于的关系（如果不考虑 $R_n$ ），即可推出 $f'(x_{k+1}) \approx f'(x_k)+f''(x_k)(x_{k+1}-x_k)$ ，严格上来说，此时求得的 $\Delta x$ 并不能使参数更新到最优点，只能说是“相对最优”的点，所以多次迭代还是有必要的。