完整代码地址:https://github.com/apachecn/MachineLearning/blob/master/src/python/5.Logistic/logistic.py
项目案例2: 从疝气病症预测病马的死亡率
项目概述
使用 Logistic 回归来预测患有疝病的马的存活问题。疝病是描述马胃肠痛的术语。然而,这种病不一定源自马的胃肠问题,其他问题也可能引发马疝病。这个数据集中包含了医院检测马疝病的一些指标,有的指标比较主观,有的指标难以测量,例如马的疼痛级别。
开发流程
收集数据: 给定数据文件
准备数据: 用 Python 解析文本文件并填充缺失值
分析数据: 可视化并观察数据
训练算法: 使用优化算法,找到最佳的系数
测试算法: 为了量化回归的效果,需要观察错误率。根据错误率决定是否回退到训练阶段,
通过改变迭代的次数和步长的参数来得到更好的回归系数
使用算法: 实现一个简单的命令行程序来手机马的症状并输出预测结果并非难事,
这可以作为留给大家的一道习题
收集数据: 给定数据文件
病马的训练数据已经给出来了,如下形式存储在文本文件中:
1.000000 1.000000 39.200000 88.000000 20.000000 0.000000 0.000000 4.000000 1.000000 3.000000 4.000000 2.000000 0.000000 0.000000 0.000000 4.000000 2.000000 50.000000 85.000000 2.000000 2.000000 0.000000
2.000000 1.000000 38.300000 40.000000 24.000000 1.000000 1.000000 3.000000 1.000000 3.000000 3.000000 1.000000 0.000000 0.000000 0.000000 1.000000 1.000000 33.000000 6.700000 0.000000 0.000000 1.000000
准备数据: 用 Python 解析文本文件并填充缺失值
处理数据中的缺失值
假设有100个样本和20个特征,这些数据都是机器收集回来的。若机器上的某个传感器损坏导致一个特征无效时该怎么办?此时是否要扔掉整个数据?这种情况下,另外19个特征怎么办? 它们是否还可以用?答案是肯定的。因为有时候数据相当昂贵,扔掉和重新获取都是不可取的,所以必须采用一些方法来解决这个问题。
下面给出了一些可选的做法:
- 使用可用特征的均值来填补缺失值;
- 使用特殊值来填补缺失值,如 -1;
- 忽略有缺失值的样本;
- 使用有相似样本的均值添补缺失值;
- 使用另外的机器学习算法预测缺失值。
现在,我们对下一节要用的数据集进行预处理,使其可以顺利地使用分类算法。在预处理需要做两件事:
- 所有的缺失值必须用一个实数值来替换,因为我们使用的 NumPy 数据类型不允许包含缺失值。我们这里选择实数 0 来替换所有缺失值,恰好能适用于 Logistic 回归。这样做的直觉在于,我们需要的是一个在更新时不会影响系数的值。回归系数的更新公式如下:
weights = weights + alpha * error * dataMatrix[randIndex]
如果 dataMatrix 的某个特征对应值为 0,那么该特征的系数将不做更新,即:
weights = weights
另外,由于 Sigmoid(0) = 0.5 ,即它对结果的预测不具有任何倾向性,因此我们上述做法也不会对误差造成任何影响。基于上述原因,将缺失值用 0 代替既可以保留现有数据,也不需要对优化算法进行修改。此外,该数据集中的特征取值一般不为 0,因此在某种意义上说它也满足 “特殊值” 这个要求。
- 如果在测试数据集中发现了一条数据的类别标签已经缺失,那么我们的简单做法是将该条数据丢弃。这是因为类别标签与特征不同,很难确定采用某个合适的值来替换。采用 Logistic 回归进行分类时这种做法是合理的,而如果采用类似 kNN 的方法就可能不太可行。
原始的数据集经过预处理后,保存成两个文件: horseColicTest.txt 和 horseColicTraining.txt 。
分析数据: 可视化并观察数据
将数据使用 MatPlotlib 打印出来,观察数据是否是我们想要的格式
训练算法: 使用优化算法,找到最佳的系数
下面给出 原始的梯度上升算法,随机梯度上升算法,改进版随机梯度上升算法 的代码:
# 正常的处理方案
# 两个参数:第一个参数==> dataMatIn 是一个2维NumPy数组,每列分别代表每个不同的特征,每行则代表每个训练样本。
# 第二个参数==> classLabels 是类别标签,它是一个 1*100 的行向量。为了便于矩阵计算,需要将该行向量转换为列向量,做法是将原向量转置,再将它赋值给labelMat。
def gradAscent(dataMatIn, classLabels):
# 转化为矩阵[[1,1,2],[1,1,2]....]
dataMatrix = mat(dataMatIn) # 转换为 NumPy 矩阵
# 转化为矩阵[[0,1,0,1,0,1.....]],并转制[[0],[1],[0].....]
# transpose() 行列转置函数
# 将行向量转化为列向量 => 矩阵的转置
labelMat = mat(classLabels).transpose() # 首先将数组转换为 NumPy 矩阵,然后再将行向量转置为列向量
# m->数据量,样本数 n->特征数
m,n = shape(dataMatrix)
# print m, n, '__'*10, shape(dataMatrix.transpose()), '__'*100
# alpha代表向目标移动的步长
alpha = 0.001
# 迭代次数
maxCycles = 500
# 生成一个长度和特征数相同的矩阵,此处n为3 -> [[1],[1],[1]]
# weights 代表回归系数, 此处的 ones((n,1)) 创建一个长度和特征数相同的矩阵,其中的数全部都是 1
weights = ones((n,1))
for k in range(maxCycles): #heavy on matrix operations
# m*3 的矩阵 * 3*1 的单位矩阵 = m*1的矩阵
# 那么乘上单位矩阵的意义,就代表:通过公式得到的理论值
# 参考地址: 矩阵乘法的本质是什么? https://www.zhihu.com/question/21351965/answer/31050145
# print 'dataMatrix====', dataMatrix
# print 'weights====', weights
# n*3 * 3*1 = n*1
h = sigmoid(dataMatrix*weights) # 矩阵乘法
# print 'hhhhhhh====', h
# labelMat是实际值
error = (labelMat - h) # 向量相减
# 0.001* (3*m)*(m*1) 表示在每一个列上的一个误差情况,最后得出 x1,x2,xn的系数的偏移量
weights = weights + alpha * dataMatrix.transpose() * error # 矩阵乘法,最后得到回归系数
return array(weights)
# 随机梯度上升
# 梯度上升优化算法在每次更新数据集时都需要遍历整个数据集,计算复杂都较高
# 随机梯度上升一次只用一个样本点来更新回归系数
def stocGradAscent0(dataMatrix, classLabels):
m,n = shape(dataMatrix)
alpha = 0.01
# n*1的矩阵
# 函数ones创建一个全1的数组
weights = ones(n) # 初始化长度为n的数组,元素全部为 1
for i in range(m):
# sum(dataMatrix[i]*weights)为了求 f(x)的值, f(x)=a1*x1+b2*x2+..+nn*xn,此处求出的 h 是一个具体的数值,而不是一个矩阵
h = sigmoid(sum(dataMatrix[i]*weights))
# print 'dataMatrix[i]===', dataMatrix[i]
# 计算真实类别与预测类别之间的差值,然后按照该差值调整回归系数
error = classLabels[i] - h
# 0.01*(1*1)*(1*n)
print weights, "*"*10 , dataMatrix[i], "*"*10 , error
weights = weights + alpha * error * dataMatrix[i]
return weights
# 随机梯度上升算法(随机化)
def stocGradAscent1(dataMatrix, classLabels, numIter=150):
m,n = shape(dataMatrix)
weights = ones(n) # 创建与列数相同的矩阵的系数矩阵,所有的元素都是1
# 随机梯度, 循环150,观察是否收敛
for j in range(numIter):
# [0, 1, 2 .. m-1]
dataIndex = range(m)
for i in range(m):
# i和j的不断增大,导致alpha的值不断减少,但是不为0
alpha = 4/(1.0+j+i)+0.0001 # 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(dataMatrix[randIndex]*weights))
error = classLabels[randIndex] - h
# print weights, '__h=%s' % h, '__'*20, alpha, '__'*20, error, '__'*20, dataMatrix[randIndex]
weights = weights + alpha * error * dataMatrix[randIndex]
del(dataIndex[randIndex])
return weights
测试算法: 为了量化回归的效果,需要观察错误率。根据错误率决定是否回退到训练阶段,通过改变迭代的次数和步长的参数来得到更好的回归系数
Logistic 回归分类函数
# 分类函数,根据回归系数和特征向量来计算 Sigmoid的值
def classifyVector(inX, weights):
'''
Desc:
最终的分类函数,根据回归系数和特征向量来计算 Sigmoid 的值,大于0.5函数返回1,否则返回0
Args:
inX -- 特征向量,features
weights -- 根据梯度下降/随机梯度下降 计算得到的回归系数
Returns:
如果 prob 计算大于 0.5 函数返回 1
否则返回 0
'''
prob = sigmoid(sum(inX*weights))
if prob > 0.5: return 1.0
else: return 0.0
# 打开测试集和训练集,并对数据进行格式化处理
def colicTest():
'''
Desc:
打开测试集和训练集,并对数据进行格式化处理
Args:
None
Returns:
errorRate -- 分类错误率
'''
frTrain = open('input/5.Logistic/horseColicTraining.txt')
frTest = open('input/5.Logistic/horseColicTest.txt')
trainingSet = []
trainingLabels = []
# 解析训练数据集中的数据特征和Labels
# trainingSet 中存储训练数据集的特征,trainingLabels 存储训练数据集的样本对应的分类标签
for line in frTrain.readlines():
currLine = line.strip().split(' ')
lineArr = []
for i in range(21):
lineArr.append(float(currLine[i]))
trainingSet.append(lineArr)
trainingLabels.append(float(currLine[21]))
# 使用 改进后的 随机梯度下降算法 求得在此数据集上的最佳回归系数 trainWeights
trainWeights = stocGradAscent1(array(trainingSet), trainingLabels, 500)
errorCount = 0
numTestVec = 0.0
# 读取 测试数据集 进行测试,计算分类错误的样本条数和最终的错误率
for line in frTest.readlines():
numTestVec += 1.0
currLine = line.strip().split(' ')
lineArr = []
for i in range(21):
lineArr.append(float(currLine[i]))
if int(classifyVector(array(lineArr), trainWeights)) != int(currLine[21]):
errorCount += 1
errorRate = (float(errorCount) / numTestVec)
print "the error rate of this test is: %f" % errorRate
return errorRate
# 调用 colicTest() 10次并求结果的平均值
def multiTest():
numTests = 10
errorSum = 0.0
for k in range(numTests):
errorSum += colicTest()
print "after %d iterations the average error rate is: %f" % (numTests, errorSum/float(numTests))
使用算法: 实现一个简单的命令行程序来手机马的症状并输出预测结果并非难事,这可以作为留给大家的一道习题
完整代码地址:https://github.com/apachecn/MachineLearning/blob/master/src/python/5.Logistic/logistic.py
优质内容筛选与推荐>>
1、谷歌Apps vs.微软Office2、fcitx输入法在wps、wineqq中失灵问题的解决3、08-urllib4、便签应用和笔记应用5、微信小程序demo解读(二)