逻辑斯蒂回归(Logistic Regression) | 算法实现_logistic regression算法的实现步骤

01 起

在这篇文章中，我们学习了逻辑斯谛回归模型的算法原理：统计学习方法|logistic回归

今天我们基于算法原理，给出利用随机梯度上升算法求解逻辑斯蒂回归模型参数的过程。

我们先来回顾一下逻辑斯蒂回归模型，

• logistic回归的目的是寻找一个非线性函数sigmoid函数的最佳拟合参数w， sigmoid(wx)=1/(1+exp(-wx))，找到最佳拟合参数w，使不同类别样本点的特征x输入后，被分到对应的分类中，sigmoid值可以认为是一个概率值，从0~1，当>0.5，分到1类，<=0.5分到0类

好了，开始实操。

---

02 梯度下降算法

逻辑斯蒂回归模型的求解过程可以由最优化算法完成，最常用的最优化算法是梯度下降/上升算法，又可以简化为随机梯度下降/上升算法。

一般地，梯度下降算法中的梯度是指模型损失函数对输入变量x的梯度，梯度方向是函数上升最快的方向，于是梯度下降算法是朝着梯度反方向更新系数，可以求得局部最小值，而梯度上升算法则每次朝着梯度方向更新系数，可以求得局部最大值。

如果优化目标是损失函数，那么应该使用梯度下降算法求解使损失函数最小的参数w，而本文的优化目标是二分类模型的似然函数L(w)，该函数表示某样本属于某个类别的概率，所以似然函数越大越好，因此要使用梯度上升算法求似然函数L(w)的极大值。

随机梯度上升算法相比梯度上升算法，每次只使用一个样本点更新参数w，占用更少的计算资源，且可以在新数据到来时就完成参数更新，是一种在线算法，各种梯度下降算法如下：

• 批量梯度下降算法(BGD,Batch Gradient Descent)
• 随机梯度下降算法(SGD,Stochastic Gradient Descent)
• 小批量梯度下降算法(MBGD,Mini-Batch Gradient Descent)

本文先给出梯度上升算法(BGA)，然后给出改进的梯度上升算法——随机梯度上升算法(SGA)

---

03 梯度上升算法实现

思路：

1. 初始化回归系数，即置回归系数w=(w1,w2,…,wn)为全1
2. 将训练集各样本值带入F(X)=sigmoid(w*x)函数计算F(X)值，认为是样本的分类值(二分类01)
3. 计算各训练样本的分类误差error(样本的真实分类值-样本的F(X))
4. 认为本次迭代的梯度=dataMat’*error，得到一个(n,1)的矩阵，用于将回归系数向似然函数增大的方向调整
5. 将回归系数w向似然函数增大的方向调整：w(k+1)=w(k)+a*dataMat’*error
6. 重复2~5步，直到迭代达到指定次数或error小于某值，停止迭代，得到回归系数w，逻辑斯蒂回归模型训练完毕

#sigmoid函数(类阶跃函数，只是0>1斜率较缓，不是瞬变的)
def sigmoid(z):
    return 1.0/(1+np.exp(-z)) # z=w*x

def gradAscent(dataMat,labelMat,numIter=150):
    dataMatrix=np.mat(dataMat) #样本特征矩阵 (m,n)
    labelMatrix=np.mat(labelMat).transpose() #样本类别矩阵 (1,m)
    m,n=np.shape(dataMatrix) #m个样本，n个特征
    
    alpha=0.001 #梯度上升的步长
    weights=np.ones((n,1)) #初始化回归系数
    
    #开始迭代，梯度上升
    for i in range(numIter):
        #将各样本值带入F(X)=sigmoid(w*x)函数计算F(X)，认为是样本在当前迭代的分类值
        Fx=sigmoid(dataMatrix*weights) #shape=(1,m)
        
        #计算各训练样本的分类误差error(样本的真实分类值-样本的F(X))
        #需要将label矩阵浮点数化，与Fx元素数值类型一致，否则报错
        """这里可以加一个if判断，用于在error小于设定值时停止迭代"""
        error=(labelMatrix.astype(float)-Fx) #shape=(1,m) 
        
        #本次迭代的梯度=dataMat'*error,意思为m个样本的n个特征 分别与 m个样本的误差相乘
        """样本误差越大，下一次迭代调整就越多"""
        grad=dataMatrix.transpose()*error #shape=(n,1),用于将回归系数向似然函数增大的方向调整
        
        #将回归系数w向似然函数增大的方向调整：w(k+1)=w(k)+a*dataMat'*error
        weights=weights+alpha*grad
  
    #直到迭代达到指定次数或error小于某值，停止迭代，得到回归系数w，逻辑斯蒂回归模型训练完毕
    return weights
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31

