机器学习算法（四）逻辑回归理论与python实现+经典示例（从疝气病预测病马的死亡率）

学习笔记更新

什么是逻辑回归？

简要介绍

    逻辑回归是用来解决线性回归问题的，它将线性回归得到的结果通过逻辑函数映射到[0,1]之间，因此称逻辑回归。逻辑回归模型主要用于解决二分类问题，是一个分类模型。

基本问题形式

    前面说的逻辑回归是用来解决线性回归问题的，什么是回归呢？最常见的例子，假设有一些数据点，需要用一条直线对这些数据点进行拟合，求得最佳拟合直线，拟合的过程就是回归。利用逻辑回归对线性问题进行分类的主要思想是：根据现有的数据对分类边界线建立回归公式，将线性回归得到的结果通过逻辑函数映射到 [0，1] 之间，并以此进行分类。

逻辑回归问题求解

逻辑函数

    逻辑函数是一类返回值为逻辑值true或逻辑值false的函数。在二分类问题中，我们希望输入数据特征之后，函数输出0或1，表示两种不同的类别，同时我们希望函数在0到1之间有一个变化过程，而不是如赫维赛德阶跃函数那样在跳跃点上瞬间到 1，没有一个变化过程，函数不连续。Sigmod函数满足上述的需求，坐标轴范围小时，可以看到函数值从0到1是有一个慢慢上升的过程的，同时，坐标轴范围大一些，可以观察到函数值在0处发生明显变化。

sigmod 计算式为

因此，问题可以这样求解：

    将现有待分类数据点的特征（X0，X1，X2……），每个特征乘以一个回归系数，然后将所 有结果值相加，得到值z，

将z带入sigmod函数中，会得到一个在0~1之内的数字。将>0.5的值划分为第一类，小于0.5的值划分为第二类。
    如果有一个样本被误分，本属于第一类，但分类结果是第二类，我们可以知道是因为z值输入所导致的。此时，需要通过调整回归系数w，进而调整z值，使其能够尽快的划分到正确的类别中。

    想法固然美好，如何实现呢？学过梯度下降法的话不难想到，上述调整z值的任务可以用梯度下降法来实现，学习传送 -->（传送至梯度下降法）。函数值在每一点沿着梯度的方向走，如果函数是凸函数，便可找到函数的最优解，如果函数是非凸函数，求得的解可能为全局最优解也可能是局部最优解。线性函数是凸函数的特例，梯度下降法非常适用。

梯度上升与梯度下降本质一个道理，只是方向相反。梯度上升法迭代公式如下：

一直执行，直到达到设定的迭代次数或算法达到允许的误差的误差范围内。
梯度上升法的目的：（调整参数w1，w2，使得计算得到的z值输入sigmod后能将大多数样本正确分类）。

算法步骤

1、输入数据
2、初始化：参数w1，w2……=1，步长α=0.01。
3、设定梯度下降迭代次数，进行梯度下降，计算偏导，更新w。 (训练模型)
4、返回最终w，得到决策边界。
5、w值确定后，输入测试数据集，计算分类错误率。 (测试模型)

程序实现

加载数据

######加载数据###########
def loadDataSet():
    dataMat=[];datalabel=[]   #设置list，存储数据和标签
    fr=open('testSet.txt')     
    for line in fr.readlines():  
        lineClean=line.strip().split()  #对每一行数据进行处理，去掉前后空白字符，以空字符分割，空字符包括空格、制表符、回车符、换行符等
        dataMat.append([1,float(lineClean[0]),float(lineClean[1])])
        datalabel.append(int(lineClean[2]))
    return dataMat,datalabel

梯度上升

###################梯度上升#################
def gradAcent(dataMat,dataLabel):
    dataMatrix=numpy.mat(dataMat)
    dataLabel_Matrix=mat(dataLabel).transpose()  #转换为numpy矩阵，便于进行矩阵操作
    m,n=shape(dataMatrix)
    #print(n)
    alpha=0.001   #初始化步长
    maxCycles=500  #设置最大迭代次数
    weights=numpy.ones((n,1)) # #创建nX1的矩阵，矩阵元素全为1
    for k in range(maxCycles):     #进行参数更新，循环次数为maxCycles
        h=sigmod(dataMatrix*weights)
        error=dataLabel_Matrix-h
        weights=weights+alpha*dataMatrix.transpose()*error  #梯度上升，损失函数对各变量求偏导后得到w=w+α*x*(y-h)   #y为真实值
    return weights   #返回更新后的参数
   
