Logistic回归

写在开头
　　最近在学习一些关于机器学习的基础算法，结合学习Peter Harrington的《机器学习实战》和李航老师的《统计学习方法》两本书以及网上前辈的笔记，写下了以下的学习过程。
　　代码环境：Pytharm/Python3.7
　　内容有参考也有自己的想法，由于自己的理解不足，文章肯定存在很多错误，还恳请各位批评指正。
　　个人理解的回归就是发现变量之间的关系，也就是求回归系数，经常用回归来预测目标值。回归和分类同属于监督学习，所不同的是回归的目标变量必须是连续数值型。
　　
一.关于Logistic
　　今天要学习的logistic回归的主要思想是根据现有的数据对分类边界线建立回归公式，以此进行分类。主要在流行病学中应用较多，比较常用的情形是探索某疾病的危险因素，根据危险因素预测某疾病发生的概率等等。logistic回归的因变量可以是二分类的，也可以是多分类的，但是二分类的更为常用，也更加容易解释，所以实际中最为常用的就是二分类的logistic回归。、，本次的代码也是二分类的。
　　今天我们就二分类进行分析，我们在回归分析中需要一个函数可以接受所有的输入然后预测出类别，假定用0和1分别表示两个类别，logistic函数曲线很像S型，故此我们可以联系sigmoid函数：σ = 1/(1/(1+e-z))。【这个函数，当x足够大时，x>0,值逼近于1，当想x<0，值逼近于0。强烈建议去看看李航老师的《统计学习方法》Logistic这一章的讲解】为了实现logistic回归分类器，我们可以在每个特征上乘以一个回归系数，将所有的乘积相加，将和值代入sigmoid函数中，得到一个范围为0-1之间的数，如果该数值大于0.5则被归入1类，否则被归为0类。
　　基于之前的分析，需要找到回归系数，首先我们可以将sigmoid函数的输入形式记为：z = w0x0 + w1x1 +…+wnxn,其中x为输入数据，相应的w就是我们要求的系数，为了求得最佳系数，结合最优化理论，我们可以选取梯度上升法优化算法。梯度上升法的基本思想是:要找到函数的最大值，最好的方法是沿着该函数的梯度方向寻找。

二.Logistic源代码

# Logistic回归梯度上升优化算法

import numpy


# -------------数据获取及格式化----------------
def loadDataSet():
    dataMat = []
    labelMat = []
    lineArray = []
    # read()和readline()是有区别的
    txt = open("testSet.txt", "r", encoding="utf-8").readlines()  # 打开文件
    for line in txt:
        # strip()
        # 方法用于移除字符串头尾指定的字符（默认为空格或换行符）或字符序列。
        # 注意：该方法只能删除开头或是结尾的字符，不能删除中间部分的字符。
        lines = line.strip('\n')
        lineArray = lines.split()  # lineArray返回的是一个含有三个元素的list
        # 进行向量扩充，简化回归，由此计算出的weight[0]就是对应的b.
        dataMat.append([1.0, float(lineArray[0]), float(lineArray[1])])  # append()能够添加列表作为元素，区分extend（）
        labelMat.append(int(lineArray[2]))
    return dataMat, labelMat


# ------------sigmoid函数模拟了一个阶跃函数-------------
# X大于0为1，小于0为0，比较粗略，只有离0很远的时候才比较准确，导致分类结果不是很准确
def sigmoid(inX):
    import math
    return 1.0 / (1 + math.exp(-inX))


# -----------------梯度上升算法-----------------
def gradAscent(dataMatIn, labelMat):
    dataMatrix = numpy.mat(dataMatIn)  # 转换成矩阵形式
    labelMatrix = numpy.mat(labelMat).transpose()  # transpose()矩阵的转置，将1*100矩阵转置为100*1矩阵
    row, column = numpy.shape(dataMatrix)  # 获取矩阵行数100和列数3
    alpha = 0.001  # 设置步长
    maxCycles = 1000  # 最大迭代次数，迭代次数的增加意味着更小的误差，到一定程度后（error逼近0）效果就不明显了
    weight = numpy.ones((column, 1))  # 生产一个3*1的权重列矩阵，一定要打括号!
    for i in range(maxCycles):
        # 计算真实类别与预测类别的差值，按照差值的方向调整回归系数，使error逼近0
        # 将线性函数转换成概率
        h = numpy.vectorize(sigmoid)(dataMatrix * weight)  # h是100*1列矩阵
        # 计算所得分类与实际分类的误差
        error = labelMatrix - h  # 列矩阵n*1
        weight = weight + alpha * dataMatrix.transpose() * error  # weight是一个3*1的矩阵
    return weight


# ----------------随机梯度上升算法-------------
def stocGradAscent0(dataMatrix, classLabels):
    row, column = numpy.shape(dataMatrix)
    alpha = 0.01  # 设置步长
    weight = numpy.ones(column)  # 只有一个参数则默认为行矩阵
    maxT = numpy.zeros((column, 1))
    for i in range(row):
        a = dataMatrix[i]
        h = sigmoid(sum(dataMatrix[i] * weight))
        error = classLabels[i] - h
        # 不能直接浮点乘以列表。
        for n in range(len(a)):
            a[n] = alpha * error * a[n]
            weight[n] = weight[n] + a[n]
    # for i in range(column):
    #     maxT[i,0] =weight[i]
    return weight


