logistic模型matlab代码_机器学习实战个人学习分享（四）：Logistic回归

点击上方蓝字  关注我们 这次的学习是有关Logistic回归的内容，简单的分类，算法不难，但是涉及到算法的优化，其中的梯度上升算法在许多课程里都有涉及，也比较经典，所以这次学习可以系统的归纳，算是基础，老样子，这里用到了Jack Cui 博客和黄海安的视频讲解，内容得到了二位的授权，相关链接如下： Python3《机器学习实战》学习笔记(六)：Logistic回归基础篇之梯度上升算法：https://blog.csdn.net/c406495762/article/details/77723333

Python3《机器学习实战》学习笔记(七)：Logistic回归实战篇之预测病马死亡率：https://blog.csdn.net/c406495762/article/details/77851973

【机器学习实战】【python3版本】【代码讲解】：https://www.bilibili.com/video/BV16t411Q7TM

相关的运行平台和环境

运行平台:  win10

python:  3.7.6

Anaconda:  4.8.3

IDE:  pycharm community

01 前期知识点Logistic回归一般处理的是二分类问题。Logistic回归进行分类的主要思想是：根据现有数据对分类边界线建立回归公式，以此进行分类。其实，Logistic本质上是一个基于条件概率的判别模型(Discriminative Model)。所以要想了解Logistic回归，我们必须先看一看Sigmoid函数 ，我们也可以称它为Logistic函数。它的公式如下：

整合成一个公式，就变成了如下公式：

这个就是sigmoid函数，这个在吴恩达的机器学习里已经学习过，很经典的函数，它的出现主要是可以分类，类似的函数还有tan函数，在MATLAB里也很常见，介绍完sigmoid函数，下面就是一些公式推导，推导出梯度上升算法，之前推导过，这次再推导一下，由于打公式比较累，就借用Jack Cui 博客里的内容，公式推导如下：

这是以x和theta作为参数时，y 为1或者0的概率值，可以将两个式子合并成一个，这个可以看作目标函数，这样合并后，可以看出来如果我获得的数值是y=1，那么我希望Htheta(x)越大越好，这样也就是目标函数越大越好，当y=0，那么我希望1-Htheta(x)越大越好，也就是Htheta(x)越小越好，这是符合逻辑的，这就将两个问题合并起来，来求P的最大值，也就是目标函数的最大值，就可以用最优化理论来求解了。

由于后面涉及到复合求导，所以这里将上式对数化处理，如下：

这是对于一个样本而言的，如果是对于整体而言的话，所以累加即可：

所以问题就转化成求J的最大值问题，这就引出了梯度上升算法。 那么现在开始求解J(θ)对θ的偏导，求解如下( 数学推导 )：

其中：

再由：

可得：

接下来，就剩下第三部分：

综上所述：

因此，梯度上升迭代公式为：

以上就是梯度上升推导的公式，比较简单，主要就是复合链式求导，找出J与theta的关系，然后利用sigmoid函数求导的特殊性，这个推导在BP算法里也有，但是在这里需要注意的是最后的梯度上升公式第二项并没有累加，这个地方累加的代表是批量梯度上升法，没有累加的是随机梯度上升法，在后面跑代码的时候会对比分析。02牛刀初试下面开始跑代码，先来个二次函数理解下梯度上升算法：

def Gradient_Ascent_test():    def f_prime(x_old):  # f(x)的导数        return -2 * x_old + 4    x_old = -1  # 初始值，给一个小于x_new的值    x_new = 0  # 梯度上升算法初始值，即从(0,0)开始    alpha = 0.01  # 步长，也就是学习速率，控制更新的幅度    presision = 0.00000001  # 精度，也就是更新阈值    while abs(x_new - x_old) > presision:        x_old = x_new        x_new = x_old + alpha * f_prime(x_old)  # 上面提到的公式    print(x_new)  # 打印最终求解的极值近似值

这个就是固定的步长，不断地沿着梯度上升的方向去寻找函数最大值，然后更新，直到达到精度要求，我尝试F8不断更新，发现速度还是比较慢的，当我尝试增大步长，发现快了很多，这就是增加了迈步的大小，事实上，太大也会导致振荡，这个都是很常见的问题，废话不多说，运行查看结果：

达到了精度要求，并且接近实际数值。03实战1下面开始针对书上的数据集进行分析，主要就是三个内容，第一是将文本分割后，向量化后画出散点图，第二就是根据梯度上升法求出权重，然后画出决策边界。先看第一部分处理数据：

