机器学习之旅(四):Logistic回归
Logistic回归工作原理
logistic回归是一种广义线性回归(generalized linear model),因此与多重线性回归分析有很多相同之处。它们的模型形式基本上相同,都具有 w‘x+b,其中w和b是待求参数,其区别在于他们的因变量不同,多重线性回归直接将w‘x+b作为因变量,即y =w‘x+b,而logistic回归则通过函数L将w‘x+b对应一个隐状态p,p =L(w‘x+b),然后根据p 与1-p的大小决定因变量的值。如果L是logistic函数,就是logistic回归,如果L是多项式函数就是多项式回归。
logistic回归的因变量可以是二分类的,也可以是多分类的,但是二分类的更为常用,也更加容易解释,多类可以使用softmax方法进行处理。实际中最为常用的就是二分类的logistic回归。
Logistic回归的优缺点
- 优点:计算代价不高,易于理解和实现。
- 缺点:容易欠拟合,分类精度可能不高。
- 使用数据范围:数值型和标称型数据。
Logistic回归的一般流程
(1)收集数据:采用任意方法收集数据。
(2)准备数据:由于需要进行距离计算,因此要求数据类型为数值型。另外,结构化数据格式则最佳。
(3)分析数据:采用任意方法对数据进行分析。
(4)训练算法:大部分时间将用于训练,训练的目的是为了找到最佳的分类回归系数。
(5)测试算法:一旦训练步骤完成,分类将会很快。
(6)使用算法:首先,我们需要输入一些数据,并将其转换成对应的结构化数值;接着,基于训练好的回归系数就可以对这些数值进行简单的回归计算,半丁他们属于哪个类别;在这之后,我们就可以在输出的类别上做一些其他分析工作。
编写代码实现Logistic回归算法
首先导入一个数据集,看下数据集的分布情况
import numpy as np
import matplotlib.pyplot as plt
"""
函数说明:加载数据
Parameters:
无
Returns:
dataMat - 数据列表
labelMat - 标签列表
Author:
Jack Cui
"""
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 # 返回
"""
函数说明:绘制数据集
Parameters:
无
Returns:
无
Author:
Jack Cui
"""
def plotDataSet():
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]) # 1为正样本
ycord1.append(dataArr[i, 2])
else:
xcord2.append(dataArr[i, 1]) # 0为负样本
ycord2.append(dataArr[i, 2])
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) # 绘制负样本
plt.title('DataSet') # 绘制title
plt.xlabel('x') # 绘制label
plt.ylabel('y')
plt.show()
if __name__ == '__main__':
plotDataSet()
运行结果如下
从上图可以看出数据的分布情况。假设Sigmoid函数的输入记为z,那么z=w0x0 + w1x1 + w2x2,即可将数据分割开。其中,x0为全是1的向量,x1为数据集的第一列数据,x2为数据集的第二列数据。另z=0,则0=w0 + w1x1 + w2x2。横坐标为x1,纵坐标为x2。这个方程未知的参数为w0,w1,w2,也就是我们需要求的回归系数(最优参数)。。
编写代码训练算法
根据矢量化的公式,编写代码如下
import numpy as np
"""
函数说明:加载数据
Parameters:
无
Returns:
dataMat - 数据列表
labelMat - 标签列表
Author:
Jack Cui
"""
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
"""
def sigmoid(inX):
return 1.0 / (1 + np.exp(-inX))
"""
函数说明:梯度上升算法
Parameters:
dataMatIn - 数据集
classLabels - 数据标签
Returns:
weights.getA() - 求得的权重数组(最优参数)
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.001 # 移动步长,也就是学习速率,控制更新的幅度
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() # 将矩阵转换为数组,返回权重数组
if __name__ == '__main__':
dataMat, labelMat = loadDataSet()
print(gradAscent(dataMat, labelMat))
运行结果如下图所示
可以看出,我们已经求解出回归系数[w0,w1,w2]。 通过求解出的参数,我们就可以确定不同类别数据之间的分隔线,画出决策边界。
编写代码绘制决策边界
我们已经解出了一组回归系数,它确定了不同类别数据之间的分隔线。现在开始绘制这个分隔线,编写代码如下
#!/usr/bin/env python
# -*- coding:utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
"""
函数说明:加载数据
Parameters:
无
Returns:
dataMat - 数据列表
labelMat - 标签列表
Author:
Jack Cui
"""
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
"""
def sigmoid(inX):
return 1.0 / (1 + np.exp(-inX))
"""
函数说明:梯度上升算法
Parameters:
dataMatIn - 数据集
classLabels - 数据标签
Returns:
weights.getA() - 求得的权重数组(最优参数)
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.001 # 移动步长,也就是学习速率,控制更新的幅度
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() # 将矩阵转换为数组,返回权重数组
"""
函数说明:绘制数据集
Parameters:
无
Returns:
无
Author:
Jack Cui
"""
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]) # 1为正样本
ycord1.append(dataArr[i, 2])
else:
xcord2.append(dataArr[i, 1]) # 0为负样本
ycord2.append(dataArr[i, 2])
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') # 绘制label
plt.ylabel('Y1')
plt.show()
if __name__ == '__main__':
dataMat, labelMat = loadDataSet()
weights = gradAscent(dataMat, labelMat)
plotBestFit(weights)
运行结果如下
这个分类结果相当不错,从上图可以看出,只分错了几个点而已。但是,尽管例子简单切数据集很小,但是这个方法却需要大量的计算(300次乘法)。
假设,我们使用的数据集一共有100个样本。那么,dataMatrix就是一个100*3的矩阵。每次计算h的时候,都要计算dataMatrix*weights这个矩阵乘法运算,要进行100*3次乘法运算和100*2次加法运算。同理,更新回归系数(最优参数)weights时,也需要用到整个数据集,要进行矩阵乘法运算。总而言之,该方法处理100个左右的数据集时尚可,但如果有数十亿样本和成千上万的特征,那么该方法的计算复杂度就太高了。因此,需要对算法进行改进,我们每次更新回归系数(最优参数)的时候,能不能不用所有样本呢?一次只用一个样本点去更新回归系数(最优参数)?这样就可以有效减少计算量了,这种方法就叫做随机梯度上升算法。
编写代码随机梯度上升算法
该算法第一个改进之处在于,alpha在每次迭代的时候都会调整,并且,虽然alpha会随着迭代次数不断减小,但永远不会减小到0,因为这里还存在一个常数项。必须这样做的原因是为了保证在多次迭代之后新数据仍然具有一定的影响。如果需要处理的问题是动态变化的,那么可以适当加大上述常数项,来确保新的值获得更大的回归系数。另一点值得注意的是,在降低alpha的函数中,alpha每次减少1/(j+i),其中j是迭代次数,i是样本点的下标。第二个改进的地方在于跟新回归系数(最优参数)时,只使用一个样本点,并且选择的样本点是随机的,每次迭代不使用已经用过的样本点。这样的方法,就有效地减少了计算量,并保证了回归效果。
编写代码如下,看下改进的随机梯度上升算法分类效果如何
import random
import numpy as np
import matplotlib.pyplot as plt
"""
函数说明:加载数据
Parameters:
无
Returns:
dataMat - 数据列表
labelMat - 标签列表
Author:
Jack Cui
"""
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
"""
def sigmoid(inX):
return 1.0 / (1 + np.exp(-inX))
"""
函数说明:改进的随机梯度上升算法
Parameters:
dataMatrix - 数据数组
classLabels - 数据标签
numIter - 迭代次数
Returns:
weights - 求得的回归系数数组(最优参数)
Author:
Jack Cui
"""
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 # 返回
"""
函数说明:绘制数据集
Parameters:
无
Returns:
无
Author:
Jack Cui
"""
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]) # 1为正样本
ycord1.append(dataArr[i, 2])
else:
xcord2.append(dataArr[i, 1]) # 0为负样本
ycord2.append(dataArr[i, 2])
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') # 绘制label
plt.ylabel('Y1')
plt.show()
if __name__ == '__main__':
dataMat, labelMat = loadDataSet()
weights = stocGradAscent1(np.array(dataMat), labelMat)
plotBestFit(weights)
代码运行结果:
可以看到分类效果也是不错的。不过,从这个分类结果中,我们不好看出迭代次数和回归系数的关系,也就不能直观的看到每个回归方法的收敛情况。因此,我们编写程序,绘制出回归系数和迭代次数的关系曲线
import random
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.font_manager import FontProperties
"""
函数说明:加载数据
Parameters:
无
Returns:
dataMat - 数据列表
labelMat - 标签列表
Author:
Jack Cui
"""
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
"""
def sigmoid(inX):
return 1.0 / (1 + np.exp(-inX))
"""
函数说明:梯度上升算法
Parameters:
dataMatIn - 数据集
classLabels - 数据标签
Returns:
weights.getA() - 求得的权重数组(最优参数)
Author:
Jack Cui
"""
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 - 求得的回归系数数组(最优参数)
Author:
Jack Cui
"""
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
"""
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)
运行结果如下
由于改进的随机梯度上升算法,随机选取样本点,所以每次的运行结果是不同的。但是大体趋势是一样的。我们改进的随机梯度上升算法收敛效果更好。为什么这么说呢?让我们分析一下。我们一共有100个样本点,改进的随机梯度上升算法迭代次数为150。而上图显示15000次迭代次数的原因是,使用一次样本就更新一下回归系数。因此,迭代150次,相当于更新回归系数150*100=15000次。简而言之,迭代150次,更新1.5万次回归参数。从上图左侧的改进随机梯度上升算法回归效果中可以看出,其实在更新2000次回归系数的时候,已经收敛了。相当于遍历整个数据集20次的时候,回归系数已收敛。训练已完成。
再让我们看看上图右侧的梯度上升算法回归效果,梯度上升算法每次更新回归系数都要遍历整个数据集。从图中可以看出,当迭代次数为300多次的时候,回归系数才收敛。凑个整,就当它在遍历整个数据集300次的时候已经收敛好了。
没有对比就没有伤害,改进的随机梯度上升算法,在遍历数据集的第20次开始收敛。而梯度上升算法,在遍历数据集的第300次才开始收敛。
感谢代码原作者:Jack-Cui
原文链接:http://blog.csdn.net/c406495762/article/details/77723333
原文链接:http://blog.csdn.net/c406495762/article/details/77851973