# ----------------改进—随机梯度上升算法-------------
# 改进一：迭代系数随迭代次数增加逐渐减小，抑制波动
# 改进二：通过随机选取样本来更新回归系数。这种方法将减少周期性的波动
def stocGradAscent1(dataMatrix, classLabels):
    import random
    row, column = numpy.shape(dataMatrix)
    weight = numpy.ones(column)  # 只有一个参数则默认为行矩阵
    #  range() 函数返回的是一个可迭代对象（类型是对象），而不是列表类型
    # shuffle() 函数是打乱一个有序列表
    for j in range(100):
        dataIndex = list(range(row))
        random.shuffle(dataIndex)  # 无需赋值，直接函数操作
        for i in range(row):
            alpha = 4 / (1.0 + i + j) + 0.01  # 迭代调整系数
            a = dataMatrix[dataIndex[i]]
            h = sigmoid(sum(dataMatrix[dataIndex[i]] * weight))
            error = classLabels[dataIndex[i]] - h
            # 不能直接浮点乘以列表。
            for n in range(len(a)):
                weight[n] = weight[n] + alpha * error * a[n]
    return weight


# ---------画出数据集和Logistic回归最佳拟合直线的函数------------
def plotBestFit(wei):
    import matplotlib.pyplot as plt
    # getAway（）函数与mat（）函数作用相反，是将一个矩阵转换为一个数组
    # weights = wei.getA()
    weights = wei
    dataMat, labelMat = loadDataSet()
    dataArr = numpy.array(dataMat)  # 强制转换数组类型
    # shape()获取了矩阵的行数及列数，当只想获得行数时增加[0]，只想获得列数时增加[1]
    n = numpy.shape(dataArr)[0]  # 获得行数
    xcode1 = []
    ycode1 = []
    xcode2 = []
    ycode2 = []
    for i in range(n):
        if int(labelMat[i]) == 1:
            # 获取一个点的坐标
            xcode1.append(dataArr[i, 1])
            ycode1.append(dataArr[i, 2])
        else:
            xcode2.append(dataArr[i, 1])
            ycode2.append(dataArr[i, 2])
    fig = plt.figure()  # 创建一个图标
    ax = fig.add_subplot(111)  # 创建一个1*1的图表矩阵，绘制的子图为矩阵中的1号
    # 设置图表参数
    ax.scatter(xcode1, ycode1, s=30, c='red', marker='s')
    ax.scatter(xcode2, ycode2, s=30, c='green')
    # 在[-0.3，0.3）范围内生成一个步长为0.1的ndarray数组
    x = numpy.arange(-3.0, 3.0, 0.1)
    # 由于之前设置了sigmoid函数为0（x=0是分类为0-1的界限）
    # 因此我设定了0=w0*x0+w1*x1+w2*x2
    # 即0=w0+w1*x+w2*y(为什么x0要取1呢？)
    y = (-weights[0] - weights[1] * x) / weights[2]
    ax.plot(x, y)
    plt.xlabel('x1')
    plt.ylabel('y1')
    plt.show()


def main():
    dataMat, labelMat = loadDataSet()
    # weights = gradAscent(dataMat, labelMat) #梯度上升算法
    # weights = stocGradAscent0(dataMat, labelMat)  # 随机梯度上升算法
    weights = stocGradAscent1(dataMat, labelMat)  # 改进-随机梯度上升算法
    print(weights)
    plotBestFit(weights)


if __name__ == '__main__':
    main()

我对算法代码的理解（纯粹个人理解，有错还请多多指教）：
　　将数据与权重值进行相乘，再通过sigmoid函数进行分类，这样就完成了从线性函数到概率的转化，线性函数越接近于正无穷，则概率越接近于1；线性函数越接近于负无穷，则概率越接近于0。
　　那权重值怎么求？初次迭代权重值都相等取1，之后通过求取分类结果与实际结果直接的差值（error）来进行更新权重，error大的自然该权重就大，而更新权重的目的就是让error不断逼近于零（1.应该不会等于0；2.error接近0代表着我通过迭代之后，每一次的分类都是正确的，因此权值就会不便，这个时候就代表完成了分类）{不断的将全部数据与权值相乘更新权值即批处理，迭代次数自己设定}
　　而当数据较多的时候，再使用传统的梯度上升算法就显得比较吃力（因为对于有成千上万个数据或者特征的时候，传统方法的计算复杂度就显著提高，这样会大大减慢数据的处理，甚至无法处理），这样就自然的产生了随机梯度上升，随机梯度梯度上升即一次只用一个样本点来更新回归系数，由于可以在新样本带来时对分类器进行增量式更新，因而随机梯度上升算法是一个在线学习算法。{一次取一个数据与权值相乘更新权值，迭代次数为全部的数据数量}。
　　Logistic回归的目的是寻找一个非线性函数Sigmoid的最佳拟合参数，求解过程可以由最优化
算法来完成。在最优化算法中，最常用的就是梯度上升算法，而梯度上升算法又可以简化为随机
梯度上升算法。
　　随机梯度上升算法与梯度上升算法的效果相当，但占用更少的计算资源。此外，随机梯度上升是一个在线算法，它可以在新数据到来时就完成参数更新，而不需要重新读取整个数据集来进
行批处理运算。

