Logistic回归应用-预测马的死亡率
本文所有代码均来自《机器学习实战》，数据也是
本例中的数据有以下几个特征：
- 部分指标比较主观、难以很好的定量测量，例如马的疼痛级别
- 数据集中有30%的数据是缺失的，例如：
如何处理数据缺失情况
数据是很宝贵的，对于有缺失的数据我们不能直接丢弃，而应该想办法将其修补后再利用起来
常用的方法有：
- 使用可用样本中该特征的均值来填补缺失样本中缺失特征的值
- 使用特殊值来填补缺失值，如-1
- 忽略有缺失值的样本
- 使用相似样本的均值填补缺失值
- 使用另外的机器学习算法来预测缺失值
但如果是Label信息丢失的话，一个简单的办法就是将该data丢弃，因为Label和attributes不一样，很难确定采用某个合适的值来替换缺失值。但是这种方法在类似KNN这样的、样本数量和样本完整度对结果有较为明显影响的算法中就不太好了。
回到本例，这里是采用了将缺失数据置换为0的操作。因为在Logistic回归的梯度下降法中，如果某个样本的一个attribute为0的话，那么就不使用它进行weights的更新了，这显然是一种很好的结果。另外，由于sigmoid(0)=0.5,既不偏向0也不偏向1，所以这个特征设置为0的话，该特征不会对结果的预测产生偏向性。
这是未处理的数据：
这是处理后的数据：
关于每个attribute代表什么意思，可以在这个网站上查到：http://archive.ics.uci.edu/ml/datasets/Horse+Colic。反正就是关于马的一些特征，由上图我们可以看到，这里我们并没有使用原数据中的全部特征，只选了一部分。
将Logistic回归用于分类
#用于返回分类结果（0或1）
def classifyVector(inX, weights):
prob = sigmoid(sum(inX*weights))
if prob > 0.5: return 1.0
else: return 0.0
def colicTest():
#训练集和测试集
frTrain = open('horseColicTraining.txt'); frTest = open('horseColicTest.txt')
trainingSet = []; trainingLabels = []
#调用训练集进行训练
for line in frTrain.readlines():
currLine = line.strip().split(' ')
lineArr =[]
#一共有21个特征
for i in range(21):
#lineArr以数组形式保存由字符转换为浮点数的一条数据
lineArr.append(float(currLine[i]))
trainingSet.append(lineArr)
trainingLabels.append(float(currLine[21]))
trainWeights = stocGradAscent1(array(trainingSet), trainingLabels, 1000)
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
完整代码（实现+测试）如下：
'''
Created on Oct 27, 2010
Logistic Regression Working Module
@author: Peter
'''
'''
注意，本例中默认data在testSet.txt中且有两个attributes，
使用该模型的时候要注意根据实际使用的attributes个数修改代码
'''
from numpy import *
def loadDataSet():#加载数据
dataMat = []; labelMat = []
fr = open('testSet.txt')
for line in fr.readlines():
lineArr = line.strip().split()
#Python strip() 方法用于移除字符串头尾指定的字符（默认为空格）或字符序列。
#注意，在这里为data添加了1的一列
dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])
#数据有两个attributes
labelMat.append(int(lineArr[2]))
return dataMat,labelMat
def sigmoid(inX):
return 1.0/(1+exp(-inX))
#梯度下降方法
def gradAscent(dataMatIn, classLabels):
#转换为numpy矩阵类型，便于优化和并行计算
dataMatrix = mat(dataMatIn) #convert to NumPy matrix
labelMat = mat(classLabels).transpose() #convert to NumPy matrix
#将Label转置为列向量
m,n = shape(dataMatrix)#m行n列
alpha = 0.001#学习速度
maxCycles = 500#最大学习步数
weights = ones((n,1))#将来的输出
for k in range(maxCycles): #heavy on matrix operations
h = sigmoid(dataMatrix*weights) #matrix mult
error = (labelMat - h) #vector subtraction
weights = weights + alpha * dataMatrix.transpose()* error #matrix mult
return weights
def plotBestFit(weights):
import matplotlib.pyplot as plt
dataMat,labelMat=loadDataSet()#加载源数据，用于画图
dataArr = array(dataMat)#将data向量化
n = shape(dataArr)[0] #n为data的个数
xcord1 = []; ycord1 = []
xcord2 = []; ycord2 = []
#给data分类
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')
#arange([start,] stop[, step,], dtype=None)
#根据start与stop指定的范围以及step设定的步长，生成一个 ndarray。
'''
这一段是在干什么呢？其实就是想要画一条直线出来，但是只能画点，由点组成线
'''
x = arange(-3.0, 3.0, 0.1)
y = (-weights[0]-weights[1]*x)/weights[2]
ax.plot(x, y)
plt.xlabel('X1'); plt.ylabel('X2');
plt.show()
#随机梯度下降法，只对整个数据集迭代了一次，显然不够
def stocGradAscent0(dataMatrix, classLabels):
m,n = shape(dataMatrix)
alpha = 0.01
weights = ones(n) #initialize to all ones
for i in range(m):
h = sigmoid(sum(dataMatrix[i]*weights))
error = classLabels[i] - h
#这里其实不是真正的随机，样本是按顺序取的
weights = weights + alpha * error * dataMatrix[i]
return weights
#改进的随机梯度下降法
def stocGradAscent1(dataMatrix, classLabels, numIter=150):
m,n = shape(dataMatrix)
weights = ones(n) #initialize to all ones
for j in range(numIter):#对样本迭代150次
dataIndex = list(range(m))#dataIndex保存了所有在本次迭代中还没有使用过的data
for i in range(m):
'''
改进1：学习速度在每次迭代后都会调整，越往后学习速度越小，这会缓解在最优点附近的波动.
由于我们在alpha之后加了一个常数，这使得alpha即使再迭代也不会减小到0.
在降低alpha的函数中，j是迭代次数，i是data下标，这样可以避免参数的严格下降
'''
alpha = 4/(1.0+j+i)+0.0001 #apha decreases with iteration, does not
#改进2：随机选择用来计算的data
randIndex = int(random.uniform(0,len(dataIndex)))#go to 0 because of the constant
h = sigmoid(sum(dataMatrix[randIndex]*weights))
error = classLabels[randIndex] - h
weights = weights + alpha * error * dataMatrix[randIndex]
del(dataIndex[randIndex])
return weights
#用于返回分类结果（0或1）
def classifyVector(inX, weights):
prob = sigmoid(sum(inX*weights))
if prob > 0.5: return 1.0
else: return 0.0
def colicTest():
#训练集和测试集
frTrain = open('horseColicTraining.txt'); frTest = open('horseColicTest.txt')
trainingSet = []; trainingLabels = []
#调用训练集进行训练
for line in frTrain.readlines():
currLine = line.strip().split(' ')
lineArr =[]
#一共有21个特征
for i in range(21):
#lineArr以数组形式保存由字符转换为浮点数的一条数据
lineArr.append(float(currLine[i]))
trainingSet.append(lineArr)
trainingLabels.append(float(currLine[21]))
trainWeights = stocGradAscent1(array(trainingSet), trainingLabels, 1000)
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
#马类问题的测试接口
def multiTest():
#重复执行10次以获取模型分类的平均误差
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)))
dataArr,labelMat=loadDataSet()
weights=gradAscent(dataArr,labelMat)
print("梯度下降法：")
#getA()是numpy的一个函数，作用是将矩阵转成一个ndarray，getA()函数和mat()函数的功能相反，是将一个矩阵转化为数组。
plotBestFit(weights.getA())
print("随机梯度下降法")
weights2=stocGradAscent0(array(dataArr),labelMat)
plotBestFit(weights2)
print("一次随机梯度下降很可能是得不到最优结果的")
print("改进的随机梯度下降法")
weights3=stocGradAscent1(array(dataArr),labelMat)
plotBestFit(weights3)
print("马类问题")
multiTest()