文章原创,最近更新:2018-09-19
前言:
本文介绍机器学习分类算法中的Logistic回归算法并给出伪代码,Python代码实现。
学习参考链接:
1.第5章 Logistic回归
本章节的主要内容是:
项目案例1: 使用 Logistic 回归在简单数据集上的分类。
1.Logistic回归项目案例介绍:
项目案例1:
使用 Logistic 回归在简单数据集上的分类
项目概述:
在一个简单的数据集上,采用梯度上升法找到 Logistic 回归分类器在此数据集上的最佳回归系数.
Logistic 回归 工作原理
每个回归系数初始化为 1
重复 R 次:
计算整个数据集的梯度
使用 步长 x 梯度 更新回归系数的向量
返回回归系数
开发流程:
- 收集数据: 可以使用任何方法
- 准备数据: 由于需要进行距离计算,因此要求数据类型为数值型。另外,结构化数据格式则最佳
- 分析数据: 画出决策边界
- 训练算法: 使用梯度上升找到最佳参数
- 测试算法: 使用 Logistic 回归进行分类
- 使用算法: 对简单数据集中数据进行分类
Logistic 回归 算法特点:
- 优点: 计算代价不高,易于理解和实现。
- 缺点: 容易欠拟合,分类精度可能不高。
适用数据类型: 数值型和标称型数据。
数据集介绍
我们采用存储在 TestSet.txt 文本文件中的数据,存储格式如下:
-0.017612 14.053064 0
-1.395634 4.662541 1
-0.752157 6.538620 0
-1.322371 7.152853 0
0.423363 11.054677 0
2.改善前-项目相关代码
2.1 loadDataSet()
这段代码主要是解析数据
现在需要将截图1中的数据集,如何转换成截图2两组列表数据集?其中列表1的数据类型是float,列表2的数据类型是int.
具体方法见代码如下:
def loadDataSet(file_name):
'''
Desc:
加载并解析数据
Args:
file_name -- 要解析的文件路径
Returns:
dataMat -- 原始数据的特征
labelMat -- 原始数据的标签,也就是每条样本对应的类别。即目标向量
'''
# dataMat为原始数据, labelMat为原始数据的标签
dataMat = []
labelMat = []
fr = open(file_name)
for line in fr.readlines():
lineArr = line.strip().split()
# 为了方便计算,我们将 X0 的值设为 1.0 ,也就是在每一行的开头添加一个 1.0 作为 X0
dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])
labelMat.append(int(lineArr[2]))
return dataMat, labelMat
测试代码及其结果如下:
filename="testSet.txt"
dataMat,labelMat=loadDataSet(filename)
print(dataMat)
[[1.0, -0.017612, 14.053064], [1.0, -1.395634, 4.662541], [1.0, -0.752157, 6.53862], [1.0, -1.322371, 7.152853], [1.0, 0.423363, 11.054677], [1.0, 0.406704, 7.067335], [1.0, 0.667394, 12.741452], [1.0, -2.46015, 6.866805], [1.0, 0.569411, 9.548755], [1.0, -0.026632, 10.427743], [1.0, 0.850433, 6.920334], [1.0, 1.347183, 13.1755], [1.0, 1.176813, 3.16702], [1.0, -1.781871, 9.097953], [1.0, -0.566606, 5.749003], [1.0, 0.931635, 1.589505], [1.0, -0.024205, 6.151823], [1.0, -0.036453, 2.690988], [1.0, -0.196949, 0.444165], [1.0, 1.014459, 5.754399], [1.0, 1.985298, 3.230619], [1.0, -1.693453, -0.55754], [1.0, -0.576525, 11.778922], [1.0, -0.346811, -1.67873], [1.0, -2.124484, 2.672471], [1.0, 1.217916, 9.597015], [1.0, -0.733928, 9.098687], [1.0, -3.642001, -1.618087], [1.0, 0.315985, 3.523953], [1.0, 1.416614, 9.619232], [1.0, -0.386323, 3.989286], [1.0, 0.556921, 8.294984], [1.0, 1.224863, 11.58736], [1.0, -1.347803, -2.406051], [1.0, 1.196604, 4.951851], [1.0, 0.275221, 9.543647], [1.0, 0.470575, 9.332488], [1.0, -1.889567, 9.542662], [1.0, -1.527893, 12.150579], [1.0, -1.185247, 11.309318], [1.0, -0.445678, 3.297303], [1.0, 1.042222, 6.105155], [1.0, -0.618787, 10.320986], [1.0, 1.152083, 0.548467], [1.0, 0.828534, 2.676045], [1.0, -1.237728, 10.549033], [1.0, -0.683565, -2.166125], [1.0, 0.229456, 5.921938], [1.0, -0.959885, 11.555336], [1.0, 0.492911, 10.993324], [1.0, 0.184992, 8.721488], [1.0, -0.355715, 10.325976], [1.0, -0.397822, 8.058397], [1.0, 0.824839, 13.730343], [1.0, 1.507278, 5.027866], [1.0, 0.099671, 6.835839], [1.0, -0.344008, 10.717485], [1.0, 1.785928, 7.718645], [1.0, -0.918801, 11.560217], [1.0, -0.364009, 4.7473], [1.0, -0.841722, 4.119083], [1.0, 0.490426, 1.960539], [1.0, -0.007194, 9.075792], [1.0, 0.356107, 12.447863], [1.0, 0.342578, 12.281162], [1.0, -0.810823, -1.466018], [1.0, 2.530777, 6.476801], [1.0, 1.296683, 11.607559], [1.0, 0.475487, 12.040035], [1.0, -0.783277, 11.009725], [1.0, 0.074798, 11.02365], [1.0, -1.337472, 0.468339], [1.0, -0.102781, 13.763651], [1.0, -0.147324, 2.874846], [1.0, 0.518389, 9.887035], [1.0, 1.015399, 7.571882], [1.0, -1.658086, -0.027255], [1.0, 1.319944, 2.171228], [1.0, 2.056216, 5.019981], [1.0, -0.851633, 4.375691], [1.0, -1.510047, 6.061992], [1.0, -1.076637, -3.181888], [1.0, 1.821096, 10.28399], [1.0, 3.01015, 8.401766], [1.0, -1.099458, 1.688274], [1.0, -0.834872, -1.733869], [1.0, -0.846637, 3.849075], [1.0, 1.400102, 12.628781], [1.0, 1.752842, 5.468166], [1.0, 0.078557, 0.059736], [1.0, 0.089392, -0.7153], [1.0, 1.825662, 12.693808], [1.0, 0.197445, 9.744638], [1.0, 0.126117, 0.922311], [1.0, -0.679797, 1.22053], [1.0, 0.677983, 2.556666], [1.0, 0.761349, 10.693862], [1.0, -2.168791, 0.143632], [1.0, 1.38861, 9.341997], [1.0, 0.317029, 14.739025]]
print(labelMat)
[0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0]
2.2 sigmoid()
这个函数主要是定义sigmoid阶跃函数.
什么叫Sigmoid 函数呢?
首先我们先来看看这个函数的具体计算公式,如下:
下图给出了 Sigmoid 函数在不同坐标尺度下的两条曲线图。当 x 为 0 时,Sigmoid 函数值为 0.5 。随着 x 的增大,对应的 Sigmoid 值将逼近于 1 ; 而随着 x 的减小, Sigmoid 值将逼近于 0 。如果横坐标刻度足够大, Sigmoid 函数看起来很像一个阶跃函数。
那么如何用程序代码实现这个函数呢?
下面来看看具体的方法,如下:
# sigmoid阶跃函数
def sigmoid(inX):
# return 1.0 / (1 + exp(-inX))
return 1.0/(1+np.exp(-inX))
涉及到相关知识点
知识点1:
Python exp() 函数描述: exp() 函数返回 x 的指数,
1)语法
import math
math.exp(x)
注意:exp() 是不能直接访问的,需导入 math 模块,通过静态对象调用该方法。
- 参数
x 数值表达式 - 返回值
返回 x 的指数
2)案例:
# -*- coding: UTF-8 -*-
import numpy as np # 导入 math 模块
print ("np.exp(2) : ", np.exp(2))
print ("np.exp(100.12) : ", np.exp(100.12))
print ("np.exp(100.72) : ", np.exp(100.72))
print( "np.exp(119L) : ", np.exp(119L))
输出结果如下:
np.exp(2) : 7.38905609893
np.exp(100.12) : 3.03084361407e+43
np.exp(100.72) : 5.52255713025e+43
np.exp(119L) : 4.7978133273e+51
2.3 gradAscent()
这个函数主要是定义Logistic 回归梯度上升优化算法
具体方法见代码如下:
def gradAscent(dataMatIn, classLabels):
#用mat函数将列表转换成numpy矩阵
dataMatrix = np.mat(dataMatIn)
#用mat函数将列表转换成numpy矩阵,并进行转置
labelMat = np.mat(classLabels).transpose()
#返回dataMatrix的大小。m为行数,n为列数。
m,n = np.shape(dataMatrix)
#移动步长,也就是学习速率,控制更新的幅度。
alpha = 0.001
#最大迭代次数
maxCycles = 500
# weights 代表回归系数, 此处的 ones((n,1)) 创建一个长度和特征数相同的矩阵,其中的数全部都是 1
weights = np.ones((n,1))
for k in range(maxCycles):
# m*3 的矩阵 * 3*1 的矩阵 = m*1的矩阵
h = sigmoid(dataMatrix*weights)
# labelMat是实际值,labelMat - 是两个向量相减
error = (labelMat - h)
# 0.001* (3*m)*(m*1) 表示在每一个列上的一个误差情况,最后得出 x1,x2,xn的系数的偏移量
# 矩阵乘法,最后得到回归系数
weights = weights + alpha * dataMatrix.transpose() * error
return np.array(weights)
测试代码及其结果如下:
filename="testSet.txt"
dataMat,labelMat=loadDataSet(filename)
a=gradAscent(dataMat,labelMat)
print(a)
[[ 4.12414349]
[ 0.48007329]
[-0.6168482 ]]
相关知识点
- 知识点1:这个代码weights = weights + alpha * dataMatrix.transpose() * error为什么这样写呢?
一起看一下证明的过程,具体如下:
1)Sigmoid函数求导
这里
2)迭代更新公式
proof:
这里我们先规定logsitic函数的形式:
然后分析分类问题,y有两种取值0和1,根据这个概率p ( y=1 | x,w ) = h(x) , 所以p ( y=0 | x,w ) = 1-h(x),我们把两个概率结合起来,得到:
现在我们可以开始求Logistic回归系数了,因为logistic回归可以看做是一种概率模型,y发生的概率与w有关,因而我们可以对w使用极大似然估计(极大似然估计的知识推荐看《统计学完全教程》,讲的比较全面),使y发生的概率最大,此时的w便是我们寻找的最优回归参数,现根据概率密度写出极大似然函数:
取对数得到对数极大似然函数:
现在要使用梯度上升法w(即迭代更新公式),则我们需要求出对数极大似然函数的梯度:
这里yi对应书中labelMat的标签值,h(xi)对应sigmoid函数的输出值,将梯度带入公式:
这就得到了书中梯度上升法的最后一步迭代,这里alpha是步长,决定梯度上升的快慢。
学习参考链接:Logistic回归-数学原理(1)机器学习实战
- 知识点2:矩阵乘法
小明今天要做饭,消耗2斤肉,1斤蔬菜。肉每斤20元,蔬菜每斤5元,则一共需多少花费?
我们用向量相乘的方法写出来:
如果小明第二天有另一种做饭的方法,需要消耗1斤肉,4斤蔬菜,那么这两种方法的花费各是多少呢?我们显然需要另算这第二种方法的花费。把这个做饭方式写在第二个矩阵(向量是宽度或长度为1的矩阵)里:
小明家附近还有另一个菜市场,那里肉每斤15元,蔬菜每斤10元。那么,小明如果去这个菜市场,花费又是多少呢(分别计算上述两种做饭方式)?我们把这另外的一种价格写进第一个矩阵里:
这样我们看到了一个矩阵乘法的例子。在左边的这个矩阵的每一行,都代表了一种价目表;在右边的矩阵的每一列,都代表了一种做饭方式。那么所有可能的组合所最终产生的花费,则在结果矩阵中表示出来了。
学习参考链接:矩阵乘法的本质是什么?
2.4 plotBestFit()
这个函数主要是画出数据集和 Logistic 回归最佳拟合直线的函数
# 画出数据集和 Logistic 回归最佳拟合直线的函数
def plotBestFit(dataArr, labelMat, weights):
n = np.shape(dataArr)[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])
else:
xcord2.append(dataArr[i][1]); ycord2.append(dataArr[i][2])
fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(xcord1, ycord1, s=30, c='red', marker='s')
ax.scatter(xcord2, ycord2, s=30, c='green')
x = np.arange(-3.0, 3.0, 0.1)
y = (-weights[0]-weights[1]*x)/weights[2]
ax.plot(x, y)
plt.xlabel('X'); plt.ylabel('Y')
plt.show()
测试代码及其结果如下:
filename="testSet.txt"
dataArr,labelMat=loadDataSet(filename)
weights=gradAscent(dataArr,labelMat)
plotBestFit(dataArr, labelMat, weights)
相关知识点:
- 知识点1:
y = (-weights[0]-weights[1]*x)/weights[2]这段代码的由来?
首先理论上是这个样子的。
dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])
w0*x0+w1*x1+w2*x2=f(x)
x0最开始就设置为1, x2就是我们画图的y值,而f(x)被我们磨合误差给算到w0,w1,w2身上去了,所以:
w0+w1*x+w2*y=0 => y = (-w0-w1*x)/w2
2.5 testLR()
这段代码主要是测试算法: 使用 Logistic 回归进行分类
备注:这一段相当于2.4plotBestFit()的测试代码,只不过简化成了函数.
def testLR():
dataMat, labelMat = loadDataSet("TestSet.txt")
dataArr = np.array(dataMat)
weights = gradAscent(dataArr, labelMat)
plotBestFit(dataArr, labelMat, weights)
测试代码及其结果如下:
testLR()
2.6 完整的代码如下:
import matplotlib.pyplot as plt
import numpy as np
#解析数据
def loadDataSet(file_name):
'''
Desc:
加载并解析数据
Args:
file_name -- 要解析的文件路径
Returns:
dataMat -- 原始数据的特征
labelMat -- 原始数据的标签,也就是每条样本对应的类别。即目标向量
'''
# dataMat为原始数据, labelMat为原始数据的标签
dataMat = []
labelMat = []
fr = open(file_name)
for line in fr.readlines():
lineArr = line.strip().split()
# 为了方便计算,我们将 X0 的值设为 1.0 ,也就是在每一行的开头添加一个 1.0 作为 X0
dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])
labelMat.append(int(lineArr[2]))
return dataMat, labelMat
# sigmoid阶跃函数
def sigmoid(inX):
# return 1.0 / (1 + exp(-inX))
return 1.0/(1+np.exp(-inX))
# Logistic 回归梯度上升优化算法
def gradAscent(dataMatIn, classLabels):
#用mat函数将列表转换成numpy矩阵
dataMatrix = np.mat(dataMatIn)
#用mat函数将列表转换成numpy矩阵,并进行转置
labelMat = np.mat(classLabels).transpose()
#返回dataMatrix的大小。m为行数,n为列数。
m,n = np.shape(dataMatrix)
#移动步长,也就是学习速率,控制更新的幅度。
alpha = 0.001
#最大迭代次数
maxCycles = 500
# weights 代表回归系数, 此处的 ones((n,1)) 创建一个长度和特征数相同的矩阵,其中的数全部都是 1
weights = np.ones((n,1))
for k in range(maxCycles):
# m*3 的矩阵 * 3*1 的矩阵 = m*1的矩阵
h = sigmoid(dataMatrix*weights)
# labelMat是实际值,labelMat - 是两个向量相减
error = (labelMat - h)
# 0.001* (3*m)*(m*1) 表示在每一个列上的一个误差情况,最后得出 x1,x2,xn的系数的偏移量
# 矩阵乘法,最后得到回归系数
weights = weights + alpha * dataMatrix.transpose() * error
return np.array(weights)
# 画出数据集和 Logistic 回归最佳拟合直线的函数
def plotBestFit(dataArr, labelMat, weights):
n = np.shape(dataArr)[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])
else:
xcord2.append(dataArr[i][1]); ycord2.append(dataArr[i][2])
fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(xcord1, ycord1, s=30, c='red', marker='s')
ax.scatter(xcord2, ycord2, s=30, c='green')
x = np.arange(-3.0, 3.0, 0.1)
y = (-weights[0]-weights[1]*x)/weights[2]
ax.plot(x, y)
plt.xlabel('X'); plt.ylabel('Y')
plt.show()
#测试算法: 使用 Logistic 回归进行分类
def testLR():
dataMat, labelMat = loadDataSet("TestSet.txt")
dataArr = np.array(dataMat)
weights = gradAscent(dataArr, labelMat)
plotBestFit(dataArr, labelMat, weights)
testLR()
3.改善后-项目相关代码
梯度上升算法在每次更新回归系数时都需要遍历整个数据集,该方法在处理 100 个左右的数据集时尚可,但如果有数十亿样本和成千上万的特征,那么该方法的计算复杂度就太高了。一种改进方法是一次仅用一个样本点来更新回归系数,该方法称为 随机梯度上升算法。由于可以在新样本到来时对分类器进行增量式更新,因而随机梯度上升算法是一个在线学习(online learning)算法。与 “在线学习” 相对应,一次处理所有数据被称作是 “批处理” (batch) 。
随机梯度上升算法可以写成如下的伪代码:
#随机梯度上升算法-改善前
每个回归系数初始化为 1
重复 R 次:
计算整个数据集的梯度
使用 步长 x 梯度 更新回归系数的向量
返回回归系数
#随机梯度上升算法-改善后
所有回归系数初始化为 1
对数据集中每个样本
计算该样本的梯度
使用 alpha x gradient 更新回归系数值
返回回归系数值
以下是随机梯度上升算法的实现代码:
改善前
# Logistic 回归梯度上升优化算法
def gradAscent(dataMatIn, classLabels):
#用mat函数将列表转换成numpy矩阵
dataMatrix = np.mat(dataMatIn)
#用mat函数将列表转换成numpy矩阵,并进行转置
labelMat = np.mat(classLabels).transpose()
#返回dataMatrix的大小。m为行数,n为列数。
m,n = np.shape(dataMatrix)
#移动步长,也就是学习速率,控制更新的幅度。
alpha = 0.001
#最大迭代次数
maxCycles = 500
# weights 代表回归系数, 此处的 ones((n,1)) 创建一个长度和特征数相同的矩阵,其中的数全部都是 1
weights = np.ones((n,1))
for k in range(maxCycles):
# m*3 的矩阵 * 3*1 的矩阵 = m*1的矩阵
h = sigmoid(dataMatrix*weights)
# labelMat是实际值,labelMat - 是两个向量相减
error = (labelMat - h)
# 0.001* (3*m)*(m*1) 表示在每一个列上的一个误差情况,最后得出 x1,x2,xn的系数的偏移量
# 矩阵乘法,最后得到回归系数
weights = weights + alpha * dataMatrix.transpose() * error
return np.array(weights)
测试及其代码如下:
dataMat, labelMat = loadDataSet("TestSet.txt")
dataArr = np.array(dataMat)
weights = gradAscent(dataArr, labelMat)
plotBestFit(dataArr, labelMat, weights)
改善后1
# Logistic 回归梯度上升优化算法
def gradAscent0(dataMatIn, classLabels):
#返回dataMatrix的大小。m为行数,n为列数。
m,n = np.shape(dataMatIn)
#移动步长,也就是学习速率,控制更新的幅度。
alpha = 0.01
#函数ones创建一个n*1的数组
weights = np.ones(n)
for i in range(m):
#sum(dataMatrix[i]*weights)为了求 f(x)的值, f(x)=a1*x1+a2*x2+..+an*xn,此处求出的 h 是一个具体的数值,而不是一个矩阵
h = sigmoid(np.sum(dataMatIn[i]*weights))
# 计算真实类别与预测类别之间的差值,然后按照该差值调整回归系数
error = (classLabels[i] - h)
weights = weights + alpha * error * dataMatIn[i]
return weights
测试及其代码如下:
dataMat, labelMat = loadDataSet("TestSet.txt")
dataArr = np.array(dataMat)
weights = gradAscent0(dataArr, labelMat)
plotBestFit(dataArr, labelMat, weights)
可以看到,随机梯度上升算法与梯度上升算法在代码上很相似,但也有一些区别: 第一,后者的变量 h 和误差 error 都是向量,而前者则全是数值;第二,前者没有矩阵的转换过程,所有变量的数据类型都是 NumPy 数组。
判断优化算法优劣的可靠方法是看它是否收敛,也就是说参数是否达到了稳定值,是否还会不断地变化?下图展示了随机梯度上升算法在 200 次迭代过程中回归系数的变化情况。其中的系数2,也就是 X2 只经过了 50 次迭代就达到了稳定值,但系数 1 和 0 则需要更多次的迭代。如下图所示:
针对这个问题,我们改进了之前的随机梯度上升算法,如下:
改善后2
# 随机梯度上升算法(随机化)
def stocGradAscent1(dataMatIn, classLabels, numIter=150):
m,n = np.shape(dataMatIn)
# 创建与列数相同的矩阵的系数矩阵,1行3列
weights = np.ones(n)
# 随机梯度, 循环150,观察是否收敛
for j in range(numIter):
# [0, 1, 2 .. m-1]
dataIndex = list(range(m))
for i in range(m):
# i和j的不断增大,导致alpha的值不断减少,但是不为0
alpha = 4/(1.0+j+i)+0.01 # alpha 会随着迭代不断减小,但永远不会减小到0,因为后边还有一个常数项0.0001
# 随机产生一个 0~len()之间的一个值
# random.uniform(x, y) 方法将随机生成下一个实数,它在[x,y]范围内,x是这个范围内的最小值,y是这个范围内的最大值。
randIndex = int(random.uniform(0,len(dataIndex)))
# sum(dataMatrix[i]*weights)为了求 f(x)的值, f(x)=a1*x1+b2*x2+..+nn*xn
h = sigmoid(sum(dataMatIn[dataIndex[randIndex]]*weights))
error = classLabels[dataIndex[randIndex]] - h
weights = weights + alpha * error * dataMatIn[dataIndex[randIndex]]
del(dataIndex[randIndex])
return weights
测试及其代码如下:
dataMat, labelMat = loadDataSet("TestSet.txt")
dataArr = np.array(dataMat)
weights = gradAscent1(dataArr, labelMat)
plotBestFit(dataArr, labelMat, weights)
上面的改进版随机梯度上升算法,我们修改了两处代码。
第一处改进为 alpha 的值。alpha 在每次迭代的时候都会调整,这回缓解上面波动图的数据波动或者高频波动。另外,虽然 alpha 会随着迭代次数不断减少,但永远不会减小到 0,因为我们在计算公式中添加了一个常数项。
第二处修改为 randIndex 更新,这里通过随机选取样本拉来更新回归系数。这种方法将减少周期性的波动。这种方法每次随机从列表中选出一个值,然后从列表中删掉该值(再进行下一次迭代)。
程序运行之后能看到类似于上图的结果图。
改善后完整的代码如下:
import matplotlib.pyplot as plt
import numpy as np
import random
#解析数据
def loadDataSet(file_name):
'''
Desc:
加载并解析数据
Args:
file_name -- 要解析的文件路径
Returns:
dataMat -- 原始数据的特征
labelMat -- 原始数据的标签,也就是每条样本对应的类别。即目标向量
'''
# dataMat为原始数据, labelMat为原始数据的标签
dataMat = []
labelMat = []
fr = open(file_name)
for line in fr.readlines():
lineArr = line.strip().split()
# 为了方便计算,我们将 X0 的值设为 1.0 ,也就是在每一行的开头添加一个 1.0 作为 X0
dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])
labelMat.append(int(lineArr[2]))
return dataMat, labelMat
# sigmoid阶跃函数
def sigmoid(inX):
# return 1.0 / (1 + exp(-inX))
return 1.0/(1+np.exp(-inX))
# Logistic 回归梯度上升优化算法
def gradAscent(dataMatIn, classLabels):
#用mat函数将列表转换成numpy矩阵
dataMatrix = np.mat(dataMatIn)
#用mat函数将列表转换成numpy矩阵,并进行转置
labelMat = np.mat(classLabels).transpose()
#返回dataMatrix的大小。m为行数,n为列数。
m,n = np.shape(dataMatrix)
#移动步长,也就是学习速率,控制更新的幅度。
alpha = 0.001
#最大迭代次数
maxCycles = 500
# weights 代表回归系数, 此处的 ones((n,1)) 创建一个长度和特征数相同的矩阵,其中的数全部都是 1
weights = np.ones((n,1))
for k in range(maxCycles):
# m*3 的矩阵 * 3*1 的矩阵 = m*1的矩阵
h = sigmoid(dataMatrix*weights)
# labelMat是实际值,labelMat - 是两个向量相减
error = (labelMat - h)
# 0.001* (3*m)*(m*1) 表示在每一个列上的一个误差情况,最后得出 x1,x2,xn的系数的偏移量
# 矩阵乘法,最后得到回归系数
weights = weights + alpha * dataMatrix.transpose() * error
return np.array(weights)
# Logistic 回归梯度上升优化算法
def gradAscent0(dataMatIn, classLabels):
#返回dataMatrix的大小。m为行数,n为列数。
m,n = np.shape(dataMatIn)
#移动步长,也就是学习速率,控制更新的幅度。
alpha = 0.01
#函数ones创建一个n*1的数组
weights = np.ones(n)
for i in range(m):
#sum(dataMatrix[i]*weights)为了求 f(x)的值, f(x)=a1*x1+a2*x2+..+an*xn,此处求出的 h 是一个具体的数值,而不是一个矩阵
h = sigmoid(np.sum(dataMatIn[i]*weights))
# 计算真实类别与预测类别之间的差值,然后按照该差值调整回归系数
error = (classLabels[i] - h)
weights = weights + alpha * error * dataMatIn[i]
return weights
# 画出数据集和 Logistic 回归最佳拟合直线的函数
def plotBestFit(dataArr, labelMat, weights):
n = np.shape(dataArr)[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])
else:
xcord2.append(dataArr[i][1]); ycord2.append(dataArr[i][2])
fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(xcord1, ycord1, s=30, c='red', marker='s')
ax.scatter(xcord2, ycord2, s=30, c='green')
x = np.arange(-3.0, 3.0, 0.1)
y = (-weights[0]-weights[1]*x)/weights[2]
ax.plot(x, y)
plt.xlabel('X'); plt.ylabel('Y')
plt.show()
#测试算法: 使用 Logistic 回归进行分类
def testLR():
dataMat, labelMat = loadDataSet("TestSet.txt")
dataArr = np.array(dataMat)
weights = stocGradAscent1(dataArr, labelMat)
plotBestFit(dataArr, labelMat, weights)
# 随机梯度上升算法(随机化)
def stocGradAscent1(dataMatIn, classLabels, numIter=150):
m,n = np.shape(dataMatIn)
# 创建与列数相同的矩阵的系数矩阵,1行3列
weights = np.ones(n)
# 随机梯度, 循环150,观察是否收敛
for j in range(numIter):
# [0, 1, 2 .. m-1]
dataIndex = list(range(m))
for i in range(m):
# i和j的不断增大,导致alpha的值不断减少,但是不为0
alpha = 4/(1.0+j+i)+0.01 # alpha 会随着迭代不断减小,但永远不会减小到0,因为后边还有一个常数项0.0001
# 随机产生一个 0~len()之间的一个值
# random.uniform(x, y) 方法将随机生成下一个实数,它在[x,y]范围内,x是这个范围内的最小值,y是这个范围内的最大值。
randIndex = int(random.uniform(0,len(dataIndex)))
# sum(dataMatrix[i]*weights)为了求 f(x)的值, f(x)=a1*x1+b2*x2+..+nn*xn
h = sigmoid(sum(dataMatIn[dataIndex[randIndex]]*weights))
error = classLabels[dataIndex[randIndex]] - h
weights = weights + alpha * error * dataMatIn[dataIndex[randIndex]]
del(dataIndex[randIndex])
return weights
testLR()