LR(Logistic Regression)算法详解

Logistic Regression本质上还是Linear Regression的一种，只是用了一个Logistic Function将线性回归的连续值映射到了 { 0 , 1 } \{0, 1\} {0,1}空间。因此Linear Regression只能对具有线性边界的分类问题有很好的预测效果，对于非线性的边界是无能为力的。至于下面这张很经典的大家常用的图，只是做了一个feature mapping，根据已有的特征构造了其他的很多新的特征，跟Logistic Regression没有任何关系。而feature mapping的工作则是一项专门的工作，绝没有这么简单。

因此说白了，Logistic Regression就是试图找到不同类别之间的线性决策边界。

1. Logistic Regression的引出

对于 { 0 , 1 } \{0,1\} {0,1}二分类问题，使用Linear Regression难以模拟出一个合适的模型来拟合数据，因此我们试图采用另一种方法来拟合。对于二分类，其实就是某事件发生或者不发生。因此我们希望能计算出某事件发生的概率和不发生的概率。假设某事件发生的概率为$p$，则不发生的概率为$1-p$，定义一个让步值$odds$为：
$$odds=\frac{P(occuring)}{P(notoccuring)}=\frac{p}{1-p}$$
再定义一个让步比$OR(oddsratio)$为两个事件$s_1$和$s_2$的$odds$的比值：
$$OR=\frac{\frac{p_1}{1-p_1}}{\frac{p_0}{1-p_0}}$$
其中$p_1$是实验组的事件发生概率，而$p_0$是对照组(一般事先已知)的事件发生概率。
在Logistic Regression问题中，因变量$y$服从伯努利分布(又称二项分布)，事件发生概率为$p$。在Logistic Regression中，我们需要根据给定的一堆自变量，模拟出一个现象模型，来预测事件发生的概率$p$。因此我们需要让自变量的分布和伯努利分布联系起来，这被称为$logit$。Logistic Regression的目标就是在统计学上预估出事件发生的概率$\hat{p}$。

基于以上思想，为了将自变量的线性模型和伯努利分布联系起来，我们需要一个函数来将线性模型映射到一个 { 0 , 1 } \{0, 1\} {0,1}伯努利分布上。因此上面提到的$OR$的自然对数，就是我们需要的这样的一个完美的映射，也被称为$logit$。
$$\ln (odds)=\ln (\frac{p}{1-p})$$
此时$p$的范围是$(0,1)$，$\ln (odds)$的范围是$(-\infty ,+\infty )$。我们需要求一下反函数：
$$p=logi{t}^{-1}(\alpha )=\frac{1}{1+e^{-\alpha }}$$
$\alpha$范围$(-\infty ,+\infty )$。此时的$\alpha$就是自变量的线性组合。图像为：

也就是著名的$sigmod$函数，此函数极有利于数学计算。然后令：
$$\ln \frac{p}{1-p}={\theta }_0+{\theta }_1{x}_1+{\theta }_2{x}_2+\cdots +{\theta }_j{x}_j$$
其中$j$为特征的总数量。选择这个式子的原因是：假设我们的样本拥有线性决策边界(这是Logistics Regression的前提)，当一个样本的特征$x^{(j)}$的特征给出时，能够通过${\theta }_0+{\sum }_{i=1}^j{\theta }_j{x}_j$很强烈的判断出该样本属于1类还是0类。因此此时$p$最好无限接近1或者0，因此$\frac{p}{1-p}$无限接近$+\infty$或者$-\infty$，而$\ln \frac{p}{1-p}$无限接近1或者0，刚好能够对应起来。另外选择$e$为底数，是为了方便计算。

2. 极大似然估计($MLE$)

对于给定的一堆样本点，若已知其为第1类，但是我们的模型并不会算出是1，而是会根据自变量得出一个概率值。此概率为：$P(Y=y_i=1|X=x_i)$。若已知其为第0类，模型预测的概率值为：$P(Y=y_i=0|X=x_i)$。好的模型当然是两个概率尽可能越大越好，因此我们把它们乘起来，求其最大值对应的参数，求出来的就是最佳拟合模型。因此极大似然估计其实就是使得我们的模型正确地分类数据集的可能性到达最大。
$$P=\prod _{y_i=1}P(Y=y_i=1|X=x_i)\times \prod _{y_i=0}P(Y=y_i=0|X=x_i)=\prod _{i=1}^{n}P(Y=y_i=1|X=x_i{)}^{y_i}(1-P(Y=y_i=1|X=x_i){)}^{1-y_i}$$
因此逻辑回归就是要求$P$的最大值。这就是极大似然估计。
假设概率$P$个自变量$X$之间存在关系：
$$\begin{gathered} P(y=1|x;\theta )=h_{\theta }(x) \\ P(y=0|x;\theta )=1-h_{\theta }(x) \end{gathered}$$
合并为：
$$P(y|x;\theta )=h_{\theta }^y(x)(1-h_{\theta }(x){)}^{1-y}$$
损失函数为:
$$L(\theta )=\prod _{i=1}^m{h}_{\theta }(x^{(i)}{)}^{y^{(i)}}(1-h_{\theta }(x^{(i)}){)}^{1-y^{(i)}}$$
这个式子计算不方便，我们对其进行求对数：
$$\begin{array}{rl}l(\theta )& =\ln (L(\theta ))\\ & =\sum _{i=1}^m[y^{(i)}\ln (h_{\theta }(x^{(i)}))+(1-y^{(i)})\ln (1-h_{\theta }(x^{(i)}))]\end{array}$$
然后对$l(\theta )$求偏导：
$$\frac{\partial (l(\theta ))}{\partial {\theta }_j}=\sum _{i=1}^m(y^{(i)}-h_{\theta }(x^{(i)})\cdot x_j^{(i)})$$
由于这里求最大值，因此是加号：
$${\theta }_j:={\theta }_j+\alpha \sum _{i=1}^m(y^{(i)}-h_{\theta }(x(i)))\cdot x_j^{(i)}$$
通过用类似线性回归的不断迭代，最终可以求出$\vec{\theta }$的值(具体查看线性回归的求解，此处略去)。求解方法其实有很多，最常用的用是牛顿法。

