逻辑回归

1. 历史背景

1830 年代和 1840 年代，在 Adolphe Quetelet 的指导下，Pierre François Verhulst 将逻辑函数（logistic function）作为人口增长模型并命名为“逻辑函数”。Verhulst 在 1830 年代中期首先设计了这个函数，但没有具体说明他如何将曲线拟合到数据中。随后，在 1838 年发表了一篇简短的笔记。最后在 1844 年进行了扩展分析并命名了这个函数（1845 年出版），Verhulst 通过使曲线通过三个观察点来确定模型的三个参数，这产生了较差的预测。逻辑函数最早应用就是作为人口增长的通用模型，满足微分方程式 $\displaystyle {\frac {dP}{dt}}=rP\left(1-{\frac {P}{K}}\right)$，其中，$P$ 代表种群大小，$t$ 代表时间，常数 $r$ 定义为增长率，$K$ 是承载容量，微分方程结果为 $\displaystyle P(t)={\frac {K}{1+\left({\frac {K-P_{0}}{P_{0}}}\right)e^{-rt}}}$，$\displaystyle \frac {K-P_{0}}{P_{0}}$ 是常数可以移动到指数上，相当于对 $t$ 进行平移，就可以得到逻辑函数标准形式。

1883 年，逻辑函数在化学领域中独立开发出来作为自催化模型（Wilhelm Ostwald）。自催化反应是指反应中的一种产物本身就是同一反应的催化剂，而其中一种反应物的供应是固定的。出于与人口增长相同的原因，这自然会产生逻辑方程：反应是自我强化的，但受到限制。

1920 年，逻辑函数被 Raymond Pearl 和 Lowell Reed 再次独立地发现为人口增长的模型，这导致它在现代统计学中的应用。Pearl 和 Reed 将该模型首次应用于美国人口，并通过三个点对曲线进行了初步拟合，与 Verhulst 一样，这再次产生了糟糕的结果。他们最初不知道 Verhulst 的工作，大概是从 L. Gustave du Pasquier 那里了解到的，但他们没有给他多少信任，也没有采用他的术语。在 1925 年，Verhulst 的优先权得到承认，Udny Yule 重新使用了“logistic”一词，并一直沿用至今。

1930 年代，概率模型（probit model）由 Chester Ittner Bliss 开发和系统化，并且 John Gaddum 通过 Ronald A. Fisher 的最大似然估计拟合该模型。probit 模型主要用于生物测定，早在 1860 年的工作就已经开始了。probit 模型影响了 logit 模型的后续发展。

1943 年，逻辑模型（logistic model）可能首先被 Edwin Bidwell Wilson 和他的学生 Jane Worcester 用作生物测定中概率模型的替代方案。然而，logistic 模型作为 probit 模型的一般替代方案的发展主要归功于 Joseph Berkson 数十年来的工作，从 1944年开始，他通过与“probit”类比创造了“logit”。logit 模型最初被认为不如 probit 模型，但“逐渐达到了与 probit 平等的地位”，特别是在 1960 年至 1970 年间。到 1970 年，统计期刊中使用logit 模型与 probit 模型数量相当，然后前者超越了后者。这种相对流行是由于在生物测定之外采用了 logit，而不是在生物测定中取代 probit，以及它在实践中的非正式使用；logit 的受欢迎程度归功于 logit 模型的计算简单性、数学特性和通用性，允许其在各种领域中使用。

1966年，多项式logit模型由 David Cox 引入，大大增加了logit模型的应用范围和普及度。1973 年，Daniel McFadden 将多项式 logit 与离散选择理论，特别是 Luce 的选择公理联系起来，表明多项式 logit 遵循 independence of irrelevant alternatives (IIA) 假设并将备选方案的几率解释为相对偏好，这给出了一个逻辑回归的理论基础。

2. 算法模型

逻辑回归实际上是分类模型，之所以名字中带有回归，我想是因为它是对概率值进行回归的。逻辑回归是对事件发生的概率进行建模，根据分类事件的数量分为二元逻辑回归和多元逻辑回归，

2.1 逻辑函数

阐述逻辑回归模型之前，先介绍一下逻辑函数（logistic function）。

逻辑函数的公式如下