代码中有比较多的地方与《机器学习实战》的源码不一样，有的是为了克服python2与python3之间的语法区别，有的也根据自己的理解进行了一些优化，但是结果都是一样的（主要是理解回归思想感觉通过代码可以更好的理解理论，李航老师的书看起来对我来说有点吃力，但是看懂代码之后又能够进一步理解算法的推导过程，结合起来就有比较好的效果，但还是提醒我要注重理论，学好数学啊！！！）。

三. 估计马疝病的死亡率

import numpy


# 其实可以通过文件相互调用来精简代码的，但是为了清晰还是写在程序里面吧
# ------------sigmoid函数模拟了一个阶跃函数-------------
# X大于0为1，小于0为0，比较粗略，只有离0很远的时候才比较准确，导致分类结果不是很准确
def sigmoid(inX):
    import math
    try:
        ans = 1.0 / 1 + math.exp(-inX)
    except OverflowError:
        if inX > 1000:
            ans = 1.0
        else:
            ans = 0.0
    return ans


def classifyVector(inX, weights):
    a = sum(inX * weights)
    prob = sigmoid(a)
    if prob > 0.5:
        return 1.0
    else:
        return 0.0


# ----------------改进—随机梯度上升算法-------------
# 改进一：迭代系数随迭代次数增加逐渐减小，抑制波动
# 改进二：通过随机选取样本来更新回归系数。这种方法将减少周期性的波动
def stocGradAscent1(dataMatrix, classLabels):
    import random
    row, column = numpy.shape(dataMatrix)
    weight = numpy.ones(column)  # 只有一个参数则默认为行矩阵
    weights = numpy.zeros(column)
    #  range() 函数返回的是一个可迭代对象（类型是对象），而不是列表类型
    # shuffle() 函数是打乱一个有序列表
    for j in range(500):
        dataIndex = list(range(row))
        random.shuffle(dataIndex)  # 无需赋值，直接函数操作
        for i in range(row):
            alpha = 4 / (1.0 + i + j) + 0.01  # 迭代调整系数
            a = dataMatrix[dataIndex[i]]
            b = sum(dataMatrix[dataIndex[i]] * weight)
            h = sigmoid(b)
            error = classLabels[dataIndex[i]] - h
            # 不能直接浮点乘以列表。
            for n in range(len(a)):
                weights[n] = weight[n] + alpha * error * a[n]
    return weights


def coliceTest():
    # -------数据处理------------
    trainingSet = []
    trainingLabels = []
    # readlines()和readline()不一样！！！！！！！！！！
    frTrain = open("horseColicTraining.txt", "r", encoding="utf-8").readlines()
    frTest = open("horseColicTest.txt", "r", encoding="utf-8").readlines()
    for line in frTrain:
        lines = line.strip("\n")
        currLine = lines.split()
        lineArr = []
        # 测试集中有21个特征值+1个标签
        for i in range(21):
            lineArr.append(float(currLine[i]))
        trainingSet.append(lineArr)
        trainingLabels.append(float(currLine[21]))
    # ---------调用算法计算权重-----------
    trainingWeright = stocGradAscent1(numpy.array(trainingSet), trainingLabels)
    errorCount = 0.0
    numTextVec = 0.0
    # -------------测试数据-----------------
    for line in frTest:
        numTextVec += 1.0  # 计算样本数量
        lines = line.strip("\n")
        currLine = lines.split()
        lineArr = []
        # 测试集中有21个特征值+1个标签
        for i in range(21):
            lineArr.append(float(currLine[i]))
        if int(classifyVector(numpy.array(lineArr), trainingWeright)) != int(currLine[21]):
            errorCount += 1
    errorRate = (float(errorCount) / numTextVec)  # 分类错误率
    print("分类错误率为：{}".format(errorRate))
    return errorRate


def main():
    numTests = 10  # 重复测试次数
    errorsum = 0.0
    for i in range(numTests):
        errorsum += coliceTest()
    print("重复测试{}次后，平均错误率为：{:.2f}".format(numTests, errorsum / numTests))


if __name__ == '__main__':
    main()

程序同样根据自己的理解做了修改，期间遇到math range errormath range error，折腾了很久很久解决不了，尝试了限制小数位数，但是没有作用，最后还是利用try解决的。书中提到通过调整步长和迭代次数能够降低分类错误率，但是自己尝试过很多次都不能起到很好的作用，最好的一组数据能够把错误率降低3%，但是远远达不到20%，可能是自己对算法理解还不够深刻，加油吧。

参考文档：
［1］《机器学习实战》第五章
［2］《统计学习方法》李航