LR(Logistic Regression)算法详解

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

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

1. Logistic Regression的引出

对于{0,1}二分类问题，使用Linear Regression难以模拟出一个合适的模型来拟合数据，因此我们试图采用另一种方法来拟合。对于二分类，其实就是某事件发生或者不发生。因此我们希望能计算出某事件发生的概率和不发生的概率。假设某事件发生的概率为p，则不发生的概率为1−p，定义一个让步值odds为：

odds=P(occuring)P(notoccuring)=p1−p

再定义一个让步比OR(oddsratio)为两个事件s1和s2的odds的比值：

OR=p11−p1p01−p0

其中p1是实验组的事件发生概率，而p0是对照组(一般事先已知)的事件发生概率。
在Logistic Regression问题中，因变量y服从伯努利分布(又称二项分布)，事件发生概率为p。在Logistic Regression中，我们需要根据给定的一堆自变量，模拟出一个现象模型，来预测时间发生的概率p。因此我们需要让自变量的分布和伯努利分布联系起来，这被称为logit。Logistic Regression的目标就是在统计学上预估出事件发生的概率p^。

基于以上思想，为了将自变量的线性模型和伯努利分布联系起来，我们需要一个函数来将线性模型映射到一个{0,1}伯努利分布上。因此上面提到的OR的自然对数，就是我们需要的这样的一个完美的映射，也被称为logit。

ln(odds)=ln(p1−p)

此时p的范围是(0,1)，ln(odds)的范围是(−∞,+∞)。我们需要求一下反函数：

p=logit−1(α)=11+e−α

α范围(−∞,+∞)。此时的α就是自变量的线性组合。图像为：

也就是著名的sigmod函数，此函数极有利于数学计算。然后令：

lnp1−p=θ0+θ1x1+θ2x2+⋯+θjxj

其中j为特征的总数量。选择这个式子的原因是：假设我们的样本拥有线性决策边界(这是Logistics Regression的前提)，当一个样本的特征x(j)的特征给出时，能够通过θ0+∑ji=1θjxj很强烈的判断出该样本属于1类还是0类。因此此时p最好无限接近1或者0，因此p1−p无限接近+∞或者−∞，而lnp1−p无限接近1或者0，刚好能够对应起来。另外选择e为底数，是为了方便计算。

2. 极大似然估计(MLE)

对于给定的一堆样本点，若已知其为第1类，但是我们的模型并不会算出是1，而是会根据自变量得出一个概率值。此概率为：P(Y=yi=1|X=xi)。若已知其为第0类，模型预测的概率值为：P(Y=yi=0|X=xi)。好的模型当然是两个概率尽可能越大越好，因此我们把它们乘起来，求其最大值对应的参数，求出来的就是最佳拟合模型。因此极大似然估计其实就是使得我们的模型正确地分类数据集的可能性到达最大。

P=∏yi=1P(Y=yi=1|X=xi)×∏yi=0P(Y=yi=0|X=xi)=∏i=1nP(Y=yi=1|X=xi)yi(1−P(Y=yi=1|X=xi))1−yi

因此逻辑回归就是要求P的最大值。这就是极大似然估计。
假设概率P个自变量X之间存在关系：

P(y=1|x;θ)=hθ(x)P(y=0|x;θ)=1−hθ(x)

合并为：

P(y|x;θ)=hyθ(x)(1−hθ(x))1−y

损失函数为:

L(θ)=∏i=1mhθ(x(i))y(i)(1−hθ(x(i)))1−y(i)

这个式子计算不方便，我们对其进行求对数：

l(θ)=ln(L(θ))=∑i=1m[y(i)ln(hθ(x(i)))+(1−y(i))ln(1−hθ(x(i)))]

然后对l(θ)求偏导：

∂(l(θ))∂θj=∑i=1m(y(i)−hθ(x(i))⋅x(i)j)

由于这里求最大值，因此是加号：

θj:=θj+α∑i=1m(y(i)−hθ(x(i)))⋅x(i)j

通过用类似线性回归的不断迭代，最终可以求出θ⃗ 的值(具体查看线性回归的求解，此处略去)。求解方法其实有很多，最常用的用是牛顿法。

3. 正则化

为了避免模型的过拟合，需要加上正则项，则：

l′(θ)=l(θ)+regularization

常用的正则化方法有2个，分别是
+ L1

θ=argminθ∑i=1m[y(i)ln(hθ(x(i)))+(1−y(i))ln(1−hθ(x(i)))]+λ∑i=0k|θi|

• L2
  
  θ=argminθ∑i=1m[y(i)ln(hθ(x(i)))+(1−y(i))ln(1−hθ(x(i)))]+λ∑i=0kθ2i

适用情况如下:

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

4. 牛顿法求解

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

4.1 Hessian

Hessian是一个多变量实值函数的二阶偏导数组成的方块矩阵，假设有一实数函数f(x1,x2,⋯,xm)。如果 f 所有的二阶偏导数都存在，那么f的海森矩阵的第ij项即：

H(f)ij(x)=DiDjf(x)

其中x=(x1,x2,⋯,xm)，即

H(f)=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢∂2f∂x21∂2f∂x1x2⋯∂2f∂x1xm∂2f∂x2x1∂2f∂x22⋯∂2f∂x2xm⋮∂2f∂xmx1∂2f∂xmx2⋯∂2f∂x2m⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥

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

4.2 牛顿法步骤：

不断进行下面的迭代:

βnew=βold−Hf(β)−1∇f(β)

其中Hf(β)=−XTWX，为函数f(x)的海森矩阵。∇f(β)=XT(y−p)是函数的梯度。W:=diag(p(1−p))是权重，p是在当前β下的概率预测值。代码为：

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_newbeta_old, beta = np.ones((len(X.columns),1)), np.zeros((len(X.columns),1))coefs_converged = Falsewhile not coefs_converges:    beta_old = beta    beta = newton_step(beta, X, lam=lam)    coef_converged = check_coefs_convergence(beta_old, beta, tolerance, itercount)