《Machine Learning in action》,机器学习实战(笔记)之Logistic regression
使用工具
- Python3.7
- pycharm
- anaconda
- jupyter notebook
一、前言
- 在上一篇中,我们介绍了梯度上升算法,,而在这一篇中。我们将对前面的算法进行改进,随机梯度上升算法,我们同时将会进行,对比,分析其中的优缺点,并对sklearn.linear_model.LogisticRegression进行了详细介绍。
二、改进的梯度上升算法-随机梯度上升算法
- 梯度上升算法在每次更新回归系数(最优参数)时,都需要遍历整个数据集。可以看一下我们之前写的梯度上升算法:
def gradAscent(dataMatIn, classLabels):
# 将数据集转换为numpy的mat
dataMatrix = np.mat(dataMatIn)
# 将标签转换成numpy的mat。并且进行转置
# transpose 就是表示转置
labelMat = np.mat(classLabels).transpose()
# 返回dataMatrix的大小。m为行数,n为列数。
m, n = np.shape(dataMatrix)
# 移动步长,也就是学习速率,控制更新的幅度。
alpha = 0.001
# 最大迭代次数
maxCycles = 500
# 权重的初始化为0
weights = np.ones((n,1))
for k in range(maxCycles):
# h就是代表z的含义,梯度上升矢量化公式
h = sigmoid(dataMatrix * weights)
error = labelMat - h
weights = weights + alpha * dataMatrix.transpose() * error
# 将矩阵转换为数组,返回权重数组
# 这里的最优权重数组就是所求的最佳的解
return weights.getA()
- 我们使用的数据集一共有100个样本。那么,dataMatrix就是一个1003的矩阵。每次计算h的时候,都要计算dataMatrixweights这个矩阵乘法运算,要进行1003次乘法运算和1002次加法运算。同理,更新回归系数(最优参数)weights时,也需要用到整个数据集,要进行矩阵乘法运算。总而言之,使用这种方法时,会让我们的计算量非常大,当数据非常大时,那么该方法的计算复杂度就太高了。因此,需要对算法进行改进,我们每次更新回归系数(最优参数)的时候,能不能不用所有样本呢?一次只用一个样本点去更新回归系数(最优参数)?这样就可以有效减少计算量了,这种方法就叫做随机梯度上升算法 。
1. 随机梯度上升算法
- 代码如下:
def stocGradAscent1(dataMatrix, classLabels, numIter=150):
# 返回dataMatrix的大小。m为行数,n为列数。
m,n = np.shape(dataMatrix)
# 参数初始化
weights = np.ones(n)
for j in range(numIter):
dataIndex = list(range(m))
for i in range(m):
# 降低alpha的大小,每次减小1/(j+i)。
# 相比如前面所做的改变:Alpha在每一次迭代中,都进行了改变
alpha = 4/(1.0+j+i)+0.01
# 随机选取样本
# 这里也是做出的改变:是随机的选择样本,而不是遍历所有的样本
randIndex = int(random.uniform(0,len(dataIndex))
# 随机选取的一个样本,计算h
# h就是代表z的含义,梯度上升矢量化公式
h = sigmoid(sum(dataMatrix[randIndex]*weights))
# 计算误差
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是样本点的下标。第二个改进的地方在于跟新回归系数(最优参数)时,只使用一个样本点,并且选择的样本点是随机的,每次迭代不使用已经用过的样本点。这样的方法,就有效地减少了计算量,并保证了回归效果。
2. 利用随机梯度上升算法分类
- 直接是上代码:
# -*- coding:UTF-8 -*-
import matplotlib.pyplot as plt
import numpy as np
import random
def loadDataSet():
# 创建数据列表,用来放入我们的数据
dataMat = []
# 该列表就是用来存放我们的标签的
labelMat = []
# 读取我们的数据集txt文件
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
# 定义sigmiod函数
def sigmoid(inx):
return 1.0/(1 + np.exp(-inx))
'''
下面就是我们的随机梯度上升算法梯度上升算法的实现
stocGradAscent
'''
def stocGradAscent1(dataMatrix, classLabels, numIter=150):
# 返回dataMatrix的大小。m为行数,n为列数。
m,n = np.shape(dataMatrix)
# 参数初始化
weights = np.ones(n)
for j in range(numIter):
dataIndex = list(range(m))
for i in range(m):
# 降低alpha的大小,每次减小1/(j+i)。
# 相比如前面所做的改变:Alpha在每一次迭代中,都进行了改变
alpha = 4/(1.0+j+i)+0.01
# 随机选取样本
# 这里也是做出的改变:是随机的选择样本,而不是遍历所有的样本
randIndex = int(random.uniform(0,len(dataIndex)))
# 随机选取的一个样本,计算h
# h就是代表z的含义,梯度上升矢量化公式
h = sigmoid(sum(dataMatrix[randIndex]*weights))
# 计算误差
error = classLabels[randIndex] - h
# 更新回归系数
weights = weights + alpha * error * dataMatrix[randIndex]
# 将已经使用过得样本进行删除
del(dataIndex[randIndex])
# 返回权重
return weights
def plotBestFit(weights):
# 数据集的实例化,也就是加载数据
dataMat, labelMat = loadDataSet()
# 将数据转换为有序列的numpy中的array数组
dataArr = np.array(dataMat)
# 获取数据的个数,也就是行数
n = np.shape(dataArr)[0]
# 下面两个不同的列表就是分别计算正样本和负样本的
# 正样本
xcord1 = []; ycord1 = []
# 负样本
xcord2 = []; ycord2 = []
# 根据数据集的标签进行分类
for i in range(n):
# 如果标签是1
if int(labelMat[i]) == 1:
# 将为1 的放入我们的正样本中
xcord1.append(dataArr[i,1]); ycord1.append(dataArr[i,2])
else:
# 其他,既是为0时放入负样本中
xcord2.append(dataArr[i,1]); ycord2.append(dataArr[i,2])
# 绘制绘图框
fig = plt.figure()
# 添加子图subplot
ax = fig.add_subplot(111)
# 实现正样本的可视化
ax.scatter(xcord1, ycord1, s = 20, c = 'red', marker = 's')
# 实现负样本的可视化
ax.scatter(xcord2, ycord2, s = 20, c = 'blue')
x = np.arange(-3.0, 3.0, 0.1)
y = (-weights[0]- weights[1] *x)/ weights[2]
ax.plot(x,y)
# 设置绘图框的title
plt.title('BestFit')
plt.xlabel('x'); plt.ylabel('y')
plt.show()
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
def loadDataSet():
# 创建数据列表,用来放入我们的数据
dataMat = []
# 该列表就是用来存放我们的标签的
labelMat = []
# 读取我们的数据集txt文件
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
def sigmoid(inX):
return 1.0 / (1 + np.exp(-inX))
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 #将矩阵转换为数组,并返回
def stocGradAscent1(dataMatrix, classLabels, numIter=150):
# 返回dataMatrix的大小。m为行数,n为列数。
m,n = np.shape(dataMatrix)
# 参数初始化
weights = np.ones(n)
# 存储每次更新的回归系数
weights_array = np.array([])
for j in range(numIter):
dataIndex = list(range(m))
for i in range(m):
# 降低alpha的大小,每次减小1/(j+i)。
# 相比如前面所做的改变:Alpha在每一次迭代中,都进行了改变
alpha = 4/(1.0+j+i)+0.01
# 随机选取样本
# 这里也是做出的改变:是随机的选择样本,而不是遍历所有的样本
randIndex = int(random.uniform(0,len(dataIndex)))
# 随机选取的一个样本,计算h
# h就是代表z的含义,梯度上升矢量化公式
h = sigmoid(sum(dataMatrix[randIndex]*weights))
# 计算误差
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
"""
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'W2',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'W2',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)
- 运行的结果如下:
- 分析:
- 改进后的随机梯度上升,随机的选取样本点,所以造成每一次的结果都会有一定的差别的,但是大体上面都是一样的,通过观察会发现我们的改进后的算法,的收敛效果会更好一些,所以我们一般情况下都是会采取这样的改进后的随机梯度上升的算法。
三、estimating hores fatalities from colic-(根据疝气病症状预测马的死亡率)
1. 背景介绍:
-
本次实战内容,将使用Logistic回归来预测患疝气病的马的存活问题。原始数据集下载地址:http://archive.ics.uci.edu/ml/datasets/Horse+Colic
-
这里的数据包含了368个样本和28个特征。这种病不一定源自马的肠胃问题,其他问题也可能引发马疝病。该数据集中包含了医院检测马疝病的一些指标,有的指标比较主观,有的指标难以测量,例如马的疼痛级别。另外需要说明的是,除了部分指标主观和难以测量外,该数据还存在一个问题,数据集中有30%的值是缺失的。下面将首先介绍如何处理数据集中的数据缺失问题,然后再利用Logistic回归和随机梯度上升算法来预测病马的生死。
2. 准备数据
数据中的缺失值是一个非常棘手的问题,很多文献都致力于解决这个问题。那么,数据缺失究竟带来了什么问题?假设有100个样本和20个特征,这些数据都是机器收集回来的。若机器上的某个传感器损坏导致一个特征无效时该怎么办?它们是否还可用?答案是肯定的。因为有时候数据相当昂贵,扔掉和重新获取都是不可取的,所以必须采用一些方法来解决这个问题。下面给出了一些可选的做法:
- 使用可用特征的均值来填补缺失值;
- 使用特殊值来填补缺失值,如-1;
- 忽略有缺失值的样本;
- 使用相似样本的均值添补缺失值;
- 使用另外的机器学习算法预测缺失值。
数据的预处理将会做一下的两件事:
- 1.如果测试集中一条数据的特征值已经确实,那么我们选择实数0来替换所有缺失值,因为本文使用Logistic回归。因此这样做不会影响回归系数的值。sigmoid(0)=0.5,即它对结果的预测不具有任何倾向性。
- 如果测试集中一条数据的类别标签已经缺失,那么我们将该类别数据丢弃,因为类别标签与特征不同,很难确定采用某个合适的值来替换。
原始的数据集经过处理,保存为两个文件:horseColicTest.txt和horseColicTraining.txt。已经处理好的“干净”可用的数据集下载地址: - https://github.com/Tupeng2019/Machine_jichu/blob/master/Machine-Learning/Machine-Learning_note/Logistic regression/2_ horse fatalities from colic/horseColicTest.txt
- https://github.com/Tupeng2019/Machine_jichu/blob/master/Machine-Learning/Machine-Learning_note/Logistic regression/2_ horse fatalities from colic/horseColicTraining.txt
有了这些数据集,我们只需要一个Logistic分类器,就可以利用该分类器来预测病马的生死问题了。
3.构建分类器-Logistic regression classification function
我们先用我们自己写的随机梯度上升算法来实现我们的逻辑回归
# -*- coding:UTF-8 -*-
import numpy as np
import random
def sigmoid(inX):
return 1.0 / (1 + np.exp(-inX))
'''
下面就是我们的随机梯度上升算法梯度上升算法的实现
stocGradAscent
'''
def stocGradAscent1(dataMatrix, classLabels, numIter=150):
# 返回dataMatrix的大小。m为行数,n为列数。
m,n = np.shape(dataMatrix)
# 参数初始化
weights = np.ones(n)
for j in range(numIter):
dataIndex = list(range(m))
for i in range(m):
# 降低alpha的大小,每次减小1/(j+i)。
# 相比如前面所做的改变:Alpha在每一次迭代中,都进行了改变
alpha = 4/(1.0+j+i)+0.01
# 随机选取样本
# 这里也是做出的改变:是随机的选择样本,而不是遍历所有的样本
randIndex = int(random.uniform(0,len(dataIndex)))
# 随机选取的一个样本,计算h
# h就是代表z的含义,梯度上升矢量化公式
h = sigmoid(sum(dataMatrix[randIndex]*weights))
# 计算误差
error = classLabels[randIndex] - h
# 更新回归系数
weights = weights + alpha * error * dataMatrix[randIndex]
# 将已经使用过得样本进行删除
del(dataIndex[randIndex])
# 返回权重
return weights
'''
这就是:分类函数
inX - 特征向量
weights - 回归系数(就是w,相当于权重)
'''
def classifyVector(inX, weights):
# 计算Sigmoid函数的值
# 当值是>0.5时,就是分为1类
# 当值是<0.5时,就是分为0类
prob = sigmoid(sum(inX*weights))
if prob > 0.5:
return 1.0
else:
return 0.0
'''
这就是我们的colic的测试函数
'''
def colicTest():
# 这就是读取训练集
frTrain = open('horseColicTraining.txt')
# 读取测试集
frTest = open('horseColicTest.txt')
trainingSet = []; trainingLabels = []
for line in frTrain.readlines():
currLine = line.strip().split('\t')
lineArr = []
for i in range(len(currLine)-1):
lineArr.append(float(currLine[i]))
trainingSet.append(lineArr)
trainingLabels.append(float(currLine[-1]))
# 使用改进后的随机梯度上升算法
trainWeights = stocGradAscent1(np.array(trainingSet), trainingLabels, 500)
errorCount = 0; numTestVec = 0.0
for line in frTest.readlines():
numTestVec += 1.0
currLine = line.strip().split('\t')
lineArr =[]
for i in range(len(currLine)-1):
lineArr.append(float(currLine[i]))
if int(classifyVector(np.array(lineArr), trainWeights))!= int(currLine[-1]):
errorCount += 1
# 计算错误率
errorRate = (float(errorCount)/numTestVec) * 100
print("测试集错误率为: %.2f%%" % errorRate)
if __name__ == '__main__':
colicTest()
运行结果:
从上面可以看出,我们正确的分辨率为35.82%,本身就有30%的数据缺失,所以来说这个结果还是不错的,但是我们可以尝试一下用直接梯度上升算法来试试,看看是否有所不同的结果。
# -*- coding:UTF-8 -*-
import numpy as np
import random
def sigmoid(inX):
return 1.0 / (1 + np.exp(-inX))
def gradAscent(dataMatIn, classLabels):
dataMatrix = np.mat(dataMatIn)
labelMat = np.mat(classLabels).transpose()
m, n = np.shape(dataMatrix)
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()
"""
函数说明:使用Python写的Logistic分类器做预测
"""
def colicTest():
frTrain = open('horseColicTraining.txt')
frTest = open('horseColicTest.txt')
trainingSet = []; trainingLabels = []
for line in frTrain.readlines():
currLine = line.strip().split('\t')
lineArr = []
for i in range(len(currLine)-1):
lineArr.append(float(currLine[i]))
trainingSet.append(lineArr)
trainingLabels.append(float(currLine[-1]))
trainWeights = gradAscent(np.array(trainingSet), trainingLabels)
errorCount = 0; numTestVec = 0.0
for line in frTest.readlines():
numTestVec += 1.0
currLine = line.strip().split('\t')
lineArr =[]
for i in range(len(currLine)-1):
lineArr.append(float(currLine[i]))
if int(classifyVector(np.array(lineArr), trainWeights[:,0]))!= int(currLine[-1]):
errorCount += 1
errorRate = (float(errorCount)/numTestVec) * 100
print("测试集错误率为: %.2f%%" % errorRate)
"""
函数说明:分类函数
Parameters:
inX - 特征向量
weights - 回归系数
"""
def classifyVector(inX, weights):
prob = sigmoid(sum(inX*weights))
if prob > 0.5:
return 1.0
else:
return 0.0
if __name__ == '__main__':
colicTest()
运行结果如下:
从上面的结果可以看出,实验的分类结果还是较好于使用随机梯度上升的结果。
所以我们可以得出这样的总结:
- 当数据集较小时,我们使用梯度上升算法
- 当数据集较大时,我们使用改进的随机梯度上升算法
对应的,在Sklearn中(将在下一篇当中 介绍),我们就可以根据数据情况选择优化算法,比如数据较小的时候,我们使用liblinear,数据较大时,我们使用sag和saga。