Logistic回归

1、什么是回归

已知数据集，求这些数据集的函数表达式的过程

2、logistic回归数据类型：数值型和标称型

3、优点：计算代价不高，易于理解和实现

缺点：容易欠拟合，分类精度可能不高

4、实现原理：将每个特征值乘以一个回归系数，然后将这些值相加，将总和带入到sigmoid函数里面，得到一个0-1的数值，小于0.5的归为0类，大于0.5的归为1类

接下来的任务就是确定最佳的回归系数，我们采用基于最优化方法的最佳回归系数的确定

sigmoid函数：

sigmoid函数的输入记为z

这里的w为回归系数矩阵，x为输入的特征矩阵

我们需要确定的就是w

为了寻找最佳参数，这里我们会用到最优化理论的相关知识

首先我们采用梯度上升的最优化方法：要找到函数的最大值，沿着该函数的梯度方向探寻

函数的梯度表示：

梯度算子总是指向函数增长最快的方向

将梯度算子移动的步长记为a,则梯度算法的迭代公式为：

该公式一直被迭代执行，直到迭代次数达到某个指定值或者算法的精度满足误差范围

5、logistic回归代码实现

伪代码：

初始化回归系数为1

重复n次：

计算数据集的梯度

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

返回回归系数

1>处理文本数据，得到数据集和标签集

#处理文本数据，获取数据集和标签集
def loadDataSet():
    #定义两个数据矩阵和标签矩阵
    dataMat = []; labelMat = []
    fr = open('testSet.txt')
    #遍历每行数据
    for line in fr.readlines():
        #去除首尾空格之后，按照空格分隔
        lineArr = line.strip().split()
        #向矩阵dataMat中添加前两列的数据
        dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])
        #向矩阵labelMat中添加第三列数据
        labelMat.append(int(lineArr[2]))
    return dataMat,labelMat

2>定义sigmiod分类函数

#sigmod函数
def sigmoid(inX):
    return 1.0/(1+exp(-inX))

3>梯度上升算法，找到最好的回归系数w

#梯度上升算法
#计算预测值和真实值的差值，并按照差值的方向调整回归系数
#将loadDataSet返回的数据集和标签集作为参数
def gradAscent(dataMatIn, classLabels):
    #将他们先转化为numpy矩阵形式
    dataMatrix = mat(dataMatIn)             #convert to NumPy matrix
    #对标签矩阵进行转置
    labelMat = mat(classLabels).transpose() #convert to NumPy matrix
    #获取矩阵的行数和列数 m=100 n =2
    m,n = shape(dataMatrix)
    #向目标移动的步长
    alpha = 0.001
    #迭代次数
    maxCycles = 500
    #权重矩阵：3行1列
    weights = ones((n,1))
    for k in range(maxCycles):
        #heavy on matrix operations
        #dataMatrix：100行*3列
        #weights：3行*1列
        #dataMatrix*weights：100行 * 1列
        #h 为100行*1列
        h = sigmoid(dataMatrix*weights)     #matrix mult
        #error 100*1
        error = (labelMat - h)              #vector subtraction
        #dataMatrix.transpose() 3行*100列
        #dataMatrix.transpose()* error 3行1列
        weights = weights + alpha * dataMatrix.transpose()* error #matrix mult
    #返回训练好的回归系数
    return weights

命令行测试得到回归系数矩阵

>>> dataArr,labelMat = logist.loadDataSet()

>>> logist.gradAscent(dataArr,labelMat)
matrix([[ 4.12414349],
        [ 0.48007329],
        [-0.6168482 ]]) 

4>找到了回归系数我们就找到了分隔数据集的函数方程，下面我们的任务就是利用matplotlib画出这条拟合曲线

