一、逻辑回归代码
逻辑回归其实是求一个分类器,是二分类问题,其利用一个sigmoid函数去定义样本属于正类的概率,sigmoid函数的输入值z利用了线性回归的wx,所得出的数值是在[0,1]内的概率,即样本属于正类的概率值。在实践中,我们需要从训练样本中学习出一个w和b,这就是需要的逻辑回归模型,w和b包含在一个权重系数向量w中,w的更新是依据对损失函数即负对数似然函数求偏导数得梯度,然后每轮迭代中w根据负梯度来更新,梯度是一个向量,表示函数值上升最快的方向,由于这里的损失函数(负对数似然函数)是凸函数,因此在训练中求使其值最小的w值,就利用到梯度下降法(负梯度方向),当迭代到一定次数后所得到的w即为模型,预测时将利用sigmoid函数(此时输入值即wx,w是所得模型,x是预测样本)求得x属于正样本的概率(回归),若值大于阈值,即判定为正类(分类)。
关于回归预测模型和损失函数的定义,以及梯度更新公式的推导,
参考:https://mp.weixin.qq.com/s/BJxdDz7DQg5QIWzzgNQrgA
其梯度公式最后少了一个负号
其中每轮迭代中权重系数矩阵weights的更新运用到训练样本数据矩阵dataMatrix,训练样本标签矩阵labelMat和训练样本预测值矩阵h,这种向量化或矩阵化运算提高计算效率,避免for循环带来的计算开销,其矩阵计算的原理如下:这里样本特征是2,加上偏置项1,所以有3个特征,故其权重向量的维度也是3
完整的逻辑回归代码如下:
# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
def loadDataSet(file_name):
"""
:desc: 加载并解析数据
:param file_name: 文件名称,要解析的文件所在路径
:return:
dataMat -- 原始数据的特征
labelMat -- 原始数据的标签,也就是每条样本对应的类别
"""
dataMat = []
labelMat = []
fr = open(file_name)
for line in fr.readlines():
lineArr = line.strip().split()
if len(lineArr) == 1:
continue
dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])]) # 此处为每个样本添加一个特征1.0
labelMat.append(int(lineArr[2]))
return dataMat, labelMat
# sigmoid跃阶函数
def sigmoid(x):
return 1.0 / (1 + np.exp(-x))
# 批量梯度下降
def gradDescent(dataMatIn, classLabels):
"""
:param dataMatIn: 是一个2维Numpy数组,每行代表每个训练样本,每列分别代表样本不同特征
:param classLabels: 是一个List,是训练样本的类别标签,它是一个1*100的行向量。为了便于矩阵计算,需要将该行向量转换为列向量,
即将原向量转置后赋值给labelMat
:return: Numpy数组,权重系数
"""
# 将numpy数组dataMatIn转化为numpy矩阵[[1 1 2]
# [1 1 2]
# ... ]
dataMatrix = np.mat(dataMatIn) # 转换为Nupmy矩阵
# 将标签List转化为行矩阵[[0 1 0 1 0 1 ...]],并转置成列矩阵[[0]
# [1]
# [0]
# ..]
labelMat = np.mat(classLabels).transpose()
# m是训练样本个数,n是特征数
m, n = np.shape(dataMatrix)
# alpha代表移动步长
alpha = 0.001
# 迭代次数
maxCycles = 500
# 初始权重系数矩阵,长度与特征数相同[[1.]
# [1.]
# [1.]]
weights = np.ones((n, 1)) # numpy数组
for k in range(maxCycles):
h = sigmoid(dataMatrix*weights) # 对应于训练样本的预测值 100*1的列矩阵
errors = labelMat - h # 100*1的列矩阵
weights = weights + alpha * dataMatrix.transpose() * errors
return np.array(weights)
# 随机梯度下降
# 梯度下降算法在每次更新权重系数时都需要遍历整个数据集,计算复杂度较高
# 随机梯度下降每次更新权重只需要一个样本
def stocGradDescent(dataMatrix, classLabels):
"""
:param dataMatrix: Numpy2维数组
:param classLabels: List
:return: numpy数组,权重系数
"""
m, n = np.shape(dataMatrix) # m是第一维度的元素个数,即样本数,n是第二维度元素个数,即特征数
dataMat = np.mat(dataMatrix)
alpha = 0.001
weights = np.ones((n, 1))
print(np.shape(weights))
print(np.shape(dataMat[1]))
for i in range(m):
h = sigmoid(np.sum(dataMat[i]*weights)) # 取一条样本与权重系数相乘,这里的h是一个数值而不再是矩阵了
error = classLabels[i] - h
weights = weights + alpha * error * dataMat[i].transpose()
return np.array(weights)
# 改进随机梯度下降
def stocGradDescent1(dataMatrix, classLabels, numIter=150):
"""
:param dataMatrix: Numpy2维数组
:param classLabels: List
:param numIter: 迭代次数
:return:
"""
m, n = np.shape(dataMatrix)
dataMat = np.mat(dataMatrix)
weights = np.ones((n, 1))
for j in range(numIter):
dataIndex = list(range(m)) # 存放样本索引[0 1 2 .. m-1]
for i in range(m):
alpha = 4/(1.0+j+i)+0.0001 # alpha 会随着迭代不断减小
# 随机取一个样本的索引,其值在[0,len(dataIndex))左闭右开区间范围内
randIndex = int(np.random.uniform(0, len(dataIndex)))
h = sigmoid(np.sum(dataMat[dataIndex[randIndex]]*weights))
error = classLabels[dataIndex[randIndex]] - h
weights = weights + alpha * error * dataMat[dataIndex[randIndex]].transpose()
# 删除该随机索引
del dataIndex[randIndex]
return np.array(weights)
def plotBestFit(dataArr, labelMat, weights):
"""
Desc: 将训练样本和训练出的权重系数可视化
:param dataArr: numpy2维数组,样本数据的特征,[[1 1 2] [1 2 1] [1 1 2]]
:param labelMat: List,样本数据的类别标签
:param weights: numpy数组,权重系数
:return: None
"""
n = np.shape(dataArr)[0] # 样本个数
xcord1 = []; ycord1 = [] # 分别存放类别为1的样本特征
xcord2 = []; ycord2 = [] # 分别存放类别为2的样本特征
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.savefig("D://fig_2.png")
plt.show()
def testLR():
# 1.准备训练数据
dataMat, labelMat = loadDataSet("D:\\testSet.TXT")
# 2.训练模型
dataArr = np.array(dataMat)
# weights = gradDescent(dataArr, labelMat)
weights = stocGradDescent1(dataArr, labelMat)
print("The final weight is: \n", weights, "\n", type(weights), np.shape(weights))
# 数据可视化
plotBestFit(dataArr, labelMat, weights)
# 预测样本
def classifyVector(x, weights):
"""
:param x: numpy matrix, 待预测的样本特征向量
:param weights: numpy array, 根据训练样本学习由梯度下降得到的权重
:return: 如果prob大于0.5返回1,否则返回0
"""
prob = sigmoid(np.sum(x*weights))
if prob > 0.5:
return 1.0
else:
return 0.0
if __name__ == "__main__":
testLR()
其中用于学习权重系数有三种方法:分别是批量梯度下降gradDescent、随机梯度下降stocGradDescent和改进随机梯度下降stocGradDescent1。批量梯度下降也就是常说的梯度下降法,梯度下降算法在每次更新回归系数时都需要遍历整个数据集,该方法在处理 100 个左右的数据集时尚可,但如果有数十亿样本和成千上万的特征,那么该方法的计算复杂度就太高了。一种改进方法是一次仅用一个样本点来更新回归系数,该方法称为 随机梯度下降算法。相比之下,批量梯度下降比随机梯度下降具有更好的学习效果,学习得到的回归系数用作预测的效果更准确,下面是两种算法寻优后的分类超平面:
由于w0*x0+w1*x1+w2*x2=f(x), x0最开始就设置为1, x2就是我们画图的y值,令f(x)=0,
所以: w0+w1*x+w2*y=0 => y = (-w0-w1*x)/w2 即分类超平面直线
批量梯度下降:
随机梯度下降:
随机梯度下降每次迭代只用到一条样本数据,其下降方向并不是全局最优的,也就是不一定向着最低点的方向,因此虽然迭代速度快,但是很容易陷入局部最优,最终的学习效果就没有批量梯度下降好。
因此对随机梯度下降作一次改进,代码里是stocGradDescent1,第一处改进为 alpha 的值。alpha 在每次迭代的时候都会调整,会随着迭代次数不断减少,但永远不会减小到 0,因为在计算公式中添加了一个常数项。
第二处修改为 randIndex 更新,这里通过随机选取样本来更新回归系数。这种方法每次随机从列表中选出一个值,然后从列表中删掉该值(再进行下一次迭代)。。
改进随机梯度下降:
在机器学习领域,梯度下降分为3类:
1、批量梯度下降BGD:每轮迭代计算所有样本数据,所有样本共同决定最优下降方向
优点:逼近全局最优
缺点:计算量大,批量完成,速度慢,样本数量大时不适用
2、随机梯度下降SGD:每轮迭代从训练样本中抽取一个样本进行计算更新,每次都不用遍历所有数据集
优点:速度快,逐条样本进行计算
缺点:陷入局部最优,下降幅度相对慢,需要迭代更多次数,因为每次选取的方向不一定是最优
3、小批量的梯度下降MBGD:每轮迭代计算所有数据中的子集,是1和2的中和方案,节省时间,并且寻优也准确
二、从疝气病症预测病马的死亡率
使用逻辑回归来预测患有疝病的马的存活问题。疝病是描述马胃肠痛的术语。然而,这种病不一定源自马的胃肠问题,其他问题也可能引发马疝病。这个数据集中包含了医院检测马疝病的一些指标,分为训练集和测试集,每条数据 共22个字段,最后一个字段是其类别标签0或1,项目代码如下
病马数据集:
链接:https://pan.baidu.com/s/19qLbQ2TL757r0zyJjbLTUQ
提取码:gw3j
# 预测样本
def classifyVector(x, weights):
"""
:param x: numpy matrix, 待预测的样本特征向量
:param weights: numpy array, 根据训练样本学习由梯度下降得到的权重
:return: 如果prob大于0.5返回1,否则返回0
"""
prob = sigmoid(np.sum(x*weights))
if prob > 0.5:
return 1.0
else:
return 0.0
# 从疝气病症预测病马的死亡率
def colicTest():
"""
Desc:读入测试集和训练集,并对数据进行格式化处理
:return: errorRate 分类错误率
"""
frTrain = open('D:\\horseColicTraining.txt') # 训练样本
frTest = open('D:\\horseColicTest.txt') # 测试样本
trainingSet = [] # [[],[],[],...]
trainingLabels = []
for line in frTrain.readlines():
currLine = line.strip().split('\t') # 得到每一条数据字符串列表(共22个字段,最后一个字段为标签)
lineArr = [] # 存储每条样本特征
for i in range(len(currLine)-1):
lineArr.append(float(currLine[i]))
trainingSet.append(lineArr)
trainingLabels.append(float(currLine[len(currLine)-1]))
# 使用改进后的随机梯度下降算法求最佳权重系数
# trainWeights = stocGradDescent1(np.array(trainingSet), trainingLabels, 500) # 改进后随机梯度下降
# trainWeights = stocGradDescent(np.array(trainingSet), trainingLabels) # 随机梯度下降
trainWeights = gradDescent(np.array(trainingSet), trainingLabels) # 梯度下降
errorCount = 0
numTestVec = 0.0
# 读取测试样本进行测试,计算分类错误的样本条数和最终的错误率
for line in frTest.readlines():
numTestVec += 1.0 # 测试样本数量加1
currLine = line.strip().split('\t')
lineArr = [] # 存储每条样本特征
for i in range(len(currLine)-1):
lineArr.append(float(currLine[i]))
if int(classifyVector(np.mat(lineArr), trainWeights)) != int(currLine[len(currLine)-1]):
errorCount += 1 # 预测错误,错误数加1
errorRate = float(errorCount) / numTestVec
print("测试的错误率是: ", errorRate)
return errorRate
if __name__ == "__main__":
colicTest()
运行结果:
批量梯度下降的错误率:0.29850746268656714
随机梯度下降的错误率:0.47761194029850745
改进随机梯度的错误率:0.3283582089552239