$\displaystyle f(x)={\frac {L}{1+e^{-k(x-x_{0})}}}$

式中，$x_0$ 是曲线的中点坐标，$L$ 是曲线的最大值，$k$ 是逻辑增长率或曲线的陡度。逻辑函数最早用于人口增长预测，$f(x)$ 是人口数量，$x$ 代表时间。

当 $\displaystyle L=1,k=1,x_{0}=0$ 时就得到了 sigmoid 函数

$\displaystyle f(x)={\frac {1}{1+e^{-x}}}$

与正态分布函数进行类比，逻辑函数改为

$\displaystyle f(x)={\frac {1}{1+e^{-(x-\mu )/s}}}$

该函数图形如图所示，是一条 S 形曲线（sigmoid curve），可以发现形状与正态分布函数相似。该函数曲线以点 $\displaystyle(\mu,\frac{1}{2})$ 为中心对称，满足

$f(x)=1-f(-x)$

逻辑函数还有一个非常好的性质，就是其导数易求

$f'(x)=f(x)(1-f(x))$

逻辑函数的导函数满足概率密度公式的要求，因此逻辑函数就是一种概率分布函数，至于为什么在分类问题中比正态分布好用还需考证。

logistic model 与 probit model 的区别就是在于使用逻辑函数取代了正态分布函数。

2.2 二元逻辑回归

假设有数据集

$D=\{(x_1,y_1),(x_2,y_2),\dots,(x_k,y_k),\dots,(x_N,y_N)\}$

二元逻辑回归就是在分类数是 2 的情况下，也就是二项分布的情况下的逻辑回归。回忆一下伯努利分布，其中的 $y=1$ 事件对应的概率分布就是由逻辑函数 $f(x)$ 给出，由于与 $x$ 有关，因此很明显是条件概率，如下所示

$p(y=1|\, x)={\frac {1}{1+e^{-(x-\mu )/s}}}$

由于是二项分布，因此

$p(y=0|\, x)=1 - p(y=1|\, x)$

我们用 $p(x)$ 对 $p(y=1|\, x)$ 进行简写，同时把上式标准化

$\displaystyle p(x)={\frac {1}{1+e^{-(\beta _{0}+\beta _{1}x)}}}$

其中 $\beta_0=-\mu/s$ 称为截距，$\beta_1=1/s$。那么对特征 $x_k$ 的类别 $y_k$ 概率拟合就是

$p_k=p(x_k)$

有了概率分布之后，我们只需要用最大似然估计确定参数的值即可。如前所说，二元逻辑回归就是在伯努利分布的情况下，那么伯努利分布的似然函数就是

$\begin{aligned} \displaystyle L&= \prod _{k:y_{k}=1}p_{k}\,\prod _{k:y_{k}=0}(1-p_{k})\\ &= \prod^N_{k=1}p_{k}^{y_k}(1-p_{k})^{1-y_k}\\ \end{aligned}$

对 $L$ 取对数，得到对数似然函数

$\begin{aligned} \displaystyle \ell&=\sum _{k:y_{k}=1}\ln(p_{k})+\sum _{k:y_{k}=0}\ln(1-p_{k})\\ &= \sum _{k=1}^{N}\left(\,y_{k}\ln(p_{k})+(1-y_{k})\ln(1-p_{k})\right)\\ \end{aligned}$

上式是有关 $\beta_0,\beta_1$ 的函数，为了使对数似然函数最大，分别对 $\beta_0,\beta_1$ 求偏导数

$\displaystyle {\frac {\partial \ell }{\partial \beta _{0}}}=\sum _{k=1}^{K}(y_{k}-p_{k})$ $\displaystyle {\frac {\partial \ell }{\partial \beta _{1}}}=\sum _{k=1}^{K}(y_{k}-p_{k})x_{k}$

结果和线性回归太相似了，$p_k$ 如果是线性函数，那么上式就是线性回归。事实上，$p(x)$ 的作用就是将线性函数 $\beta _{0}+\beta _{1}x$ 的结果映射到 $[0,1]$ 概率区间。如果让上式分别等于 0，则得到似然方程。

