1. 逻辑回归的模型函数
前面我们讲了线性回归模型,将线性模型用于回归问题中。这篇我们讲一下线性模型用于分类任务。在二分类问题中,对于线性回归所产生的预测值
我们需将这个预测z转化为0/1值,最理想的是“单位阶跃函数”,即若预测z大于零就判为1,若预测z小于零则判为反例,预测值为临界值则可以任意判别。但是由于单位阶跃函数不连续,我们希望可以找到一个连续的,同时在一定程度上近似“单位阶跃函数“的函数,那么逻辑斯特函数(logistic function)则是一个很不错的替代函数:
如下图所示,这是一种"Sigmoid"函数,它z值转化为一个接近0或1的y值。
![](https://i-blog.csdnimg.cn/blog_migrate/a87fabfa9eca2f0e092b90145e25069c.png)
我们将
代入上式,最终得到逻辑回归的模型函数,如下
作者个人理解:“单位阶跃函数”是将预测值直接转换为分类,而“逻辑斯特函数”是将预测值转换为0-1之间的一个数y,这个值正好可以当成某一类别的概率值,从而实现分类。
2. 逻辑回归的目标函数
逻辑回归的模型函数已经知道,求逻辑回归的模型函数,其实就是确定该模型函数theta这个参数的值。在讲如何求这个参数之前,需要求得逻辑回归的目标函数,而在求解目标函数的过程中,需先讲解一下以下知识点:
一个事件的几率(odds):指该事件发生与不发生的概率比值,若事件发生概率p,那么事件发生的几率就是
那么该事件的对数几率(log odds,亦称logit)就是:
将式子
转化为如下
也就是说,输出y的对数几率是由输入x的线性函数表示的模型,这就是逻辑回归模型。当
的值越接近正无穷,y的值是越接近1的,那么我们将y视为样本x为正例,即y=1的概率。那么上述公式可转化为如下
从而得
进一步可转化为
在统计学中,常常使用极大似然估计法来求解,即找到一组参数,使得在这组参数下,我们的数据的似然度(概率)最大,似然函数如下:
取对数似然函数得
我们要最大化上式,也就等价于最小化
上式即可当做逻辑回归的目标函数
在逻辑回归模型中,我们最大化似然函数和最小化对数似然损失函数实际上是等价的。
3. 求解目标函数
上一块知识讲解中已经得到了目标函数,这一块知识讲解中,主要是讲解如何得到目标函数的最小值。自然而来的还是想到了梯度下降算法。我们求上述目标函数的梯度得
其中
那么对于
的一次迭代如下:
从而推出对于theta的一次迭代如下(并进行了向量化)
其中
上面将求和公式转化为了向量的形式,具体推导的过程已省略,可自行推导,也可以参考文末的梯度下降算法的向量化推导的参考链接。
4. 代码的实现
上面讲了逻辑回归的理论,下面我们通过代码来实现逻辑回归,例子主要是机器学习实战书上的例子。首先来看一下部分数据集(数据集的地址见文末给的链接地址),如下
-0.017612 14.053064 0
-1.395634 4.662541 1
-0.752157 6.538620 0
-1.322371 7.152853 0
...
-2.168791 0.143632 1
1.388610 9.341997 0
0.317029 14.739025 0
我们需要从文件中读取数据,并返回读取之后的列表结果,实现的代码如下
'''
从文件中读取数据
'''
def loadDataSet(filePath):
dataList = []
labelList = []
f = open(filePath)
for line in f.readlines():
lineArr = line.strip().split()
dataList.append([1.0, float(lineArr[0]), float(lineArr[1])])
labelList.append(int(lineArr[2]))
return dataList, labelList
之后我们采用梯度下降算法进行求解,在迭代过程中使用向量化,实现的代码如下:
'''
sigmoid函数
'''
def sigmoid(inX):
return 1.0 / (1 + np.exp(-inX))
'''
梯度下降函数
'''
def gradAscent(dataList, labelList):
dataMat = np.mat(dataList) # x数据转化为矩阵
labelMat = np.mat(labelList).T # y数据转换为矩阵并转置为列向量
m, n = np.shape(dataMat) # 返回矩阵的大小
alpha = 0.001 # 步长
maxCycles = 500 # 迭代次数
weights = np.ones((n, 1)) # 权重列向量
# 梯度下降算法
for k in range(maxCycles):
h = sigmoid(dataMat * weights)
error = h - labelMat
weights = weights - alpha * dataMat.T * error
return weights
通过函数gradAscent
可以得到一组回归系数,它确定了不同类别数据之间的分隔线,下面我们将这个分割线画出来,代码如下:
'''
画分类的分割线
'''
def plotBestFit(dataList, labelList, weights):
dataArr = np.array(dataList)
m = np.shape(dataArr)[0] # 获取样本的数量
xcord0 = []
ycord0 = []
xcord1 = []
ycord1 = []
# 对不同标记的数据进行分类
for i in range(m):
if int(labelList[i]) == 0:
xcord0.append(dataArr[i, 1])
ycord0.append(dataArr[i, 2])
else:
xcord1.append(dataArr[i, 1])
ycord1.append(dataArr[i, 2])
# 画散点图
plt.scatter(xcord0, ycord0, s=30, c='red', marker='s')
plt.scatter(xcord1, ycord1, s=30, c='green')
# 画分割线
x1 = np.arange(-3.0, 3.0, 0.1)
x2 = (-weights[0] - weights[1] * x1) / weights[2]
plt.plot(x1, x2)
plt.xlabel('X1')
plt.ylabel('X2')
plt.show()
上述代码中,比较令人不解的是x2 = (-weights[0] - weights[1] * x1) / weights[2]
这个是怎么来的?因为对于样本数据来说
是两个分类(类别1和类别0)的分界处,根据该式代入x_1即可求x_2从而这条线可以在图中表示出来。运行如下代码:
dataList, labelList = loadDataSet("testSet.txt")
weights = gradAscent(dataList, labelList)
plotBestFit(dataList, labelList, weights.getA())
输出结果如下图所示,从图中可以看到分类效果相当不错,从图上看值只分错了几个点。
![](https://i-blog.csdnimg.cn/blog_migrate/62dc01488030490ddcae803c5f089c90.png)
然而跟线性回归中使用梯度下降算法求最小值类似,这种方法需要大量的计算,每次更新回归系数都需要遍历整个数据集,假如数据量一多,该方法的计算复杂度就很高了。所以我们对上述的批量梯度下降算法进行改变,改进方法是一次仅用一个样本点来更新回归系数,该方法称为随机梯度下降算法。下面是随机梯度下降算法的实现代码
'''
随机梯度下降算法
'''
def stocGradAscent0(dataList, labelList):
dataArr = np.array(dataList) # x数据转化为矩阵
m, n = np.shape(dataArr)
alpha = 0.01
weights = np.ones(n)
for i in range(m):
h = sigmoid(np.sum(dataArr[i]*weights))
error = (h-labelList[i])
weights = weights-alpha*error*dataArr[i]
return weights
运行如下代码,得到了该方法的分类效果图:
dataList, labelList = loadDataSet("testSet.txt")
weights = stocGradAscent0(dataList, labelList)
plotBestFit(dataList, labelList, weights)
![](https://i-blog.csdnimg.cn/blog_migrate/309c7daf20068bbf7c3823647d9fbbe6.png)
相比批量梯度下降算法而言,效果较差。直接比较效果是有点不公平的,因为两者使用的数据集次数不同。而判断一个优化算法优劣的可靠方法是看它是否收敛,也就是说参数是否达到稳定值。通过下面这段代码,来绘制三个回归系数的变化情况:
'''
记录随机梯度下降过程中的历史权重值
'''
def stocGradAscent0WeightHistory(dataList, labelList):
dataMat = np.array(dataList) # x数据转化为矩阵
m, n = np.shape(dataMat)
alpha = 0.01
times = 100
weights = np.ones(n)
weightsHistory = np.zeros((times*m, n) )
for j in range(times):
for i in range(m):
h = sigmoid(np.sum(dataMat[i] * weights))
error = (h - labelList[i])
weights = weights - alpha * error * dataMat[i]
weightsHistory[j*m+i, :] = weights
return weightsHistory
'''
绘制历史权值图
'''
def plotSDGError(weightsHistory):
plt.rcParams['font.sans-serif'] = ['SimHei'] # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False # 用来正常显示负号
fig = plt.figure()
ax = fig.add_subplot(311)
ax.plot(weightsHistory[:,0])
plt.ylabel('X0')
ax = fig.add_subplot(312)
ax.plot(weightsHistory[:,1])
plt.ylabel('X1')
ax = fig.add_subplot(313)
ax.plot(weightsHistory[:,2])
plt.ylabel('X2')
plt.xlabel(u'迭代次数')
plt.show()
def main():
dataList, labelList = loadDataSet("testSet.txt")
weightsHistory = stocGradAscent0WeightHistory(dataList, labelList)
plotSDGError(weightsHistory)
if __name__ == '__main__':
main()
最终绘制的三个回归系数的变化情况如下:
![](https://i-blog.csdnimg.cn/blog_migrate/46e4f13c3e63e68dd342bc47360d0d87.png)
从图中可以看到在大的波动之后,有些系数还会存在一些小的周期性波动。这种现象的原因是因为存在一些不能正确分类的样本点,在每次迭代时会引发系数的剧烈改变。我们希望算法能避免来回波动,从而收敛于某个值,同时是加快收敛速度,因此我们将上述的随机梯度下降算法进行更改。更改为如下:
'''
改进的随机梯度下降算法
'''
def stocGradAscent1(dataList, labelList, numIter=150):
dataArr = np.array(dataList)
m,n = np.shape(dataArr)
weights = np.ones(n)
for j in range(numIter):
dataIndex = list(range(m))
for i in range(m):
alpha = 4/(1.0+j+i)+0.01 # 步长为动态变化
rand = int(np.random.uniform(0, len(dataIndex)))
choseIndex = dataIndex[rand]
h = sigmoid(np.sum(dataArr[choseIndex]*weights))
error = h-labelList[choseIndex]
weights = weights-alpha*error*dataArr[choseIndex]
del(dataIndex[rand])
return weights
相对第一个的随机梯度下降算法,上述算法一方面使alpha在每次迭代的时候都会进行调整,缓解了数据波动或者高频波动。虽然alpha会随着迭代次数不断减小,但是不会减少至0,而这样子做也保证了在多次迭代之后新数据仍然具有一定的影响。另一方面,通过随机选取样本来更新回归系数,从而减少周期性的波动。同时每次随机从列表中选出一个值之后,就从列表中删除该值。
上述算法对机器学习实战中的使用的算法有进行更改,作者个人认为书上在选出一个值并删除该值之后的实现方式不对。按照书上实现的算法,比如我们dataIndex序列是[0, 1, 2, 3, 4, 5],那么假如第一次随机产生的数是0,那么使用的则是dataArr[0]这个样本,之后使用del删除了0下标这个数。那么dataIndex序列变成了[1, 2, 3, 4, 5]。之后再随机产生一个数,比如还是0,那么使用的还是dataArr[0]这个样本,之后del删除了下标为0的这个数,那么序列就变成了[2, 3, 4, 5]了。那么这种实现方式跟书的作者意思就有点不一致了。书中作者的意思应该是从[0, 1, 2, 3, 4, 5]这个列表中随机选取一个值,比如选取了0,然后使用第0个样本,之后删除0这个数,那么下次就变成从[1, 2, 3, 4, 5]这个列表中随机选取一个值,比如选取了1,然后使用第1个样本,之后删掉1这个数。
运行如下代码,得到的效果如下图所示,该效果达到了跟批量梯度下降算法差不多的效果,但是所使用的计算机更少
dataList, labelList = loadDataSet("testSet.txt")
weights = stocGradAscent1(dataList, labelList)
plotBestFit(dataList, labelList, weights)
![](https://i-blog.csdnimg.cn/blog_migrate/1e05246eca372a6fd289d78e8f1900cd.png)
使用更改之后的算法,来显示每次迭代时各个回归系数的变化情况,实现的代码如下
'''
记录改进的随机梯度下降过程中的历史权重值
'''
def stocGradAscent1WeightHistory(dataList, labelList, numIter=150):
dataArr = np.array(dataList)
m, n = np.shape(dataArr)
weights = np.ones(n)
weightsHistory = np.zeros((numIter*m, n))
for j in range(numIter):
dataIndex = list(range(m))
for i in range(m):
alpha = 4 / (1.0 + j + i) + 0.01 # 步长为动态变化
rand = int(np.random.uniform(0, len(dataIndex)))
choseIndex = dataIndex[rand]
h = sigmoid(np.sum(dataArr[choseIndex] * weights))
error = h - labelList[choseIndex]
weights = weights - alpha * error * dataArr[choseIndex]
weightsHistory[j*m+i, :] = weights
del (dataIndex[rand])
return weightsHistory
运行如下代码,得到的效果图如图所示。从图中可以看到没有出现周期性的波动,这归功于样本随机选择机制,同时收敛的速度也更快。
dataList, labelList = loadDataSet("testSet.txt")
weightsHistory = stocGradAscent1WeightHistory(dataList, labelList, 40)
plotSDGError(weightsHistory)
![](https://i-blog.csdnimg.cn/blog_migrate/a974d7a16f3f9ebd56aaf4c04e7b7c7d.png)
附:数据集及代码
本文所有数据集及代码,均可在https://github.com/DawnGuoDev/MachineLearningStudy/tree/master/LogisticRegression中获取。
本文参考与作者推荐
机器学习,周志华
机器学习实战,[美]Peter Harrington
CSDN_梯度下降法的向量化推导
CSDN_梯度上升算法与随机梯度上升算法的实现
不甘于「本该如此」,「多选参数 」值得关注
「多选参数 」的参数们正等着您加入,然后一起学习!