我们使用一组二维分类数据集来测试一下算法，得到的分类边界如下图所示，其中蓝点红点分别表示数据集中划分好的不同类别，橙色线条为决策边界，是训练得到的逻辑斯蒂回归模型，

---

04 随机梯度上升算法实现

以上计算可以看到，梯度上升法计算量较大，即使是简单数据集计算量也很大（100条3维度数据需要进行300次乘法，再算上1000此迭代，计算量就更大了），当面对上万上亿条数据，计算复杂度就太大了。

于是我们给出改进的梯度上升算法——随机梯度上升算法，

• 一次仅用一个样本点来更新回归系数w
• 此方法可以在新样本到来时进行增量式的更新，因此随机梯度上升法是一种在线学习算法

随机梯度上升法思路：与梯度上升法类似，不同之处在于：

1. 随机梯度上升法一次只用一个样本点来更新回归系数w
2. 需要遍历所有的样本点来更新回归系数w，然后才相当于梯度上升法的一次迭代
3. 相比于梯度上升法，随机梯度上升法每次迭代只有100次(样本数)乘法，复杂度大大降低，特别对于高维样本

def stocGradAscent(dataMat,labelMat):
    m,n=np.shape(dataMat) #m个样本，n个维度(特征)，dataMat为list
    dataMatrix=np.mat(dataMat) #list转换为矩阵，便于后续计算
    alpha=0.01 #步子迈大点
    weights=np.ones((n,1)) #初始化回归系数，一维数组
    
    for i in range(m): #遍历所有样本点
        #第i个样本点的F(x)值,float内为第i个样本点的w*x
        Fxi=sigmoid(float((dataMatrix[i]*weights)))
        
        #第i个样本的真实值与预测值误差，用于回归系数调整
        error=float(labelMat[i])-Fxi
        
        #回归系数调整，只用第i个样本点更新回归系数w
        weights=weights+alpha*error*dataMatrix[i].transpose()
    return weights
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

我们使用相同的数据集测试一下随机梯度上升算法训练的逻辑斯蒂回归模型，如下图，

咦，怎么感觉分类能力更弱了？

不不不，再仔细看看，随机梯度上升算法只迭代了1次，而刚才我们迭代了1000次，这样对比好像不太公平呐！

我们再改进一下，

1. 加入多次迭代
2. 每次迭代的步长(alpha)会被调整,越往后迭代，步长越小
   • 可缓解迭代过程中系数的波动，前几次迭代中，步长较大，即，alpha不是严格下降的，便于系数快速收敛
3. 每次迭代中，遍历样本不是按顺序遍历的，而是随机遍历的，即，随机选取样本来更新回归系数w，每次随机从样本集中选出一个样本，更新系数，然后删除该样本点，进行下一次i迭代
   • 可减少迭代过程中回归系数的周期性波动，因为周期性波动(波形图中的小锯齿)来自于一些不能被正确分类的样本点（数据集并非线性可分），这些样本点每次被用于更新系数时会引起系数的剧烈波动

def stocGradAscentIter(dataMat,labelMat,numIter=150):
    m,n=np.shape(dataMat) #m个样本，n个维度(特征)，dataMat为list
    dataMatrix=np.mat(dataMat) #list转换为矩阵，便于后续计算
    weights=np.ones((n,1)) #初始化回归系数，一维数组

    for j in range(numIter):
        dataIndex=list(range(m))#样本点索引，用于存储尚未用于更新系数的样本点的索引
        for i in range(m): #遍历所有样本点
            """第j次迭代的第i个小迭代"""
            alpha=4/(1.0+j+i)+0.01 #每次迭代的步长(alpha)会被调整,越往后迭代，步长越小

            #随机选取样本点进行系数更新
            randIndex=np.random.randint(0,len(dataIndex)) #从尚未用于更新系数的样本点中随机选取一个样本
            Fxi=sigmoid(float((dataMatrix[randIndex]*weights)))
            error=float(labelMat[randIndex])-Fxi
            weights=weights+alpha*error*dataMatrix[randIndex].transpose()
            
            #删除已经用于更新系数的样本点的索引
            del(dataIndex[randIndex])
    return weights
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

再来看看训练效果，

对比一下使用梯度上升算法迭代1000次后的模型，改进的随机梯度上升算法只迭代了10次就已经展现出较好的分类能力了，溜~~

---

05 总结

本文给出了随机梯度下降/上升算法的原理和方法，并基于随机梯度上升算法求解逻辑斯蒂回归模型的参数w，完成模型训练，旨在加深对逻辑斯谛回归模型和随机梯度下降算法的进一步理解，希望对你有帮助=.=

---

06 参考

1. 《统计学习方法》 李航 Chapter6
2. 《机器学习实战》 Peter Harrington Chapter5