def plotBestFit(wei):
    import matplotlib.pyplot as plt
    #将自身转为一个多维数组
    #参考博客http://www.jianshu.com/p/37a9cb41401f
    weights = wei.getA()
    dataMat,labelMat=loadDataSet()
    dataArr = array(dataMat)
    #取出数据集的行数
    n = shape(dataArr)[0]
    #xcord1 = []; ycord1 = []用来存储类别标签为1的x,y坐标
    #xcord2 = []; ycord2 = []用来存储类别标签为其他的x,y坐标
    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.0, 3.0, 0.1)
    #由于sigmoid的分类边界是z=0的时候
    #拟合曲线为z = 0 = w0*x0+w1*x1+w2*x2, 故x2 = (-w0*x0-w1*x1)/w2, x0为1,x1为x, x2为y,故有
    y = (-weights[0]-weights[1]*x)/weights[2]
    ax.plot(x, y)
    plt.xlabel('X1'); plt.ylabel('X2');
    plt.show()

>>> reload(logist)
<module 'logist' from 'E:\Python36\logist.py'>
>>> dataArr,labelMat = logist.loadDataSet()
>>> wei = logist.gradAscent(dataArr,labelMat)
>>> logist.plotBestFit(wei)

曲线如下：

5>训练算法

由于每次训练回归系数的时候都要遍历所有的样本计算所有样本的梯度变换，导致计算复杂度偏高

改进：一次仅用一个样本点来更新回归系数（随机梯度上升算法）

伪代码：

初始化所有的回归系数为1

对数据中的每个样本：

计算每个样本的梯度

使用alpha和gradient更新回归系数的值

返回回归系数

代码实现：

#随机梯度上升
def stocGradAscent0(dataMatrix, classLabels):
    m,n = shape(dataMatrix)
    #梯度算子移动步长0.01
    alpha = 0.01
    #初始化回归系数为1
    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

>>> dataArr,labelMat = logist.loadDataSet()
>>> logist.stocGradAscent0(array(dataArr),labelMat)
array([ 1.01702007,  0.85914348, -0.36579921])
>>> weights = logist.stocGradAscent0(array(dataArr),labelMat)
>>> logist.plotBestFit(weights)

效果如下

6>算法改进，增加收敛速度

#改进随机梯度上升算法
#增加了一个参数：迭代次数
def stocGradAscent1(dataMatrix, classLabels, numIter=150):
    m,n = shape(dataMatrix)
    #初始化回归系数
    weights = ones(n)   #initialize to all ones
    #迭代150次
    for j in range(numIter):
        dataIndex = list(range(m))
        #遍历所有样本
        for i in range(m):
            #梯度算子的步长随着迭代次数的增加以及样本变换而变化
            #总体的变化趋势是不断减小，但是不会等于0
            #如果处理的数据是动态变化的，那么可以适当的增大常数项0.0001
            alpha = 4/(1.0+j+i)+0.0001    #apha decreases with iteration, does not
            #random.uniform（）产生从0到len(dataIndex)的一个随机数
            #随机选取样本更新回归系数
            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

第二处是样本的选取是随机选取，并且选取了之后不会再重复选取

通过两处的优化，最后的结果效果也挺好的

7、应用：从疝气病症预测病马的死亡率

1>数据准备：处理数据缺失值

□使用可用特征的均值来填补缺失值；
□使用特殊值来±真补缺失值，如-1;
□ 忽略有缺失值的样本；
□使用相似样本的均值添补缺失值；
□使用另外的机器学习算法预测缺失值。

2>如果特征值缺失，用0来代替

原因如下：

回归系数更新公式：

如果某特征对应值为0，那么该特征的系数将不做更新，即：
w e i g h t s = w e i g h t s

另外，由于sigmoid ( 0 ) = 0 . 5 ，即它对结果的预测不具有任何倾向性，因此上述做法也不会对误差项造成任何影响。基于上述原因，将缺失值用0代替既可以保留现有数据，也不需要对优化算法进行修改。此外，该数据集中的特征取值一般不为0 , 因此在某种意义上说它也满足“特殊值”这个要求。