def loadDataSet():    dataMat = []  # 创建数据列表    labelMat = []  # 创建标签列表    fr = open('testSet.txt')  # 打开文件    for line in fr.readlines():  # 逐行读取        lineArr = line.strip().split()  # 去回车，放入列表        dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])  # 添加数据        labelMat.append(int(lineArr[2]))  # 添加标签    fr.close()  # 关闭文件    return dataMat, labelMat  # 返回

这里都是常规操作，需要注意的是添加了一列1，这是b0，看过吴恩达的机器学习都应该知道，属于偏置项，其他的没什么问题，处理后数据如下：

第二部分就是画散点图

def plotDataSet():    dataMat, labelMat = loadDataSet()  # 加载数据集    dataArr = np.array(dataMat)  # 转换成numpy的array数组    n = np.shape(dataMat)[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])  # 1为正样本        else:            xcord2.append(dataArr[i, 1])            ycord2.append(dataArr[i, 2])  # 0为负样本    fig = plt.figure()    ax = fig.add_subplot(111)  # 添加subplot    ax.scatter(xcord1, ycord1, s=20, c='red', marker='s', alpha=.5)  # 绘制正样本    ax.scatter(xcord2, ycord2, s=20, c='green', alpha=.5)  # 绘制负样本    plt.title('DataSet')  # 绘制title    plt.xlabel('X1')    plt.ylabel('X2')  # 绘制label    plt.show()  # 显示

由于刚才已经切分好成list，直接转化成np的array数组即可，如下：

这就可以直接画散点图了，比较简单，根据标签的不同，将数据分为两类，采取不同的颜色区分，直接给结果：

第三部分就是找决策边界，也是代码的核心部分，贴代码如下：

def gradAscent(dataMatIn, classLabels):    dataMatrix = np.mat(dataMatIn)  # 转换成numpy的mat    labelMat = np.mat(classLabels).transpose()  # 转换成numpy的mat,并进行转置    m, n = np.shape(dataMatrix)  # 返回dataMatrix的大小。m为行数,n为列数。    alpha = 0.001  # 移动步长,也就是学习速率,控制更新的幅度。    maxCycles = 500  # 最大迭代次数    weights = np.ones((n, 1))    for k in range(maxCycles):        h = sigmoid(dataMatrix * weights)  # 梯度上升矢量化公式        error = labelMat - h        weights = weights + alpha * dataMatrix.transpose() * error    return weights.getA()  # 将矩阵转换为数组，返回权重数组

这里先将list转化成np的array方便后续处理，核心的是for循环的内容，将特征向量与权重相乘后，经过sigmoid函数，求出值与实际值(这里是二分类的标签值0和1)差值，然后去更新权重。这里需要注意的是返回的时候将矩阵转换成数组，方便下一步直接取数画出决策边界，求出的权重结果如下：

第三步就是画出决策边界，比较简单，因为权重已经求出直接贴代码：

def plotBestFit(weights):    dataMat, labelMat = loadDataSet()  # 加载数据集    dataArr = np.array(dataMat)  # 转换成numpy的array数组    n = np.shape(dataMat)[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])  # 1为正样本        else:            xcord2.append(dataArr[i, 1])            ycord2.append(dataArr[i, 2])  # 0为负样本    fig = plt.figure()    ax = fig.add_subplot(111)  # 添加subplot    ax.scatter(xcord1, ycord1, s=20, c='red', marker='s', alpha=.5)  # 绘制正样本    ax.scatter(xcord2, ycord2, s=20, c='green', alpha=.5)  # 绘制负样本    x = np.arange(-3.0, 3.0, 0.1)    y = (-weights[0] - weights[1] * x) / weights[2]    ax.plot(x, y)    plt.title('BestFit')  # 绘制title    plt.xlabel('X1')    plt.ylabel('X2')  # 绘制label    plt.show()

运行得到结果如下：

这就完成了第一个实例，结果还不错，下面就是改进了。04实战1改进由于当数据特别大的时候，批量梯度上升法比较慢，计算量太大，所以需要进行优化，可以随机，也可以随机批量，在保证能求出最优权重的前提下减小运算量是最好的，上代码对比，这是全批量梯度上升法：

# 全批量梯度上升法-每次全部使用样本计算，计算量非常大def gradAscent(dataMatIn, classLabels):    dataMatrix = np.mat(dataMatIn)  # 转换成numpy的mat    labelMat = np.mat(classLabels).transpose()  # 转换成numpy的mat,并进行转置    m, n = np.shape(dataMatrix)  # 返回dataMatrix的大小。m为行数,n为列数。    alpha = 0.01  # 移动步长,也就是学习速率,控制更新的幅度。    maxCycles = 500  # 最大迭代次数    weights = np.ones((n, 1))    weights_array = np.array([])    for k in range(maxCycles):        h = sigmoid(dataMatrix * weights)  # 梯度上升矢量化公式        error = labelMat - h        weights = weights + alpha * dataMatrix.transpose() * error        weights_array = np.append(weights_array, weights)    weights_array = weights_array.reshape(maxCycles, n)    return weights.getA(), weights_array  # 将矩阵转换为数组，并返回