3. 正则化

为了避免模型的过拟合，需要加上正则项，则：
$$l^{\prime }(\theta )=l(\theta )+regularization$$
常用的正则化方法有2个，分别是

• L1
  $$\theta =\arg \underset{\theta }{\min }\sum _{i=1}^m[y^{(i)}\ln (h_{\theta }(x^{(i)}))+(1-y^{(i)})\ln (1-h_{\theta }(x^{(i)}))]+\lambda \sum _{i=0}^k|{\theta }_i|$$

• L2
  $$\theta =\arg \underset{\theta }{\min }\sum _{i=1}^m[y^{(i)}\ln (h_{\theta }(x^{(i)}))+(1-y^{(i)})\ln (1-h_{\theta }(x^{(i)}))]+\lambda \sum _{i=0}^k{\theta }_i^2$$

适用情况如下:

L2	L1
由于有分析解决方案，计算效率高	数据不稀疏的时候计算效率不高
产生非稀疏的输出	产生稀疏的输出
没有特征选择	有内置的特征选择

4. 牛顿法求解

牛顿法是一种进行$MLE$的优化方法，牛顿法在求解逻辑回归时比梯度下降法收敛速度快。

4.1 Hessian

Hessian是一个多变量实值函数的二阶偏导数组成的方块矩阵，假设有一实数函数$f(x_1,x_2,\cdots ,x_m)$。如果 $f$ 所有的二阶偏导数都存在，那么$f$的海森矩阵的第$ij$项即：
$$H(f{)}_{ij}(x)=D_i{D}_{j}f(x)$$
其中$x=(x_1,x_2,\cdots ,x_m)$，即
$$H(f)=\left[\begin{array}{c}\frac{{\partial }^{2}f}{\partial x_1^2}\frac{{\partial }^{2}f}{\partial x_1{x}_2}\cdots \frac{{\partial }^{2}f}{\partial x_1{x}_m}\\ \frac{{\partial }^{2}f}{\partial x_2{x}_1}\frac{{\partial }^{2}f}{\partial x_2^2}\cdots \frac{{\partial }^{2}f}{\partial x_2{x}_m}\\ ⋮\\ \frac{{\partial }^{2}f}{\partial x_m{x}_1}\frac{{\partial }^{2}f}{\partial x_m{x}_2}\cdots \frac{{\partial }^{2}f}{\partial x_m^2}\end{array}\right]$$

Hessian matrix在用牛顿法解决大规模优化问题时经常用到。

4.2 牛顿法步骤：

不断进行下面的迭代:
$${\beta }^{new}={\beta }^{old}-Hf(\beta {)}^{-1}\nabla f(\beta )$$
其中$Hf(\beta )=-X^{T}WX$，为函数$f(x)$的海森矩阵。$\nabla f(\beta )=X^T(y-p)$是函数的梯度。$W:=diag(p(1-p))$是权重，$p$是在当前$\beta$下的概率预测值。代码为：

def new_step(curr_beta, X, lam=None):
	p = np.array(sigmod(X.dot(curr[:,0])), ndmin=2).T #计算新的概率p
    
    W = np.diag((p(1-p))[:, 0]) #计算新的权重
    
    hessian = X.T.dot(W).dot(X) #计算当前权重的hessian矩阵
    
    grad = X.T.dot(y-p) #计算梯度
    
    if lam:
    	#正则化
    	step, *_ = np.linalg.lstsq(hessian + lam*np.eye(curr_beta.shape[0]), grad)
    else:
    	step, *_ = np.linalg.lstsq(hessian, grad)
    
    beta_new = curr_beta + step
    
    return beta_new

beta_old, beta = np.ones((len(X.columns),1)), np.zeros((len(X.columns),1))

coefs_converged = False

while not coefs_converges:
	beta_old = beta
    
    beta = newton_step(beta, X, lam=lam)
    
    coef_converged = check_coefs_convergence(beta_old, beta, tolerance, itercount)