################定义sigmod函数###############
def sigmod(inX):
    return 1.0/(1+numpy.exp(-inX))   

####执行上述函数
dataArr,dataLabel=loadDataSet()
weights=gradAcent(dataArr,dataLabel)
print(weights)   #输出更新后的参数

程序输出：

[[ 4.12414349]
 [ 0.48007329]
 [-0.6168482 ]]

可视化

##画出决策边界

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.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()

###执行上述函数
dataArr,dataLabel=loadDataSet()
weights=gradAcent(dataArr,dataLabel)
plotBestFit(weights.getA())  ##通过getA()这个方法可以将weights返回成一个数组对象

程序运行结果如下：

随机梯度上升法

一次仅用一个样本点来更新回归系数

###############3随机梯度上升法#########################
def stocGradAscent(dataArray,dataLabel):
    dataArray=array(dataArray)
    m,n=shape(dataArray)
    alpha=0.01
    weights=ones(n)  #返回一个数组，[1,1,1]
    for i in range(m):
        h=sigmod(sum(dataArray[i]*weights))
        error=dataLabel[i]-h
        weights=weights + alpha * dataArray[i] * error
    return weights


######################执行函数######################
dataArr,dataLabel=loadDataSet()
weights=stocGradAscent(dataArr,dataLabel)
print(weights)
plotBestFit(weights)

程序输出：

[ 1.01702007  0.85914348 -0.36579921]    #更新后的参数值

    这里，随机梯度上升法分类效果直观上不如原梯度上升法，但不能这样直接比较。这样对第二种方法是不公平的，因为原始的梯度上升，是在整个数据集上迭代了500次才得到的结果，而随机梯度上升只遍历了一次数据集。有研究表明，将上述随机梯度上升算法在整个数据集上运行200次，特征X2的系数只经历50次迭代就达到了稳定值。而X0与X1的系数则需要更多的迭代，并且在大波动停止后，会有小范围的周期性波动。这里不难理解，参数上上下下在一定范围内波动，是因为有一些有一些样本点不能被正确分类，数据集不是完全线性可分的，每次更新参数，对这些特殊的点来说总是不能将他们正确分类。因此每次迭代都会引起系数较大变化。针对随机梯度上升算法种存在的问题，又产生了改进的随机梯度上升法，程序如下：（只想了解 解决逻辑回归问题主要步骤的伙伴，如果对这一部分不太感兴趣的话可以跳过，因为这里并不影响对整个问题的认识）

改进的随机梯度上升法

##改进的随机梯度上升法：
def Pro_stocGradAscent(dataArray,dataLabel,numIter=200):  
    dataArray=array(dataArray)
    m,n=shape(dataArray)
    alpha=0.01
    weights=ones(n)
    for j in range(numIter):
        dataIndex=list(range(m))
        for i in range(m):
            alpha=4/(1.0+j+i)+0.01     #每一次都会更新步长，随着迭代次数不断减小，但不会为0。目的：保证多次迭代后新参数仍具有影响
            randIndex=int(random.uniform(0,len(dataIndex)))  #随机选取一个来进行更新
            h=sigmod(sum(dataArray[randIndex]*weights))
            # print(type(h))
            # print(type(dataLabel[randIndex]))
            # print(dataLabel[randIndex])
            error=dataLabel[randIndex]-h
            weights=weights+alpha * dataArray[randIndex] * error
            del(dataIndex[randIndex])
    return weights

######################执行函数######################
dataArr,dataLabel=loadDataSet()
weights=stocGradAscent(dataArr,dataLabel)
print(weights)
plotBestFit(weights)

程序输出结果

[15.05877546  0.71607795 -1.97753643]    #更新后参数值

    改进后的随机梯度上升法可以解决随机梯度上升中存在的参数周期性波动的问题，因为改进后的方法里并不是对所有样本按顺序逐次遍历，依次更新参数；而是随机挑选样本来进行更新。二者虽然最后都遍历了所有样本，但是遍历的顺序不同，因此也避免了特殊点引起的周期性波动。同时改进后的随机梯度上升可以使参数收敛的更快。读者有兴趣的话，可以对Pro_stocGradAscent() 中的三个参数进行调整，以达到更好的效果。

    到目前为止，重点一直在分析如何更新回归系数，而对于具体的问题如何求解，仍是一个较为模糊的概念。通过下面这个示例，可以清晰地了解从数据输入到最后给出分类结果整个过程。