不过，尽管结果相似，但上式是非线性的，因此无法像线性回归一样，求出似然方程组的解析解。因此，需要用数值方法，如牛顿迭代法和梯度下降法等等。牛顿迭代法是求解的似然方程组，梯度下降法则是直接求似然函数的极值，区别在于牛顿迭代法需要二次导数。下面给出梯度下降法的计算公式，其实很简单，沿与梯度相反的方向移动即可，导数已经给出，梯度则为

$\bigtriangledown p(\beta_0,\beta_1)=(\frac {\partial \ell }{\partial \beta _{0}},\frac {\partial \ell }{\partial \beta _{1}})$

假设第 $t$ 次的迭代结果为 $\beta^t_0,\beta^t_1$，则

$\beta^{t+1}_0=\beta^t_0-\alpha\sum _{k=1}^{K}(y_{k}-{\frac {1}{1+e^{-(\beta^t_0+\beta^t_1x_k)}}})$ $\beta^{t+1}_1=\beta^t_1-\alpha\sum _{k=1}^{K}(y_{k}-{\frac {1}{1+e^{-(\beta^t_0+\beta^t_1x_k)}}})x_{k}$

上面求解的是简单的情况，只有一个特征。进一步地，将上面结果推广到具有多个特征的情况。

用矩阵形式表述如下：$X$ 是样本特征矩阵，大小为 $m\times n$，$n$ 是特征维度加1，$m$ 是样本数量，第 $i$ 行形如 $[1, x_1^{(i)}, x_2^{(i)}, x_3^{(i)},\dots, , x_{n-1}^{(i)}]$；$\mathbf{w}$ 是 $n\times1$ 系数向量，形如 $[\beta_1,\beta_2,\dots,\beta_n]^T$；$\mathbf y$ 是 $m\times1$ 列向量，形如 $[y_1,y_2,\dots,y_m]^T$，每个元素取值 $0$ 或 $1$，其中只有一个元素为 $1$；$\mathbf 1^T$ 为值均为 $1$ 的列向量。逻辑函数为

$p(X\mathbf{w})={\frac {1}{1+e^{-X\mathbf{w}}}}$

对数似然函数的矩阵形式为

$\ell(\mathbf{w})=\mathbf y^T\ln p(X\mathbf{w})+(\mathbf 1^T-\mathbf y^T)\ln(\mathbf 1-p(X\mathbf{w}))$

根据逐元素函数法则，逻辑函数的微分

$\begin{aligned} \mathrm{d}p(X\mathbf{w})&=p'(X\mathbf{w})\odot\mathrm{d}X\mathbf{w} \\ &= p(X\mathbf{w})\odot(\mathbf 1-p(X\mathbf{w}))\odot ((\mathrm{d}X)\mathbf{w}+X\mathrm{d}\mathbf{w}) \\ &=p(X\mathbf{w})\odot(\mathbf 1-p(X\mathbf{w}))\odot X\mathrm{d}\mathbf{w} \end{aligned}$

$\ln$ 函数的微分

$\mathrm{d}\ln Y=Y^{-1}\odot\mathrm{d}Y$

因此可以得到

$\begin{aligned} \mathrm{d}\ell(\mathbf{w})&=\mathbf y^T(p(X\mathbf{w}))^{-1}\odot p(X\mathbf{w})\odot(\mathbf 1-p(X\mathbf{w}))\odot X\mathrm{d}\mathbf{w} \\ &+(\mathbf 1^T-\mathbf y^T)(\mathbf 1-p(X\mathbf{w}))^{-1}\odot -p(X\mathbf{w})\odot(\mathbf 1-p(X\mathbf{w}))\odot X\mathrm{d}\mathbf{w} \\ &=\mathbf y^T(\mathbf 1-p(X\mathbf{w}))\odot X\mathrm{d}\mathbf{w} - (\mathbf 1^T-\mathbf y^T)p(X\mathbf{w})\odot X\mathrm{d}\mathbf{w} \\ \end{aligned}$

将上式套上迹

