机器学习——logistic回归

一，引言

　　假设我们现有一些数据点，我们用一条直线对这些点进行拟合，这个拟合的过程就称作回归。利用logistic回归进行分类的主要思想是：根据现有数据对分类边界线建立回归公式，以此进行分类。

      我们知道，logistic回归主要是进行二分类预测，也即是对于0~1之间的概率值，当概率大于0.5预测为1，小于0.5预测为0.显然，我们不能不提到一个函数，即sigmoid=1/(1+exp(-inX)),该函数的曲线类似于一个s型，在x=0处，函数值为0.5.

于是，为了实现logistic分类器，我们可以在每个特征上都乘以一个回归系数，然后所有的相乘结果进行累加，将这个总结作为输入，输入到sigmoid函数中，从而得到一个大小为0~1之间的值，当该值大于0.5归类为1，否则归类为0，这样就完成了二分类的任务。所以logistic回归可以看成是一种概率估计。

二，基于最优化方法的最佳回归系数确定

　　sigmoid函数的输入记为z，即z=w0x0+w1x1+w2x2+...+wnxn，如果用向量表示即为z=wTx，它表示将这两个数值向量对应元素相乘然后累加起来。其中，向量x是分类器的输入数据，w即为我们要拟合的最佳参数，从而使分类器预测更加准确。也就是说，logistic回归最重要的是要找到最佳的拟合参数。

1 梯度上升法

　　梯度上升法(等同于我们熟知的梯度下降法，前者是寻找最大值，后者寻找最小值)，它的基本思想是：要找到某函数的最大值，最好的方法就是沿着该函数的梯度方向搜寻。如果函数为f，梯度记为D，a为步长，那么梯度上升法的迭代公式为：w：w+a*Dwf(w)。该公式停止的条件是迭代次数达到某个指定值或者算法达到某个允许的误差范围。我们熟知的梯度下降法迭代公式为：w：w-a*Dwf(w)

2 使用梯度上升法寻找最佳参数

　　假设有100个样本点，每个样本有两个特征：x1和x2.此外为方便考虑，我们额外添加一个x0=1，将线性函数z=wTx+b转为z=wTx(此时向量w和x的维度均价1).那么梯度上升法的伪代码如下：

　　初始化每个回归系数为1

　　重复R次：

　　　　计算整个数据集梯度

　　　　使用alpha*gradient更新回归系数的向量

　　返回回归系数

#logistic回归梯度上升优化算
from numpy import *
def loadDataSet():
    dataMat =[]; labelMat =[]
    fr = open('testSet.txt')
    for line in fr.readlines():
        lineArr = line.strip().split() #对当前行去除首尾空行，并按空格进行分离
        #将每一行的两个特征x1,x2,加上x0=1,组成列表并添加到数据集列表中
        dataMat.append([1.0,float(lineArr[0]),float(lineArr[1])])
        labelMat.append(int(lineArr[2])) #将当前行标签添加到标签列表中
    return dataMat,labelMat

def sigmoid(inX):
    return 1/(1+exp(-inX))

#梯度上升法更新最优拟合参数
def gradAscent(dataMatIn,classLabels):
    dataMatrix = mat(dataMatIn) #将数据集列表转为Numpy矩阵
    labelMat = mat(classLabels).transpose()  #将行矩阵转为列矩阵
    m,n = shape(dataMatrix)  #获取矩阵行数和列数
    alpha = 0.001  #学习步长
    maxCycles = 500 #最大迭代次数
    weights = ones((n,1))  #初始化权值参数向量，每个维度均为1.0
    for k in range(maxCycles):  
        h = sigmoid(dataMatrix * weights)  #求当前sigmoid函数预测概率
        error = (labelMat - h)             #计算真实类别和预测类别的差值
        #更新权值参数
        weights = weights + alpha*dataMatrix.transpose()*error
    return weights

画出数据集:

#画出决策边界
def plotBestFit(weights):
    import matplotlib.pyplot as plt
    dataMat,labelMat = loadDataSet()
    dataArr = array(dataMat)
    n = 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 = arange(-3,3,0.1)
    y = (-weights[0]-weights[1]*x)/weights[2]
    ax.plot(x,y)
    plt.xlabel('X1')
    plt.ylabel('X2')
    plt.show()

首先，代码arange(-3,3,0.1)的意思：-3到3间以0.1为步长生成随机数。