预处理中做的第二件事是，如果在测试数据集中发现了一条数据的类别标签已经缺失，那么我们的简单做法是将该条数据丢弃。这是因为类别标签与特征不同，很难确定采用某个合适的值来替换。采用logistic回归进行分类时这种做法是合理的

3>开始分类

使 用 logistic回归方法进行分类并不需要做很多工作，所需做的只是把测试集上每个特征向量乘以最优化方法得来的回归系数，再将该乘积结果求和，最后输人到sigmoid 函数中即可0 如果对应sigmoid 值大于0.5就预测类别标签为1 , 否则为0。

#预测疝气病马病死的死亡率

#inX：特征值矩阵   weights：回归系数矩阵
#分类器
def classifyVector(inX, weights):
    #利用sigmoid函数进行分类
    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():
        #去除每行的首尾的空格，然后根据tab进行分隔
        currLine = line.strip().split('	')
        #创建一个行列表
        lineArr =[]
        #依次装取每行的特征值
        #rang(21):0-20
        for i in range(21):
            lineArr.append(float(currLine[i]))
        #读取完一行添加到训练矩阵中
        trainingSet.append(lineArr)
        #将最后一列标签添加到训练标签里面
        trainingLabels.append(float(currLine[21]))
    #使用改进的梯度上升算法进行训练1000次，返回回归系数矩阵
    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

#分类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)))

>>> reload(logist)
<module 'logist' from 'E:\Python36\logist.py'>
>>> logist.multiTest()
E:Python36logist.py:27: RuntimeWarning: overflow encountered in exp
  return 1.0/(1+exp(-inX))
the error rate of this test is: 0.462687
the error rate of this test is: 0.298507
the error rate of this test is: 0.298507
the error rate of this test is: 0.388060
the error rate of this test is: 0.358209
the error rate of this test is: 0.358209
the error rate of this test is: 0.313433
the error rate of this test is: 0.283582
the error rate of this test is: 0.388060
the error rate of this test is: 0.283582
after 10 iterations the average error rate is: 0.343284

下面我们画出分类函数

#训练得到分类器并进行分类
def colicTest():
    #读取训练样本（训练得到分类器）和测试样本（测试分类器的结果和测试样本的标签是否一致）
    frTrain = open('E:/Python36/horseColicTraining.txt');
    frTest = open('E:/Python36/horseColicTest.txt')
    #创建两个矩阵分别装取训练样本的特征值和标签值
    trainingSet = [];
    trainingLabels = []
    #循环读取训练样本的行数
    for line in frTrain.readlines():
        #去除每行的首尾的空格，然后根据tab进行分隔
        currLine = line.strip().split('	')
        #创建一个行列表
        lineArr =[]
        #依次装取每行的特征值
        #rang(21):0-20
        for i in range(21):
            lineArr.append(float(currLine[i]))
        #读取完一行添加到训练矩阵中
        trainingSet.append(lineArr)
        #将最后一列标签添加到训练标签里面
        trainingLabels.append(float(currLine[21]))
    #使用改进的梯度上升算法进行训练1000次，返回回归系数矩阵
    trainWeights = stocGradAscent1(array(trainingSet), trainingLabels, 1000)
    return trainWeights

def multiTest():
    numTests = 10;
    weightssum = 0.0
    for k in range(numTests):
        weightssum += colicTest()
        weightavr = weightssum/float(numTests)
    return weightavr

>>> reload(logist)
<module 'logist' from 'E:\Python36\logist.py'>
>>> weightavr = logist.multiTest()
>>> logist.plotBestFit(weightavr)

发现效果并不理想，以为是迭代次数的原因，但是更改了，效果还是不好，然后以为是回归系数的问题，但是我分类了10次，求取权重的平均值，但是仍然没有效果。如果读者有好的建议，还请分享一下

总结：logistc回归的目的是寻找一个非线性函数sigmoid的最佳拟合参数，求解过程可以由最优化算法来完成。在最优化算法中，最常用的就是梯度上升算法，而梯度上升算法又可以简化为随机梯度上升算法。