Logistic回归实战篇之预测病马死亡率（一）
作 者：崔家华
编 辑：李文臣
Python版本： Python3.x
运行平台 ： Windows
IDE ： Sublime text3
一、前言
本系列文章对梯度上升算法和改进的随机梯度上升算法进行了对比，总结了各自的优缺点，并对sklearn.linear_model.LogisticRegression进行了详细介绍。
二、改进的随机梯度上升算法
梯度上升算法在每次更新回归系数(最优参数)时，都需要遍历整个数据集。可以看一下我们之前写的梯度上升算法：
def gradAscent(dataMatIn, classLabels):
dataMatrix = np.mat(dataMatIn)
#转换成numpy的mat
labelMat = np.mat(classLabels).transpose()
#转换成numpy的mat,并进行转置
m, n = np.shape(dataMatrix)
#返回dataMatrix的大小。m为行数,n为列数。
alpha = 0.01
#移动步长,也就是学习速率,控制更新的幅度。
maxCycles = 500 #最大迭代次数
weights = np.ones((n,1))
for k in range(maxCycles):
h = sigmoid(dataMatrix * weights)
#梯度上升矢量化公式
error = labelMat - h
weights = weights + alpha *
dataMatrix.transpose() * error
return weights.getA(),weights_array
#将矩阵转换为数组，返回权重数组假设，我们使用的数据集一共有100个样本。那么，dataMatrix就是一个100*3的矩阵。每次计算h的时候，都要计算dataMatrix*weights这个矩阵乘法运算，要进行100*3次乘法运算和100*2次加法运算。同理，更新回归系数(最优参数)weights时，也需要用到整个数据集，要进行矩阵乘法运算。总而言之，该方法处理100个左右的数据集时尚可，但如果有数十亿样本和成千上万的特征，那么该方法的计算复杂度就太高了。因此，需要对算法进行改进，我们每次更新回归系数(最优参数)的时候，能不能不用所有样本呢？一次只用一个样本点去更新回归系数(最优参数)？这样就可以有效减少计算量了，这种方法就叫做随机梯度上升算法。
1、随机梯度上升算法
让我们直接看代码：
def stocGradAscent1(dataMatrix, classLabels, numIter=150):
m,n =np.shape(dataMatrix)
#返回dataMatrix的大小,m为行数,n为列数。
weights = np.ones(n) #参数初始化
for j in range(numIter):
dataIndex = list(range(m)
for i in range(m):
#降低alpha的大小，每次减小1/(j+i)。
randIndex = int(random.uniform(0,len(dataIndex)))
#随机选取样本
h = sigmoid(sum(dataMatrix[randIndex]*weights))
#选择随机选取的一个样本，计算h
error = classLabels[randIndex] - h
#计算误差
weights = weights + alpha
* error * dataMatrix[randIndex]
#更新回归系数
del(dataIndex[randIndex])
#删除已经使用的样本
return weights #返回该算法第一个改进之处在于，alpha在每次迭代的时候都会调整，并且，虽然alpha会随着迭代次数不断减小，但永远不会减小到0，因为这里还存在一个常数项。必须这样做的原因是为了保证在多次迭代之后新数据仍然具有一定的影响。如果需要处理的问题是动态变化的，那么可以适当加大上述常数项，来确保新的值获得更大的回归系数。另一点值得注意的是，在降低alpha的函数中，alpha每次减少1/(j+i)，其中j是迭代次数，i是样本点的下标。第二个改进的地方在于跟新回归系数(最优参数)时，只使用一个样本点，并且选择的样本点是随机的，每次迭代不使用已经用过的样本点。这样的方法，就有效地减少了计算量，并保证了回归效果。
编写代码如下，看下改进的随机梯度上升算法分类效果如何：
# -*- coding:UTF-8 -*-
from matplotlib.font_manager import FontProperties
import matplotlib.pyplot as plt
import numpy as np
import random
"""函数说明:加载数据
Parameters:
无
Returns:
dataMat - 数据列表
labelMat - 标签列表
Author:
Jack Cui
Blog:
http://blog.csdn.net/c406495762
Zhihu:
https://www.zhihu.com/people/Jack--Cui/
Modify:
2017-08-28
"""
def loadDataSet():
dataMat = [] #创建数据列表
labelMat = [] #创建标签列表
fr = open('testSet.txt') #打开文件
for line in fr.readlines(): #逐行读取
lineArr = line.strip().split()
#去回车，放入列表
dataMat.append([1.0, float(lineArr[0]),
float(lineArr[1])]) #添加数据
labelMat.append(int(lineArr[2])) #添加标签
fr.close() #关闭文件 return dataMat, labelMat #返回
"""
函数说明:sigmoid函数
Parameters:
inX - 数据
Returns:
sigmoid函数
Author:
Jack Cui
Blog:
http://blog.csdn.net/c406495762
Zhihu:
https://www.zhihu.com/people/Jack--Cui/
Modify:
2017-08-28
"""
def sigmoid(inX):
return 1.0 / (1 + np.exp(-inX))
"""
函数说明:绘制数据集
Parameters:
weights - 权重参数数组
Returns:
无
Author: Jack Cui
Blog:
http://blog.csdn.net/c406495762
Zhihu:
https://www.zhihu.com/people/Jack--Cui/
Modify:
2017-08-30
"""
def plotBestFit(weights):
dataMat, labelMat = loadDataSet() #加载数据集
dataArr = np.array(dataMat)
#转换成numpy的array数组
n = np.shape(dataMat)[0] #数据个数
xcord1 = []; ycord1 = [] #正样本
xcord2 = []; ycord2 = [] #负样本
for i in range(n): #根据数据集标签进行分类
if int(labelMat[i]) == 1:
xcord1.append(dataArr[i,1])
ycord1.append(dataArr[i,2])
#1为正样本
else:
xcord2.append(dataArr[i,1])
ycord2.append(dataArr[i,2])
#0为负样本
fig = plt.figure()
ax = fig.add_subplot(111) #添加subplot
ax.scatter(xcord1, ycord1, s = 20, c = 'red',
marker = 's',alpha=.5) #绘制正样本
ax.scatter(xcord2, ycord2, s = 20,
c = 'green',alpha=.5) #绘制负样本
x = np.arange(-3.0, 3.0, 0.1)
y = (-weights[0] - weights[1] * x) / weights[2]
ax.plot(x, y)
plt.title('BestFit') #绘制title plt.xlabel('X1');
plt.ylabel('X2') #绘制label
plt.show()
"""
函数说明:改进的随机梯度上升算法
Parameters:
dataMatrix - 数据数组
classLabels - 数据标签
numIter - 迭代次数
Returns: weights - 求得的回归系数数组(最优参数)
Author: Jack Cui
Blog: http://blog.csdn.net/c406495762
Zhihu: https://www.zhihu.com/people/Jack--Cui/
Modify:
2017-08-31
"""
def stocGradAscent1(dataMatrix, classLabels, numIter=150):
m,n= np.shape(dataMatrix)
#返回dataMatrix的大小m为行数,n为列数。
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
#降低alpha的大小，每次减小1/(j+i)。
randIndex = int(random.uniform(0,len(dataIndex)))
#随机选取样本
h = sigmoid(sum(dataMatrix[randIndex]*weights))
#选择随机选取的一个样本，计算h
error = classLabels[randIndex] - h #计算误差
weights=weights+alpha*error*dataMatrix[randIndex]
#更新回归系数
del(dataIndex[randIndex])#删除已经使用的样本 return weights #返回
if __name__ == '__main__':
dataMat, labelMat = loadDataSet()
weights = stocGradAscent1(np.array(dataMat), labelMat)
plotBestFit(weights)代码运行结果：
2、回归系数与迭代次数的关系
可以看到分类效果也是不错的。不过，从这个分类结果中，我们不好看出迭代次数和回归系数的关系，也就不能直观的看到每个回归方法的收敛情况。因此，我们编写程序，绘制出回归系数和迭代次数的关系曲线：
# -*- coding:UTF-8 -*-
from matplotlib.font_manager
import FontProperties
import matplotlib.pyplot as plt
import numpy as np
import random
"""
函数说明:加载数据
Parameters:
无
Returns:
dataMat - 数据列表
labelMat - 标签列表
Author: Jack Cui
Blog: http://blog.csdn.net/c406495762
Zhihu: https://www.zhihu.com/people/Jack--Cui/
Modify:
2017-08-28
"""
def loadDataSet():
dataMat = [] #创建数据列表
labelMat = [] #创建标签列表
fr = open('testSet.txt') #打开文件
for line in fr.readlines(): #逐行读取
lineArr = line.strip().split()#去回车，放入列表
dataMat.append([1.0, float(lineArr[0]),
float(lineArr[1])]) #添加数据
labelMat.append(int(lineArr[2])) #添加标签
fr.close() #关闭文件 return dataMat, labelMat #返回
"""
函数说明:sigmoid函数
Parameters:
inX - 数据
Returns:
sigmoid函数
Author: Jack Cui
Blog: http://blog.csdn.net/c406495762
Zhihu: https://www.zhihu.com/people/Jack--Cui/
Modify: 2017-08-28
"""
def sigmoid(inX):
return 1.0 / (1 + np.exp(-inX))
"""
函数说明:梯度上升算法
Parameters:
dataMatIn - 数据集
classLabels - 数据标签
Returns:
weights.getA() - 求得的权重数组(最优参数)
weights_array - 每次更新的回归系数
Author: Jack Cui
Blog: http://blog.csdn.net/c406495762
Zhihu: https://www.zhihu.com/people/Jack--Cui/
Modify: 2017-08-28
"""
def gradAscent(dataMatIn, classLabels):
dataMatrix = np.mat(dataMatIn)#转换成numpy的mat
labelMat = np.mat(classLabels).transpose()
#转换成numpy的mat,并进行转置
m, n = np.shape(dataMatrix)
#返回dataMatrix的大小。m为行数,n为列数。
alpha = 0.01
#移动步长,也就是学习速率,控制更新的幅度。
maxCycles = 500 #最大迭代次数
weights = np.ones((n,1))
weights_array = np.array([])
for k in range(maxCycles):
h = sigmoid(dataMatrix * weights)
#梯度上升矢量化公式
error = labelMat - h
weights = weights + alpha
* dataMatrix.transpose() * error
weights_array = np.append(weights_array,weights)
weights_array = weights_array.reshape(maxCycles,n) return weights.getA(),weights_array
#将矩阵转换为数组,并返回
"""
函数说明:改进的随机梯度上升算法
Parameters:
dataMatrix - 数据数组
classLabels - 数据标签
numIter - 迭代次数
Returns:
weights - 求得的回归系数数组(最优参数)
weights_array - 每次更新的回归系数
Author: Jack Cui
Blog: http://blog.csdn.net/c406495762
Zhihu: https://www.zhihu.com/people/Jack--Cui/
Modify: 2017-08-31
"""
def stocGradAscent1(dataMatrix, classLabels, numIter=150):
m,n= np.shape(dataMatrix)
#返回dataMatrix的大小,m为行数,n为列数。
weights = np.ones(n) #参数初始化
weights_array = np.array([])
#存储每次更新的回归系数
for j in range(numIter):
dataIndex = list(range(m))
for i in range(m):
alpha = 4/(1.0+j+i)+0.01
#降低alpha的大小,每次减小1/(j+i)。
randIndex =int(random.uniform(0,len(dataIndex)))
#随机选取样本
h = sigmoid(sum(dataMatrix[randIndex]*weights))
#选择随机选取的一个样本，计算h
error = classLabels[randIndex] - h
#计算误差
weights = weights + alpha * error
* dataMatrix[randIndex] #更新回归系数
weights_array = np.append(weights_array,
weights,axis=0) #添加回归系数到数组中
del(dataIndex[randIndex])
#删除已经使用的样本
weights_array = weights_array.reshape(numIter*m,n)
#改变维度 return weights,weights_array #返回
"""
函数说明:绘制回归系数与迭代次数的关系
Parameters:
weights_array1 - 回归系数数组1
weights_array2 - 回归系数数组2
Returns:
无
Author: Jack Cui
Blog: http://blog.csdn.net/c406495762
Zhihu: https://www.zhihu.com/people/Jack--Cui/
Modify:
2017-08-30
"""
def plotWeights(weights_array1,weights_array2):
#设置汉字格式
font =
FontProperties(fname=r"c:\windows\fonts\simsun.ttc",size=14)
#将fig画布分隔成1行1列,不共享x轴和y轴
#fig画布的大小为（13,8）当nrow=3,nclos=2时
#代表fig画布被分为个区域
#axs[0][0]表示第一行第一列
fig, axs =
plt.subplots(nrows=3, ncols=2,sharex=False,
sharey=False, figsize=(20,10))
x1 = np.arange(0, len(weights_array1), 1)
#绘制w0与迭代次数的关系
axs[0][0].plot(x1,weights_array1[:,0])
axs0_title_text =
axs[0][0].set_title(u'梯度上升算法：
回归系数与迭代次数关系',FontProperties=font)
axs0_ylabel_text = axs[0][0].
set_ylabel(u'W0', FontProperties=font)
plt.setp(axs0_title_text,
size=20,weight='bold',color='black')
plt.setp(axs0_ylabel_text,
size=20,weight='bold',color='black')
#绘制w1与迭代次数的关系
axs[1][0].plot(x1,weights_array1[:,1])
axs1_ylabel_text = axs[1][0].
set_ylabel(u'W1',FontProperties=font)
plt.setp(axs1_ylabel_text,
size=20,weight='bold',color='black')
#绘制w2与迭代次数的关系
axs[2][0].plot(x1,weights_array1[:,2])
axs2_xlabel_text = axs[2][0].
set_xlabel(u'迭代次数',FontProperties=font)
axs2_ylabel_text = axs[2][0].
set_ylabel(u'W1', FontProperties=font)
plt.setp(axs2_xlabel_text,
size=20,weight='bold',color='black')
plt.setp(axs2_ylabel_text,
size=20,weight='bold',color='black')
x2 = np.arange(0, len(weights_array2), 1)
#绘制w0与迭代次数的关系
axs[0][1].plot(x2,weights_array2[:,0])
axs0_title_text =
axs[0][1].set_title(u'改进的随机梯度上升算法：
回归系数与迭代次数关系',FontProperties=font)
axs0_ylabel_text = axs[0][1].
set_ylabel(u'W0',FontProperties=font)
plt.setp(axs0_title_text,
size=20, weight='bold',color='black')
plt.setp(axs0_ylabel_text,
size=20,weight='bold',color='black')
#绘制w1与迭代次数的关系
axs[1][1].plot(x2,weights_array2[:,1])
axs1_ylabel_text = axs[1][1].
set_ylabel(u'W1',FontProperties=font)
plt.setp(axs1_ylabel_text,
size=20,weight='bold',color='black')
#绘制w2与迭代次数的关系
axs[2][1].plot(x2,weights_array2[:,2])
axs2_xlabel_text = axs[2][1].
set_xlabel(u'迭代次数', FontProperties=font)
axs2_ylabel_text = axs[2][1].
set_ylabel(u'W1', FontProperties=font)
plt.setp(axs2_xlabel_text,
size=20,weight='bold',color='black')
plt.setp(axs2_ylabel_text,
size=20,weight='bold',color='black')
plt.show()
if __name__ == '__main__':
dataMat, labelMat = loadDataSet()
weights1,weights_array1 =
stocGradAscent1(np.array(dataMat), labelMat)
weights2,weights_array2 =
gradAscent(dataMat, labelMat)
plotWeights(weights_array1, weights_array2)运行结果如下：
由于改进的随机梯度上升算法，随机选取样本点，所以每次的运行结果是不同的。但是大体趋势是一样的。我们改进的随机梯度上升算法收敛效果更好。为什么这么说呢？让我们分析一下。我们一共有100个样本点，改进的随机梯度上升算法迭代次数为150。而上图显示15000次迭代次数的原因是，使用一次样本就更新一下回归系数。因此，迭代150次，相当于更新回归系数150*100=15000次。简而言之，迭代150次，更新1.5万次回归参数。从上图左侧的改进随机梯度上升算法回归效果中可以看出，其实在更新2000次回归系数的时候，已经收敛了。相当于遍历整个数据集20次的时候，回归系数已收敛。训练已完成。
再让我们看看上图右侧的梯度上升算法回归效果，梯度上升算法每次更新回归系数都要遍历整个数据集。从图中可以看出，当迭代次数为300多次的时候，回归系数才收敛。凑个整，就当它在遍历整个数据集300次的时候已经收敛好了。
没有对比就没有伤害，改进的随机梯度上升算法，在遍历数据集的第20次开始收敛。而梯度上升算法，在遍历数据集的第300次才开始收敛。想像一下，大量数据的情况下，谁更牛逼？
本系列篇章：
Logistic回归实战篇之预测病马死亡率(一)
Logistic回归实战篇之预测病马死亡率(二)
Logistic回归实战篇之预测病马死亡率(三)
往期精彩文章推荐
扫描燕哥微信号，
拉你进机器学习大牛群。
福利满满，名额已不多…
群里目前包括：
清华张长水教授，
清华顾险峰教授，
北大黄铁军教授，
西安电子科技大学焦李成教授，
新加坡南洋理工大学黄广斌教授，
北交李清勇教授
等等……