（学习笔记）Logistic Regression 对数几率回归

介绍

Logistic Regression 对数几率回归，虽然被称为回归，但其实际上是一种分类学习方法，并常用于二分类。
它有很多优点，例如它是直接对分类可能性进行建模，无需事先假设数据分布；它不是仅预测出类别，而是可得到近似概率预测，这对许多需利用概率辅助决策的任务很有用；此外，对率回归求解的目标函数是任意阶可导的凸函数，有很好的数学性质。

原理

考虑二分类任务，其输出标记 y ∈ { 0 , 1 } y\in\{0,1\} y∈{0,1}，而线性回归产生的预测值$z={\omega }^{T}x+b$是实值，我们需要将实值$z$转换为0/1值，最理想的是“单位阶跃函数”(unit-step function)
$$y=\{\begin{array}{ll}0& \text{z<0}\\ 0.5& \text{z=0}\\ 1& \text{z>0}\end{array}$$
即若预测值$z$大于0就判为正例，小于0则判为反例，预测值为临界值则可随意判别。
但由于单位阶跃函数不连续，不能作为广义线性模型$y=g^{-1}({\omega }^{T}x+b)$中的$g^{-1}(\ast )$，我们希望找到在一定程度上近似单位阶跃函数的替代函数，并且单调可微，对数几率函数(logistic function)正是这样一个常用的替代函数：
$$y=\frac{1}{1+e^{-z}}$$

其可变化为：
$$ln\frac{y}{1-y}={\omega }^{T}x+b$$
若将$y$视为样本$x$作为正例的可能性，$\frac{y}{1-y}$称为“几率”(odds)，取对数即得“对数几率”(log odds / logit)：
$$ln\frac{y}{1-y}$$

将$y$视为类后验概率估计$p(y=1|x)$，则式子重写为：
$$ln\frac{p(y=1|x)}{p(y=0|x)}={\omega }^{T}x+b$$
显然有：
$$\begin{gathered} p(y=1|x)=\frac{e^{{\omega }^T+b}}{1+e^{{\omega }^T+b}} \\ p(y=0|x)=\frac{1}{1+e^{{\omega }^T+b}} \end{gathered}$$
通过极大似然估计法来估计$\omega 和b$:
$$L(\omega ,b)=\prod p(y_i|x_i;\omega ,b)$$
即
$$L(\omega ,b)=\prod [{p(y=1|x_i;\omega ,b)}^{y_i}{p(y=0|x_i;\omega ,b)}^{1-y_i}]$$
令$\beta =(\omega ;b),\hat{x}=(x;1)$,$p({\hat{x}}_i)=p(y=1|{\hat{x}}_i),$两边取对数,
$$\begin{array}{rl}lnL(\omega )& =\sum [ln[p(\widehat{x_i}){]}^{y_i}+ln[(1-p(\widehat{x_i}){]}^{1-y_i}]\\ & =\sum [y_{i}lnp(\widehat{x_i})+(1-y_i)ln[1-p(\widehat{x_i})]]\\ & =\sum [y_{i}ln\frac{e^{{\beta }^T\widehat{x_i}}}{1+e^{{\beta }^T\widehat{x_i}}}+(1-y_i)ln(\frac{1}{1+e^{{\beta }^T\widehat{x_i}}})]\\ & =\sum [y_i({\beta }^T\widehat{x_i}-ln(1+e^{{\beta }^T\widehat{x_i}})+(1-y_i)(-ln(1+e^{{\beta }^T\widehat{x_i}})]\\ & =\sum [y_i{\beta }^T{x}_i-ln(1+e^{{\beta }^T\widehat{x_i}})]\end{array}$$最后要求解的问题就是
$${\beta }^{\ast }=argminJ(\omega )=\sum [-y_i{\beta }^T{\hat{x}}_i+ln(1+e^{{\beta }^T\widehat{x_i}})]$$
$J(\omega )$是关于$\beta$的高阶可导的连续凸函数，使用经典的数值优化方法，如梯度下降法、牛顿法等都可以求得其最优解。
以牛顿法为例：
Newton法步骤如下:
$$\begin{array}{rl}& 对于x^{\ast }=argminf(x)\\ (1)& 取初始点x^1,置精度要求ϵ,置k=1.\\ (2)& if||▽f(x)||\le ϵ，则停止计算；否则求解{▽}^{2}f(x^k)d=-▽f(x^k),得到d^k.\\ (3)& 置x^{k+1}=x^k+d^k,k=k+1,转步骤(2).\end{array}$$Newton法要求$f(x)$是二阶可微的。
对于$J(\omega ),$
$$\begin{array}{rl}& ▽J(\omega )=\sum {\hat{x}}_i(p({\hat{x}}_i)-y_i)\\ & {▽}^{2}J(\omega )=\sum {\hat{x}}_i{\hat{x}}_i^{T}p({\hat{x}}_i)(1-p({\hat{x}}_i))\end{array}$$
以上，就是对数几率回归的大致过程了。

Python代码实现

下面以西瓜书P89的西瓜数据集3.0$\alpha$为例，使用Python编程实现对率回归

import numpy as np

def Sigmoid(z):
    return 1 / (1 + np.exp(-z))

def Graninet(X, Y, beta):  # 求梯度
    y = np.array(Sigmoid(np.dot(beta.T, X.T)))[0]
    m, _ = X.shape
    d = np.zeros((3, 1))
    for i in range(m):
        d += X[i].T * (y[i] - Y[i])
    return d

def Hesse(X, beta):  # 黑塞矩阵
    y = np.array(Sigmoid(np.dot(beta.T, X.T)))[0]
    m, _ = X.shape
    dd = np.zeros((3, 3))
    for i in range(m):
        dd += np.dot(X[i].T, X[i]) * y[i] * (1 - y[i])
    return np.mat(dd)

def Newton(X, Y, beta, epsi):  # 牛顿法
    d = Graninet(X, Y, beta)
    dd = Hesse(X, beta)
    while np.linalg.norm(d) > epsi:
        beta = beta + np.dot(dd.I, -d)
        d = Graninet(X, Y, beta)
        dd = Hesse(X, beta)
    return beta

def main():
    X = np.mat([[0.697, 0.460, 1], [0.774, 0.376, 1], [0.634, 0.264, 1], [0.608, 0.318, 1], [0.556, 0.215, 1],
                [0.403, 0.237, 1], [0.481, 0.149, 1], [0.437, 0.211, 1], [0.666, 0.091, 1], [0.243, 0.267, 1],
                [0.245, 0.057, 1], [0.343, 0.099, 1], [0.639, 0.161, 1], [0.657, 0.198, 1], [0.360, 0.370, 1],
                [0.593, 0.042, 1], [0.719, 0.103, 1]])  # 数据集
    Y = np.array([1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0])
    epsi = 0.00001  # 精度
    beta = np.random.rand(3, 1)
    beta = Newton(X, Y, beta, epsi)
    print(beta)

if __name__ == '__main__':
    main()
   
#result:	[[ 3.15832262]
 			 [12.52116868]
 			 [-4.42885484]]

参考资料

[1]周志华.机器学习[M].北京：清华大学出版社，2016：53-60.
[2]王燕军，梁治安，崔雪婷.最优化基础理论与方法[M].上海：复旦大学出版社，2018：48-49