$\begin{aligned} \mathrm{d}\ell(\mathbf{w})&=\mathrm{tr}(\mathbf y^T(\mathbf 1-p(X\mathbf{w}))\odot X\mathrm{d}\mathbf{w}- (\mathbf 1^T-\mathbf y^T)p(X\mathbf{w})\odot X\mathrm{d}\mathbf{w})\\ &=\mathrm{tr}([\mathbf y\odot(\mathbf 1-p(X\mathbf{w}))]^TX\mathrm{d}\mathbf{w} - [(\mathbf 1-\mathbf y)\odot p(X\mathbf{w})]^T X\mathrm{d}\mathbf{w})\\ &=\mathrm{tr}([\mathbf y-\mathbf y\odot p(X\mathbf{w})]^TX\mathrm{d}\mathbf{w} - [p(X\mathbf{w})-\mathbf y\odot p(X\mathbf{w})]^T X\mathrm{d}\mathbf{w})\\ &=\mathrm{tr}(\mathbf y ^ TX\mathrm{d}\mathbf{w} - p(X\mathbf{w})^T X\mathrm{d}\mathbf{w})\\ &=\mathrm{tr}([X^T(\mathbf y- p(X\mathbf{w}) )]^T\mathrm{d}\mathbf{w})\\ \end{aligned}$

根据微分与导数的关系

$\mathrm{d}f=\mathrm{tr}(\frac {\partial f }{\partial X}^T \mathrm{d}X)$

就可以得到

$\frac {\partial f }{\partial \mathbf{w}}=X^T(\mathbf y- p(X\mathbf{w}) )$

迭代公式

$\mathbf{w}=\mathbf{w}-\alpha X^T(\mathbf y- p(X\mathbf{w}) )$

2.3 多元逻辑回归

一个事件的几率（odds）是指该事件的发生的概率与不发生的概率的比值。如果一个事件发生的概率 $p$，则不发生的概率为 $1-p$，那么该事件的几率为 $\displaystyle\frac{p}{1-p}$，该事件的对数几率或 logit 函数是

$\mathrm{logit}(p)=\ln \frac{p}{1-p}$

对于逻辑回归，一定满足

$\ln \frac{p}{1-p}=\mathbf{w}\cdot\mathbf x$

假设多元逻辑回归事件的类别分别为 $1,2,\dots,K$。多元逻辑回归相较于二元逻辑回归，概率计算公式不变

$\displaystyle p(\mathbf x)={\frac {1}{1+e^{-\mathbf{w}\cdot\mathbf x}}}$

区别在于对于同一个实例，它所有类别的概率和为 1。在二元逻辑回归中，$y=0$ 时的概率就是 $1-p$，因此只需要一组参数即可，而在多元逻辑回归中， $K$ 个类别分别对应一组参数

$\ln \frac{p(y=1|\mathbf x)}{1-p(y=0|\mathbf x)}=\mathbf w_1\cdot\mathbf x$ $\ln \frac{p(y=2|\mathbf x)}{1-p(y=0|\mathbf x)}=\mathbf w_2\cdot\mathbf x$ $\cdots$

同时满足

$\sum_{i=1}^Kp(y=i|\mathbf x)=1$

最终解得

$p(y=k|\mathbf x)=\frac{e^{\mathbf w_k\cdot\mathbf x}}{\displaystyle 1+\sum_{i=1}^{K-1}e^{\mathbf w_i\cdot\mathbf x}}$ $p(y=K|\mathbf x)=\frac{1}{\displaystyle 1+\sum_{i=1}^{K-1}e^{\mathbf w_i\cdot\mathbf x}}$

实际上就是 softmax 函数的一种特例，输出的结果为 $k-1$ 个，因为概率和为 1，自然可以求出最后剩下的一个概率值。其实也可以输出 $k$ 个结果，概率的计算公式为

$p(y=k|\mathbf x)=\frac{e^{\mathbf w_k\cdot\mathbf x}}{\displaystyle \sum_{i=1}^Ke^{\mathbf w_i\cdot\mathbf x}}$

我们采用后一种方式，似然函数为

$L=\prod_{i=0}^mp_1(\mathbf{x}_i)^{y_i^1}p_2(\mathbf{x}_i)^{y_i^2}\dots p_k(\mathbf{x}_i)^{y_i^K}$

对于一组输入 $\mathbf{x}_i$，$(y_i^1,y_i^2,\dots,y_i^K)$ 中只有一个值为 $1$，其他为 $0$。对数似然函数即为

$\ell=\sum_{i}\sum_ky_i^k\log p_k(\mathbf{x}_i)$