经典示例：从疝气病症预测病马的死亡率

    使用逻辑回归来预测患有疝病的马的存活问题。数据中包含299个样本和20个特征。数据集中的特征包含了医院检的一些指标，有的比较主观，有的难以测量，另外存在数据缺失的问题，不过不用担心，对于数据预处理，已经有前辈处理好了，可以在此基础上来学习。
数据集：包含20个特征，最后一列是标签

问题求解：

    不需要太多繁杂的步骤，所需要做的就是将测试集上每个特征向量乘以最优化方法得来的回归系数，再将该乘积结果求和，（这个过程就是文章开头说起的线性回归得到的参数后，带入特征，求得的线性回归结果z），最后输入到sigmod函数中即可。如果对应sigmod值>0.5,则预测为标签类别1（死亡率高），否则为0（死亡率低）。
程序如下：

def colicTest():
    trainArray=[];testArray=[]       #构建list,存储从文件中读取出来的训练集，测试集
    train_label=[];test_label=[]
    trainfile=open('horseColicTraining.txt')
    testfile=open('horseColicTest.txt')
    for line in trainfile.readlines():      #逐行读取训练集
        curentLine=line.strip().split()
        lineArr=[]
        for i in range(len(curentLine)):
            lineArr.append(float(curentLine[i]))  #存储训练集样本特征
        trainArray.append(lineArr)   #存储训练集标签
        train_label.append(float(curentLine[21]))   ####转为float，不然会当做字符串处理

    weights=Pro_stocGradAscent(array(trainArray),train_label,1000)   #使用改进的随机梯度上升法求参数
    numTest=0.0;errorCount=0.0                          #numTest：统计测试集样本数量，errorCount：统计分类错误的样本数量
    for line in testfile.readlines():   #逐行读取测试集
        numTest+=1
        currentLine=line.strip().split()
        lineArr=[]
        for i in range(len(curentLine)):
            lineArr.append(float(curentLine[i]))
        if classfiy(lineArr,weights)!=int(currentLine[21]):   #将分类结果与测试集样本真实标签对比
            errorCount+=1
    
    errorRate=float(errorCount)/numTest
    print("the error rate of this test is: %f" %errorCount)   #输出错误率
    return errorRate


#####执行函数
colicTest()

程序输出结果

the error rate of this test is: 0.402985

进行多次测试求平均错误率

def muliTest():
    numTest=10
    errorSum=0.0
    for i in range(numTest):
        errorSum+=colicTest()
    print("After %d iterations the average error rate is: %f" %(numTest,float(errorSum/numTest)))

####执行函数
muliTest()

程序输出：

the error rate of this test is: 0.343284
the error rate of this test is: 0.328358
the error rate of this test is: 0.283582
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.402985
the error rate of this test is: 0.343284
the error rate of this test is: 0.373134
the error rate of this test is: 0.268657
the error rate of this test is: 0.373134
After 10 iterations the average error rate is: 0.343284

总结

参考文章：逻辑回归的常见面试点总结。 写的非常好，感谢博主分享。

优点：

1、形式简单，模型的可解释性非常好。从特征的权重可以看到不同的特征对最后结果的影响，某个特征的权重值比较高，则这个特征最后对结果的影响会比较大。
2、模型效果不错。在工程上是可以接受的，如果特征工程做的好，效果不会太差，并且特征工程可以大家并行开发，大大加快开发的速度。
3、训练速度较快。分类的时候，计算量仅仅只和特征的数目相关。并且逻辑回归的分布式优化sgd发展比较成熟，训练的速度可以通过堆机器进一步提高，这样我们可以在短时间内迭代好几个版本的模型。
4、资源占用小,尤其是内存。因为只需要存储各个维度的特征值，。
5、方便输出结果调整。逻辑回归可以很方便的得到最后的分类结果，因为输出的是每个样本的概率分数，我们可以很容易的对这些概率分数进行划分。

缺点

1、准确率并不是很高。因为形式非常的简单(非常类似线性模型)，很难去拟合数据的真实分布。
2、很难处理数据不平衡的问题。举个例子：如果我们对于一个正负样本非常不平衡的问题比如正负样本比 10000:1.我们把所有样本都预测为正也能使损失函数的值比较小。但是作为一个分类器，它对正负样本的区分能力不会很好。
3、处理非线性数据较麻烦。逻辑回归在不引入其他方法的情况下，只能处理线性可分的数据，或者进一步说，处理二分类的问题 。
4、逻辑回归本身无法筛选特征。有时候，我们会用gbdt来筛选特征，然后再上逻辑回归。

学海无涯，个人整理，内容难免会有纰漏，欢迎道友指正，感激不尽！