下面是随机梯度上升法：

# 随机梯度上升法-每次随机取一个不同的样本进行更新回归系数def stocGradAscent1(dataMatrix, classLabels, numIter=150):    m, n = np.shape(dataMatrix)  # 返回dataMatrix的大小。m为行数,n为列数。    weights = np.ones(n)  # 参数初始化    weights_array = np.array([])  # 存储每次更新的回归系数    for j in range(numIter):        dataIndex = list(range(m))        for i in range(m):            alpha = 4 / (1.0 + j + i) + 0.01  # 降低alpha的大小，每次减小1/(j+i)。            randIndex = int(random.uniform(0, len(dataIndex)))  # 随机选取样本            h = sigmoid(sum(dataMatrix[randIndex] * weights))  # 选择随机选取的一个样本，计算h            error = classLabels[randIndex] - h  # 计算误差            weights = weights + alpha * error * dataMatrix[randIndex]  # 更新回归系数            weights_array = np.append(weights_array, weights, axis=0)  # 添加回归系数到数组中            del (dataIndex[randIndex])  # 删除已经使用的样本    weights_array = weights_array.reshape(numIter * m, n)  # 改变维度    return weights, weights_array  # 返回

对比可知，全批量梯度上升法每次更新需要全部数据跑一遍，而批量梯度上升法每次随机抽取样本进行梯度的计算，然后更新，每次不需要计算全部数据，每次取出的样本会删除，再次抽取，多个局部最优确定整体最优，而且步长会随着逼近不断缩小，也是为了避免振荡，下面画出对比图分析，代码如下：

def plotWeights(weights_array1, weights_array2):    # 设置汉字格式    font = FontProperties(fname=r"c:\windows\fonts\simsun.ttc", size=14)    # 将fig画布分隔成1行1列,不共享x轴和y轴,fig画布的大小为(13,8)    # 当nrow=3,nclos=2时,代表fig画布被分为六个区域,axs[0][0]表示第一行第一列    fig, axs = plt.subplots(nrows=3, ncols=2, sharex=False, sharey=False, figsize=(20, 10))    x1 = np.arange(0, len(weights_array1), 1)    # 绘制w0与迭代次数的关系    axs[0][0].plot(x1, weights_array1[:, 0])    axs0_title_text = axs[0][0].set_title(u'改进的随机梯度上升算法：回归系数与迭代次数关系', FontProperties=font)    axs0_ylabel_text = axs[0][0].set_ylabel(u'W0', FontProperties=font)    plt.setp(axs0_title_text, size=20, weight='bold', color='black')    plt.setp(axs0_ylabel_text, size=20, weight='bold', color='black')    # 绘制w1与迭代次数的关系    axs[1][0].plot(x1, weights_array1[:, 1])    axs1_ylabel_text = axs[1][0].set_ylabel(u'W1', FontProperties=font)    plt.setp(axs1_ylabel_text, size=20, weight='bold', color='black')    # 绘制w2与迭代次数的关系    axs[2][0].plot(x1, weights_array1[:, 2])    axs2_xlabel_text = axs[2][0].set_xlabel(u'迭代次数', FontProperties=font)    axs2_ylabel_text = axs[2][0].set_ylabel(u'W1', FontProperties=font)    plt.setp(axs2_xlabel_text, size=20, weight='bold', color='black')    plt.setp(axs2_ylabel_text, size=20, weight='bold', color='black')    x2 = np.arange(0, len(weights_array2), 1)    # 绘制w0与迭代次数的关系    axs[0][1].plot(x2, weights_array2[:, 0])    axs0_title_text = axs[0][1].set_title(u'梯度上升算法：回归系数与迭代次数关系', FontProperties=font)    axs0_ylabel_text = axs[0][1].set_ylabel(u'W0', FontProperties=font)    plt.setp(axs0_title_text, size=20, weight='bold', color='black')    plt.setp(axs0_ylabel_text, size=20, weight='bold', color='black')    # 绘制w1与迭代次数的关系    axs[1][1].plot(x2, weights_array2[:, 1])    axs1_ylabel_text = axs[1][1].set_ylabel(u'W1', FontProperties=font)    plt.setp(axs1_ylabel_text, size=20, weight='bold', color='black')    # 绘制w2与迭代次数的关系    axs[2][1].plot(x2, weights_array2[:, 2])    axs2_xlabel_text = axs[2][1].set_xlabel(u'迭代次数', FontProperties=font)    axs2_ylabel_text = axs[2][1].set_ylabel(u'W1', FontProperties=font)    plt.setp(axs2_xlabel_text, size=20, weight='bold', color='black')    plt.setp(axs2_ylabel_text, size=20, weight='bold', color='black')    plt.show()