将对数似然函数写成矩阵形式，对于每一组实例

$\ell=\mathbf{y}^T\log \mathrm{softmax}(W\mathbf{x})$

其中，$\mathbf{y}$ 是除一个元素为 $1$ 外其元素为 $0$ 的 $m\times1$ 列向量，$W$ 为 $m\times n$ 矩阵，$\mathbf{x}$ 为 $n\times1$ 列向量，$\displaystyle\mathrm{softmax}(\mathbf{z})=\frac{\exp(\mathbf{z})}{\mathbf{1}^T\exp(\mathbf{z})}$，其中 $\exp$ 表示逐元素求指数，$\mathbf{1}$ 表示值均为 $1$ 的列向量；$\ell$ 为标量。

将 $\mathrm{softmax}$ 带入 $\ell$

$\begin{aligned} \ell&=\mathbf{y}^T\left(\log \exp(W\mathbf{x})-\mathbf{1}\log \mathbf{1}^T\exp(W\mathbf{x})\right)\\ &=\mathbf{y}^TW\mathbf{x}-\log\left( \mathbf{1}^T\exp(W\mathbf{x})\right) \end{aligned}$

在这里用到了 $\log(\mathbf{u}/c)=\log\mathbf{u}-\mathbf{1}\log c$，以及 $\mathbf{y}^T\mathbf{1}=1$。

进行微分

$\begin{aligned} \mathrm{d}\ell&=\mathrm{d}(\mathbf{y}^TW\mathbf{x})-\mathrm{d}\left(\log( \mathbf{1}^T\exp(W\mathbf{x}))\right)\\ &=\mathbf{y}^T\mathrm{d}W\mathbf{x}-(\mathbf{1}^T\exp(W\mathbf{x}))^{-1}\odot \mathrm{d}(\mathbf{1}^T\exp(W\mathbf{x}))\\ &=\mathbf{y}^T\mathrm{d}W\mathbf{x}-(\mathbf{1}^T\exp(W\mathbf{x}))^{-1}\odot \mathbf{1}^T\mathrm{d}(\exp(W\mathbf{x}))\\ &=\mathbf{y}^T\mathrm{d}W\mathbf{x}-(\mathbf{1}^T\exp(W\mathbf{x}))^{-1}\odot \mathbf{1}^T\left(\exp(W\mathbf{x})\odot\mathrm{d}(W\mathbf{x})\right)\\ \end{aligned}$

其中

$\mathbf{1}^T\left(\exp(W\mathbf{x})\odot\mathrm{d}(W\mathbf{x})\right)=(\mathbf{1}\odot \exp(W\mathbf{x}))^T\mathrm{d}(W\mathbf{x})=\exp(W\mathbf{x})^T\mathrm{d}(W\mathbf{x})$

带入原式并且套上迹

$\begin{aligned} \mathrm{d}\ell&=\mathrm{tr}\left(\mathbf{y}^T\mathrm{d}W\mathbf{x}-\frac{\exp(W\mathbf{x})^T\mathrm{d}(W\mathbf{x})}{\mathbf{1}^T\exp(W\mathbf{x})}\right)\\ &=\mathrm{tr}\left(\mathbf{y}^T\mathrm{d}W\mathbf{x}-\mathrm{softmax}(W\mathbf{x})^T\mathrm{d}(W\mathbf{x})\right)\\ &=\mathrm{tr}\left(\mathbf{x}(\mathbf{y}-\mathrm{softmax}(W\mathbf{x}))^T\mathrm{d}W\right)\\ \end{aligned}$

根据微分与导数的关系

$\mathrm{d}f=\mathrm{tr}(\frac {\partial f }{\partial X}^T \mathrm{d}X)$

因此

$\frac {\partial f }{\partial W}=(\mathbf y- \mathrm{softmax}(W\mathbf{x}))\mathbf{x}^T$

参考

1. Logistic regression - Wikipedia
2. Cramer, J. S. (2002). The origins of logistic regression
3. An Introduction to Logistic Regression: From Basic Concepts to Interpretation with Particular Attention to Nursing Domain
4. Practical Guide to Logistic Regression.pdf)
5. Logistic function - Wikipedia
6. 逻辑回归原理小结 - 刘建平

Machine Learning

2022-02-20