• range(start, end, step)，返回一个list对象，起始值为start，终止值为end，但不含终止值，步长为step。只能创建int型list。
• arange(start, end, step)，与range()类似，但是返回一个array对象。需要引入import numpy as np，并且arange可以使用float型数据。

对于最后图片画红线的地即getA()函数（若直接输入weights则会报错）错误原因在于下面这段代码中len(x)=60，而len(y)=1

x = arange(-3,3,0.1)
y = (-weights[0]-weights[1]*x)/weights[2]
ax.plot(x,y)

getA()函数与mat()函数的功能相反，是将一个numpy矩阵转换为数组。

3. 随机梯度上升法

对数据集每个样本

　　计算该样本的梯度

　　使用alpha*gradient更新回顾系数值

 返回回归系数值

#随机梯度上升算法
def stocGradAscent(dataMatrix,classLabels):
    m,n = shape(dataMatrix)
    alpha =0.01
    weights = ones(n)
    for i in range(m):  #循环m次，每次选取数据集的一个样本更新参数
        h = sigmoid(sum(dataMatrix[i]*weights))
        error = classLabels[i]-h
        weights = weights+alpha*error*dataMatrix[i]
    return weights

4，改进的随机梯度上升算法；

#改进的随机梯度上升算法
def atocGradAscent1(dataMatrix,classLabels,numIter=150):
    m,n = shape(dataMatrix)
    weights = ones(n)
    for i in range(numIter):
        #获取数据集行下标列表
        dataIndex = list(range(m))
        for j in range(m):
            alpha = 4/(1+j+i)+0.01
            randIndex = int(random.uniform(0,len(dataIndex))) #获取随机样本
            h = sigmoid(sum(dataMatrix[randIndex]*weights))
            error = classLabels[randIndex]-h
            weights = weights+alpha*error*dataMatrix[randIndex] 
            del(dataIndex[randIndex])  #选取该样本后，将该样本下标删除，确保每次迭代只使用一次
    return weights

注意：3，4算法中传入的第一个参数应是numpy中的数组。否则报错。在4中，用到了删除del操作，所以应将dataIndex数据类型转化为list类型，否则报错（range返回的是range对象）。

三，预测马的死亡率

1. 准备数据：处理数据中的缺失值

我们可能会遇到数据缺失的情况，但有时候数据相当昂贵，扔掉和重新获取均不可取，这显然是会带来更多的成本负担，所以我们可以选取一些有效的方法来解决该类问题。比如：

　　1 使用可用特征的均值填补缺失值

　　2 使用特殊值来填补缺失的特征，如-1

　　3 忽略有缺失值的样本

　　4 使用相似样本的平均值填补缺失值

　　5 使用另外的机器学习算法预测缺失值

　　这里我们根据logstic回归的函数特征，选择实数0来替换所有缺失值，而这恰好能适用logistic回归。因此，它在参数更新时不会影响参数的值。即如果某特征对应值为0 ，那么由公式w：w+alpha*gradient，可知w不会发生改变。

　　此外，由于sigmoid(0)=0.5,表面该特征对结果的预测不具有任何倾向性，因此不会对误差造成影响。

#预测马的死亡率
def classifyVector(inX,weights):
    prob = sigmoid(sum(inX*weights))
    if prob >0.5:
        return 1
    else:
        return 0

def colicTest():
    frTrain = open('horseColicTraining.txt')
    frTest = open('horseColicTest.txt')
    trainingSet =[]; trainingLabels =[]
    for line in frTrain.readlines():
        currLine = line.strip().split('\t')
        lineArr = []
        for i in range(21): #吧0-20个特征加到列表中
            lineArr.append(float(currLine[i]))
        trainingSet.append(lineArr) #把得到的每个列表加到训练集合中
        trainingLabels.append(float(currLine[21])) #把标签加到训练标签中
    #调用算法获得权重的参数值
    trainWeights = stocGradAscent1(array(trainingSet),trainingLabels,500)
    errorCount = 0; numTestVec =0
    for line in frTest.readlines():
        numTestVec +=1
        currLine = line.strip().split('\t')
        lineArr =[]
        for i in range(21):
            lineArr.append(float(currLine[i]))
        #计算分类错误次数，currLine[21]表示真正死亡与否
        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():
    numTests = 10;errorSum =0
    for k in range(numTests):
        errorSum += colicTest()
    print('after %d iterations the average error rate is:%f'
          % (numTests,errorSum/float(numTests)))