这里画图比较简单，设置汉字格式后，画出权重与迭代次数的关系图，运行看结果：

下面分析对比结果，这里需要解释的是第二个图由于每一次迭代需要计算全部的数据集，而第一个图每次迭代仅仅计算部分的数据集，所以在看图的时候需要考虑进去，可以明显看出，改进后的收敛快，额日期额振荡小，效果很明显。05实例2 下面跑一下预测马的病死率问题，代码简单，贴代码：

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(len(currLine) - 1):            lineArr.append(float(currLine[i]))        trainingSet.append(lineArr)        trainingLabels.append(float(currLine[-1]))    trainWeights = stocGradAscent1(np.array(trainingSet), trainingLabels, 500)  # 使用改进的随即上升梯度训练    errorCount = 0    numTestVec = 0.0    for line in frTest.readlines():        numTestVec += 1.0        currLine = line.strip().split('\t')        lineArr = []        for i in range(len(currLine) - 1):            lineArr.append(float(currLine[i]))        if int(classifyVector(np.array(lineArr), trainWeights)) != int(currLine[-1]):            errorCount += 1    errorRate = (float(errorCount) / numTestVec) * 100  # 错误率计算    print("测试集错误率为: %.2f%%" % errorRate)

368个样本，28个特征，这里存在数据的缺失，处理办法是

• 使用可用特征的均值来填补缺失值；

• 使用特殊值来填补缺失值，如-1；

• 忽略有缺失值的样本；

• 使用相似样本的均值添补缺失值；

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

这里用0来填充，因为当为0时，sigmoid为0.5，中性，不具有任何倾向性。看看代码，这里先用之前的算法，再对比sklearn的，首先将训练集和测试集处理，文本处理多次提了，然后是利用随机梯度上升法，计算出结果后统计准确率，结果如下：

呃呃，准确率不太行啊，换换全批量梯度上升法再试试：

这回好了点，说明在这里随机梯度上升并没有显示出优势。06实例2sklearn测试下面借助sklearn来跑一下代码：

def colicSklearn():    frTrain = open('horseColicTraining.txt')  # 打开训练集    frTest = open('horseColicTest.txt')  # 打开测试集    trainingSet = []    trainingLabels = []    testSet = []    testLabels = []    for line in frTrain.readlines():        currLine = line.strip().split('\t')        lineArr = []        for i in range(len(currLine) - 1):            lineArr.append(float(currLine[i]))        trainingSet.append(lineArr)        trainingLabels.append(float(currLine[-1]))    for line in frTest.readlines():        currLine = line.strip().split('\t')        lineArr = []        for i in range(len(currLine) - 1):            lineArr.append(float(currLine[i]))        testSet.append(lineArr)        testLabels.append(float(currLine[-1]))    classifier = LogisticRegression(solver='sag', max_iter=5000).fit(trainingSet, trainingLabels)    test_accurcy = classifier.score(testSet, testLabels) * 100    print('正确率:%f%%' % test_accurcy)

这里大致分析一下，主要不同的代码是这两行：

 classifier = LogisticRegression(solver='sag', max_iter=5000).fit(trainingSet, trainingLabels)    test_accurcy = classifier.score(testSet, testLabels) * 100

这里有一些参数就简单说一下solver， s olver：优化算法选择参数，只有五个可选参数，即newton-cg,lbfgs,liblinear,sag,saga。默认为liblinear。solver参数决定了我们对逻辑回归损失函数的优化方法，有四种算法可以选择，分别是：

• liblinear：使用了开源的liblinear库实现，内部使用了坐标轴下降法来迭代优化损失函数。

• lbfgs：拟牛顿法的一种，利用损失函数二阶导数矩阵即海森矩阵来迭代优化损失函数。

• newton-cg：也是牛顿法家族的一种，利用损失函数二阶导数矩阵即海森矩阵来迭代优化损失函数。

• sag：即随机平均梯度下降，是梯度下降法的变种，和普通梯度下降法的区别是每次迭代仅仅用一部分的样本来计算梯度，适合于样本数据多的时候。

• saga：线性收敛的随机优化算法的的变重。

这里选择的是sag,，即是随机平均梯度下降法，迭代次数是5000次，然后运行查看结果如下：

结果还可以，试试更换solver为lbfgs，查看结果：

以上就是Logistics回归的主要学习笔记。07总结

ending

点击在看，了解更多精彩内容微信公众号：share_mojie摩羯青年 点击在看，了解更多精彩内容