5-1节 Logistic回归|使用 Logistic 回归在简单数据集上的分类项目汇总|机器学习实战-学习笔记

文章原创,最近更新:2018-09-19

前言:
本文介绍机器学习分类算法中的Logistic回归算法并给出伪代码,Python代码实现。

学习参考链接:
1.第5章 Logistic回归

本章节的主要内容是:
项目案例1: 使用 Logistic 回归在简单数据集上的分类。

1.Logistic回归项目案例介绍:

项目案例1:

使用 Logistic 回归在简单数据集上的分类

项目概述:

在一个简单的数据集上,采用梯度上升法找到 Logistic 回归分类器在此数据集上的最佳回归系数.

Logistic 回归 工作原理
每个回归系数初始化为 1
重复 R 次:
    计算整个数据集的梯度
    使用 步长 x 梯度 更新回归系数的向量
返回回归系数
开发流程:
  • 收集数据: 可以使用任何方法
  • 准备数据: 由于需要进行距离计算,因此要求数据类型为数值型。另外,结构化数据格式则最佳
  • 分析数据: 画出决策边界
  • 训练算法: 使用梯度上升找到最佳参数
  • 测试算法: 使用 Logistic 回归进行分类
  • 使用算法: 对简单数据集中数据进行分类
Logistic 回归 算法特点:
  • 优点: 计算代价不高,易于理解和实现。
  • 缺点: 容易欠拟合,分类精度可能不高。
    适用数据类型: 数值型和标称型数据。
数据集介绍

我们采用存储在 TestSet.txt 文本文件中的数据,存储格式如下:

-0.017612   14.053064   0
-1.395634   4.662541    1
-0.752157   6.538620    0
-1.322371   7.152853    0
0.423363    11.054677   0

2.改善前-项目相关代码

2.1 loadDataSet()

这段代码主要是解析数据

现在需要将截图1中的数据集,如何转换成截图2两组列表数据集?其中列表1的数据类型是float,列表2的数据类型是int.


截图1

截图2

具体方法见代码如下:

def loadDataSet(file_name):
    '''
    Desc: 
        加载并解析数据
    Args:
        file_name -- 要解析的文件路径
    Returns:
        dataMat -- 原始数据的特征
        labelMat -- 原始数据的标签,也就是每条样本对应的类别。即目标向量
    '''
    # dataMat为原始数据, labelMat为原始数据的标签
    dataMat = []
    labelMat = []
    fr = open(file_name)
    for line in fr.readlines():
        lineArr = line.strip().split()
        # 为了方便计算,我们将 X0 的值设为 1.0 ,也就是在每一行的开头添加一个 1.0 作为 X0
        dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])
        labelMat.append(int(lineArr[2]))
    return dataMat, labelMat

测试代码及其结果如下:

filename="testSet.txt"
dataMat,labelMat=loadDataSet(filename)

print(dataMat)
[[1.0, -0.017612, 14.053064], [1.0, -1.395634, 4.662541], [1.0, -0.752157, 6.53862], [1.0, -1.322371, 7.152853], [1.0, 0.423363, 11.054677], [1.0, 0.406704, 7.067335], [1.0, 0.667394, 12.741452], [1.0, -2.46015, 6.866805], [1.0, 0.569411, 9.548755], [1.0, -0.026632, 10.427743], [1.0, 0.850433, 6.920334], [1.0, 1.347183, 13.1755], [1.0, 1.176813, 3.16702], [1.0, -1.781871, 9.097953], [1.0, -0.566606, 5.749003], [1.0, 0.931635, 1.589505], [1.0, -0.024205, 6.151823], [1.0, -0.036453, 2.690988], [1.0, -0.196949, 0.444165], [1.0, 1.014459, 5.754399], [1.0, 1.985298, 3.230619], [1.0, -1.693453, -0.55754], [1.0, -0.576525, 11.778922], [1.0, -0.346811, -1.67873], [1.0, -2.124484, 2.672471], [1.0, 1.217916, 9.597015], [1.0, -0.733928, 9.098687], [1.0, -3.642001, -1.618087], [1.0, 0.315985, 3.523953], [1.0, 1.416614, 9.619232], [1.0, -0.386323, 3.989286], [1.0, 0.556921, 8.294984], [1.0, 1.224863, 11.58736], [1.0, -1.347803, -2.406051], [1.0, 1.196604, 4.951851], [1.0, 0.275221, 9.543647], [1.0, 0.470575, 9.332488], [1.0, -1.889567, 9.542662], [1.0, -1.527893, 12.150579], [1.0, -1.185247, 11.309318], [1.0, -0.445678, 3.297303], [1.0, 1.042222, 6.105155], [1.0, -0.618787, 10.320986], [1.0, 1.152083, 0.548467], [1.0, 0.828534, 2.676045], [1.0, -1.237728, 10.549033], [1.0, -0.683565, -2.166125], [1.0, 0.229456, 5.921938], [1.0, -0.959885, 11.555336], [1.0, 0.492911, 10.993324], [1.0, 0.184992, 8.721488], [1.0, -0.355715, 10.325976], [1.0, -0.397822, 8.058397], [1.0, 0.824839, 13.730343], [1.0, 1.507278, 5.027866], [1.0, 0.099671, 6.835839], [1.0, -0.344008, 10.717485], [1.0, 1.785928, 7.718645], [1.0, -0.918801, 11.560217], [1.0, -0.364009, 4.7473], [1.0, -0.841722, 4.119083], [1.0, 0.490426, 1.960539], [1.0, -0.007194, 9.075792], [1.0, 0.356107, 12.447863], [1.0, 0.342578, 12.281162], [1.0, -0.810823, -1.466018], [1.0, 2.530777, 6.476801], [1.0, 1.296683, 11.607559], [1.0, 0.475487, 12.040035], [1.0, -0.783277, 11.009725], [1.0, 0.074798, 11.02365], [1.0, -1.337472, 0.468339], [1.0, -0.102781, 13.763651], [1.0, -0.147324, 2.874846], [1.0, 0.518389, 9.887035], [1.0, 1.015399, 7.571882], [1.0, -1.658086, -0.027255], [1.0, 1.319944, 2.171228], [1.0, 2.056216, 5.019981], [1.0, -0.851633, 4.375691], [1.0, -1.510047, 6.061992], [1.0, -1.076637, -3.181888], [1.0, 1.821096, 10.28399], [1.0, 3.01015, 8.401766], [1.0, -1.099458, 1.688274], [1.0, -0.834872, -1.733869], [1.0, -0.846637, 3.849075], [1.0, 1.400102, 12.628781], [1.0, 1.752842, 5.468166], [1.0, 0.078557, 0.059736], [1.0, 0.089392, -0.7153], [1.0, 1.825662, 12.693808], [1.0, 0.197445, 9.744638], [1.0, 0.126117, 0.922311], [1.0, -0.679797, 1.22053], [1.0, 0.677983, 2.556666], [1.0, 0.761349, 10.693862], [1.0, -2.168791, 0.143632], [1.0, 1.38861, 9.341997], [1.0, 0.317029, 14.739025]]

print(labelMat)
[0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0]

2.2 sigmoid()

这个函数主要是定义sigmoid阶跃函数.

什么叫Sigmoid 函数呢?

首先我们先来看看这个函数的具体计算公式,如下:



下图给出了 Sigmoid 函数在不同坐标尺度下的两条曲线图。当 x 为 0 时,Sigmoid 函数值为 0.5 。随着 x 的增大,对应的 Sigmoid 值将逼近于 1 ; 而随着 x 的减小, Sigmoid 值将逼近于 0 。如果横坐标刻度足够大, Sigmoid 函数看起来很像一个阶跃函数。



那么如何用程序代码实现这个函数呢?

下面来看看具体的方法,如下:

# sigmoid阶跃函数
def sigmoid(inX):
    # return 1.0 / (1 + exp(-inX))
    return 1.0/(1+np.exp(-inX)) 

涉及到相关知识点
知识点1:
Python exp() 函数描述: exp() 函数返回 x 的指数,e^{x}

1)语法

import math
math.exp(x)

注意:exp() 是不能直接访问的,需导入 math 模块,通过静态对象调用该方法。

  • 参数
    x 数值表达式
  • 返回值
    返回 x 的指数 e^{x}

2)案例:

# -*- coding: UTF-8 -*-
import numpy as np   # 导入 math 模块

print ("np.exp(2) : ", np.exp(2))
print ("np.exp(100.12) : ", np.exp(100.12))
print ("np.exp(100.72) : ", np.exp(100.72))
print( "np.exp(119L) : ", np.exp(119L))

输出结果如下:

np.exp(2) :  7.38905609893
np.exp(100.12) :  3.03084361407e+43
np.exp(100.72) :  5.52255713025e+43
np.exp(119L) :  4.7978133273e+51

2.3 gradAscent()

这个函数主要是定义Logistic 回归梯度上升优化算法

具体方法见代码如下:

def gradAscent(dataMatIn, classLabels):
    #用mat函数将列表转换成numpy矩阵
    dataMatrix = np.mat(dataMatIn)
    #用mat函数将列表转换成numpy矩阵,并进行转置
    labelMat = np.mat(classLabels).transpose()
    #返回dataMatrix的大小。m为行数,n为列数。
    m,n = np.shape(dataMatrix)
    #移动步长,也就是学习速率,控制更新的幅度。
    alpha = 0.001
    #最大迭代次数
    maxCycles = 500
    # weights 代表回归系数, 此处的 ones((n,1)) 创建一个长度和特征数相同的矩阵,其中的数全部都是 1
    weights = np.ones((n,1))
    for k in range(maxCycles):
        # m*3 的矩阵 * 3*1 的矩阵 = m*1的矩阵
        h = sigmoid(dataMatrix*weights)
        # labelMat是实际值,labelMat - 是两个向量相减
        error = (labelMat - h) 
        # 0.001* (3*m)*(m*1) 表示在每一个列上的一个误差情况,最后得出 x1,x2,xn的系数的偏移量
         # 矩阵乘法,最后得到回归系数
        weights = weights + alpha * dataMatrix.transpose() * error
    return np.array(weights)

测试代码及其结果如下:

filename="testSet.txt"
dataMat,labelMat=loadDataSet(filename)
a=gradAscent(dataMat,labelMat)
print(a)

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

相关知识点

  • 知识点1:这个代码weights = weights + alpha * dataMatrix.transpose() * error为什么这样写呢?

一起看一下证明的过程,具体如下:
1)Sigmoid函数求导


这里z=w^{t}x+b

2)迭代更新公式

proof:
这里我们先规定logsitic函数的形式:


然后分析分类问题,y有两种取值0和1,根据这个概率p ( y=1 | x,w ) = h(x) , 所以p ( y=0 | x,w ) = 1-h(x),我们把两个概率结合起来,得到:



现在我们可以开始求Logistic回归系数了,因为logistic回归可以看做是一种概率模型,y发生的概率与w有关,因而我们可以对w使用极大似然估计(极大似然估计的知识推荐看《统计学完全教程》,讲的比较全面),使y发生的概率最大,此时的w便是我们寻找的最优回归参数,现根据概率密度写出极大似然函数:


取对数得到对数极大似然函数:

现在要使用梯度上升法w(即迭代更新公式),则我们需要求出对数极大似然函数的梯度:

这里yi对应书中labelMat的标签值,h(xi)对应sigmoid函数的输出值,将梯度带入公式:

这就得到了书中梯度上升法的最后一步迭代,这里alpha是步长,决定梯度上升的快慢。
学习参考链接:Logistic回归-数学原理(1)机器学习实战

  • 知识点2:矩阵乘法

小明今天要做饭,消耗2斤肉,1斤蔬菜。肉每斤20元,蔬菜每斤5元,则一共需多少花费?

这个问题的答案很简单:

我们用向量相乘的方法写出来:

如果小明第二天有另一种做饭的方法,需要消耗1斤肉,4斤蔬菜,那么这两种方法的花费各是多少呢?我们显然需要另算这第二种方法的花费。把这个做饭方式写在第二个矩阵(向量是宽度或长度为1的矩阵)里:


小明家附近还有另一个菜市场,那里肉每斤15元,蔬菜每斤10元。那么,小明如果去这个菜市场,花费又是多少呢(分别计算上述两种做饭方式)?我们把这另外的一种价格写进第一个矩阵里:

这样我们看到了一个矩阵乘法的例子。在左边的这个矩阵的每一行,都代表了一种价目表;在右边的矩阵的每一列,都代表了一种做饭方式。那么所有可能的组合所最终产生的花费,则在结果矩阵中表示出来了。
学习参考链接:矩阵乘法的本质是什么?

2.4 plotBestFit()

这个函数主要是画出数据集和 Logistic 回归最佳拟合直线的函数

# 画出数据集和 Logistic 回归最佳拟合直线的函数
def plotBestFit(dataArr, labelMat, weights):
    n = np.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 = np.arange(-3.0, 3.0, 0.1)
    y = (-weights[0]-weights[1]*x)/weights[2]
    ax.plot(x, y)
    plt.xlabel('X'); plt.ylabel('Y')
    plt.show()

测试代码及其结果如下:

filename="testSet.txt"
dataArr,labelMat=loadDataSet(filename)
weights=gradAscent(dataArr,labelMat)
plotBestFit(dataArr, labelMat, weights)

相关知识点:

  • 知识点1:
y = (-weights[0]-weights[1]*x)/weights[2]这段代码的由来?

首先理论上是这个样子的。

dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])
w0*x0+w1*x1+w2*x2=f(x)

x0最开始就设置为1, x2就是我们画图的y值,而f(x)被我们磨合误差给算到w0,w1,w2身上去了,所以:

 w0+w1*x+w2*y=0 => y = (-w0-w1*x)/w2   

2.5 testLR()

这段代码主要是测试算法: 使用 Logistic 回归进行分类

备注:这一段相当于2.4plotBestFit()的测试代码,只不过简化成了函数.

def testLR():
    dataMat, labelMat = loadDataSet("TestSet.txt")
    dataArr = np.array(dataMat)
    weights = gradAscent(dataArr, labelMat)
    plotBestFit(dataArr, labelMat, weights)

测试代码及其结果如下:

testLR()

2.6 完整的代码如下:

import matplotlib.pyplot as plt
import numpy as np

#解析数据
def loadDataSet(file_name):
    '''
    Desc: 
        加载并解析数据
    Args:
        file_name -- 要解析的文件路径
    Returns:
        dataMat -- 原始数据的特征
        labelMat -- 原始数据的标签,也就是每条样本对应的类别。即目标向量
    '''
    # dataMat为原始数据, labelMat为原始数据的标签
    dataMat = []
    labelMat = []
    fr = open(file_name)
    for line in fr.readlines():
        lineArr = line.strip().split()
        # 为了方便计算,我们将 X0 的值设为 1.0 ,也就是在每一行的开头添加一个 1.0 作为 X0
        dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])
        labelMat.append(int(lineArr[2]))
    return dataMat, labelMat

# sigmoid阶跃函数
def sigmoid(inX):
    # return 1.0 / (1 + exp(-inX))
    return 1.0/(1+np.exp(-inX)) 

# Logistic 回归梯度上升优化算法
def gradAscent(dataMatIn, classLabels):
    #用mat函数将列表转换成numpy矩阵
    dataMatrix = np.mat(dataMatIn)
    #用mat函数将列表转换成numpy矩阵,并进行转置
    labelMat = np.mat(classLabels).transpose()
    #返回dataMatrix的大小。m为行数,n为列数。
    m,n = np.shape(dataMatrix)
    #移动步长,也就是学习速率,控制更新的幅度。
    alpha = 0.001
    #最大迭代次数
    maxCycles = 500
    # weights 代表回归系数, 此处的 ones((n,1)) 创建一个长度和特征数相同的矩阵,其中的数全部都是 1
    weights = np.ones((n,1))
    for k in range(maxCycles):
        # m*3 的矩阵 * 3*1 的矩阵 = m*1的矩阵
        h = sigmoid(dataMatrix*weights)
        # labelMat是实际值,labelMat - 是两个向量相减
        error = (labelMat - h) 
        # 0.001* (3*m)*(m*1) 表示在每一个列上的一个误差情况,最后得出 x1,x2,xn的系数的偏移量
         # 矩阵乘法,最后得到回归系数
        weights = weights + alpha * dataMatrix.transpose() * error
    return np.array(weights)


# 画出数据集和 Logistic 回归最佳拟合直线的函数
def plotBestFit(dataArr, labelMat, weights):
    n = np.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 = np.arange(-3.0, 3.0, 0.1)
    y = (-weights[0]-weights[1]*x)/weights[2]
    ax.plot(x, y)
    plt.xlabel('X'); plt.ylabel('Y')
    plt.show()

#测试算法: 使用 Logistic 回归进行分类
def testLR():
    dataMat, labelMat = loadDataSet("TestSet.txt")
    dataArr = np.array(dataMat)
    weights = gradAscent(dataArr, labelMat)
    plotBestFit(dataArr, labelMat, weights)


testLR()

3.改善后-项目相关代码

梯度上升算法在每次更新回归系数时都需要遍历整个数据集,该方法在处理 100 个左右的数据集时尚可,但如果有数十亿样本和成千上万的特征,那么该方法的计算复杂度就太高了。一种改进方法是一次仅用一个样本点来更新回归系数,该方法称为 随机梯度上升算法。由于可以在新样本到来时对分类器进行增量式更新,因而随机梯度上升算法是一个在线学习(online learning)算法。与 “在线学习” 相对应,一次处理所有数据被称作是 “批处理” (batch) 。

随机梯度上升算法可以写成如下的伪代码:

#随机梯度上升算法-改善前
每个回归系数初始化为 1
重复 R 次:
    计算整个数据集的梯度
    使用 步长 x 梯度 更新回归系数的向量
返回回归系数


#随机梯度上升算法-改善后

所有回归系数初始化为 1
对数据集中每个样本
    计算该样本的梯度
    使用 alpha x gradient 更新回归系数值
返回回归系数值

以下是随机梯度上升算法的实现代码:
改善前

# Logistic 回归梯度上升优化算法
def gradAscent(dataMatIn, classLabels):
    #用mat函数将列表转换成numpy矩阵
    dataMatrix = np.mat(dataMatIn)
    #用mat函数将列表转换成numpy矩阵,并进行转置
    labelMat = np.mat(classLabels).transpose()
    #返回dataMatrix的大小。m为行数,n为列数。
    m,n = np.shape(dataMatrix)
    #移动步长,也就是学习速率,控制更新的幅度。
    alpha = 0.001
    #最大迭代次数
    maxCycles = 500
    # weights 代表回归系数, 此处的 ones((n,1)) 创建一个长度和特征数相同的矩阵,其中的数全部都是 1
    weights = np.ones((n,1))
    for k in range(maxCycles):
        # m*3 的矩阵 * 3*1 的矩阵 = m*1的矩阵
        h = sigmoid(dataMatrix*weights)
        # labelMat是实际值,labelMat - 是两个向量相减
        error = (labelMat - h) 
        # 0.001* (3*m)*(m*1) 表示在每一个列上的一个误差情况,最后得出 x1,x2,xn的系数的偏移量
         # 矩阵乘法,最后得到回归系数
        weights = weights + alpha * dataMatrix.transpose() * error
    return np.array(weights)

测试及其代码如下:

dataMat, labelMat = loadDataSet("TestSet.txt")
dataArr = np.array(dataMat)
weights = gradAscent(dataArr, labelMat)
plotBestFit(dataArr, labelMat, weights)

改善后1

# Logistic 回归梯度上升优化算法
def gradAscent0(dataMatIn, classLabels):
    #返回dataMatrix的大小。m为行数,n为列数。
    m,n = np.shape(dataMatIn)
    #移动步长,也就是学习速率,控制更新的幅度。
    alpha = 0.01
    #函数ones创建一个n*1的数组
    weights = np.ones(n)
    for i in range(m):
        #sum(dataMatrix[i]*weights)为了求 f(x)的值, f(x)=a1*x1+a2*x2+..+an*xn,此处求出的 h 是一个具体的数值,而不是一个矩阵
        h = sigmoid(np.sum(dataMatIn[i]*weights))
        # 计算真实类别与预测类别之间的差值,然后按照该差值调整回归系数
        error = (classLabels[i] - h) 
        weights = weights + alpha * error * dataMatIn[i]
    return weights

测试及其代码如下:

dataMat, labelMat = loadDataSet("TestSet.txt")
dataArr = np.array(dataMat)
weights = gradAscent0(dataArr, labelMat)
plotBestFit(dataArr, labelMat, weights)

可以看到,随机梯度上升算法与梯度上升算法在代码上很相似,但也有一些区别: 第一,后者的变量 h 和误差 error 都是向量,而前者则全是数值;第二,前者没有矩阵的转换过程,所有变量的数据类型都是 NumPy 数组。

判断优化算法优劣的可靠方法是看它是否收敛,也就是说参数是否达到了稳定值,是否还会不断地变化?下图展示了随机梯度上升算法在 200 次迭代过程中回归系数的变化情况。其中的系数2,也就是 X2 只经过了 50 次迭代就达到了稳定值,但系数 1 和 0 则需要更多次的迭代。如下图所示:

针对这个问题,我们改进了之前的随机梯度上升算法,如下:
改善后2

# 随机梯度上升算法(随机化)
def stocGradAscent1(dataMatIn, classLabels, numIter=150):
    m,n = np.shape(dataMatIn)
     # 创建与列数相同的矩阵的系数矩阵,1行3列
    weights = np.ones(n)  
    # 随机梯度, 循环150,观察是否收敛
    for j in range(numIter):
        # [0, 1, 2 .. m-1]
        dataIndex = list(range(m))
        for i in range(m):
            # i和j的不断增大,导致alpha的值不断减少,但是不为0
            alpha = 4/(1.0+j+i)+0.01    # alpha 会随着迭代不断减小,但永远不会减小到0,因为后边还有一个常数项0.0001
            # 随机产生一个 0~len()之间的一个值
            # random.uniform(x, y) 方法将随机生成下一个实数,它在[x,y]范围内,x是这个范围内的最小值,y是这个范围内的最大值。
            randIndex = int(random.uniform(0,len(dataIndex)))
            # sum(dataMatrix[i]*weights)为了求 f(x)的值, f(x)=a1*x1+b2*x2+..+nn*xn
            h = sigmoid(sum(dataMatIn[dataIndex[randIndex]]*weights))
            error = classLabels[dataIndex[randIndex]] - h
            weights = weights + alpha * error * dataMatIn[dataIndex[randIndex]]
            del(dataIndex[randIndex])
    return weights

测试及其代码如下:

dataMat, labelMat = loadDataSet("TestSet.txt")
dataArr = np.array(dataMat)
weights = gradAscent1(dataArr, labelMat)
plotBestFit(dataArr, labelMat, weights)

上面的改进版随机梯度上升算法,我们修改了两处代码。

第一处改进为 alpha 的值。alpha 在每次迭代的时候都会调整,这回缓解上面波动图的数据波动或者高频波动。另外,虽然 alpha 会随着迭代次数不断减少,但永远不会减小到 0,因为我们在计算公式中添加了一个常数项。

第二处修改为 randIndex 更新,这里通过随机选取样本拉来更新回归系数。这种方法将减少周期性的波动。这种方法每次随机从列表中选出一个值,然后从列表中删掉该值(再进行下一次迭代)。

程序运行之后能看到类似于上图的结果图。

改善后完整的代码如下:

import matplotlib.pyplot as plt
import numpy as np
import random
#解析数据
def loadDataSet(file_name):
    '''
    Desc: 
        加载并解析数据
    Args:
        file_name -- 要解析的文件路径
    Returns:
        dataMat -- 原始数据的特征
        labelMat -- 原始数据的标签,也就是每条样本对应的类别。即目标向量
    '''
    # dataMat为原始数据, labelMat为原始数据的标签
    dataMat = []
    labelMat = []
    fr = open(file_name)
    for line in fr.readlines():
        lineArr = line.strip().split()
        # 为了方便计算,我们将 X0 的值设为 1.0 ,也就是在每一行的开头添加一个 1.0 作为 X0
        dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])
        labelMat.append(int(lineArr[2]))
    return dataMat, labelMat

# sigmoid阶跃函数
def sigmoid(inX):
    # return 1.0 / (1 + exp(-inX))
    return 1.0/(1+np.exp(-inX)) 

# Logistic 回归梯度上升优化算法
def gradAscent(dataMatIn, classLabels):
    #用mat函数将列表转换成numpy矩阵
    dataMatrix = np.mat(dataMatIn)
    #用mat函数将列表转换成numpy矩阵,并进行转置
    labelMat = np.mat(classLabels).transpose()
    #返回dataMatrix的大小。m为行数,n为列数。
    m,n = np.shape(dataMatrix)
    #移动步长,也就是学习速率,控制更新的幅度。
    alpha = 0.001
    #最大迭代次数
    maxCycles = 500
    # weights 代表回归系数, 此处的 ones((n,1)) 创建一个长度和特征数相同的矩阵,其中的数全部都是 1
    weights = np.ones((n,1))
    for k in range(maxCycles):
        # m*3 的矩阵 * 3*1 的矩阵 = m*1的矩阵
        h = sigmoid(dataMatrix*weights)
        # labelMat是实际值,labelMat - 是两个向量相减
        error = (labelMat - h) 
        # 0.001* (3*m)*(m*1) 表示在每一个列上的一个误差情况,最后得出 x1,x2,xn的系数的偏移量
         # 矩阵乘法,最后得到回归系数
        weights = weights + alpha * dataMatrix.transpose() * error
    return np.array(weights)

# Logistic 回归梯度上升优化算法
def gradAscent0(dataMatIn, classLabels):
    #返回dataMatrix的大小。m为行数,n为列数。
    m,n = np.shape(dataMatIn)
    #移动步长,也就是学习速率,控制更新的幅度。
    alpha = 0.01
    #函数ones创建一个n*1的数组
    weights = np.ones(n)
    for i in range(m):
        #sum(dataMatrix[i]*weights)为了求 f(x)的值, f(x)=a1*x1+a2*x2+..+an*xn,此处求出的 h 是一个具体的数值,而不是一个矩阵
        h = sigmoid(np.sum(dataMatIn[i]*weights))
        # 计算真实类别与预测类别之间的差值,然后按照该差值调整回归系数
        error = (classLabels[i] - h) 
        weights = weights + alpha * error * dataMatIn[i]
    return weights

# 画出数据集和 Logistic 回归最佳拟合直线的函数
def plotBestFit(dataArr, labelMat, weights):
    n = np.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 = np.arange(-3.0, 3.0, 0.1)
    y = (-weights[0]-weights[1]*x)/weights[2]
    ax.plot(x, y)
    plt.xlabel('X'); plt.ylabel('Y')
    plt.show()

#测试算法: 使用 Logistic 回归进行分类
def testLR():
    dataMat, labelMat = loadDataSet("TestSet.txt")
    dataArr = np.array(dataMat)
    weights = stocGradAscent1(dataArr, labelMat)
    plotBestFit(dataArr, labelMat, weights)
    

# 随机梯度上升算法(随机化)
def stocGradAscent1(dataMatIn, classLabels, numIter=150):
    m,n = np.shape(dataMatIn)
     # 创建与列数相同的矩阵的系数矩阵,1行3列
    weights = np.ones(n)  
    # 随机梯度, 循环150,观察是否收敛
    for j in range(numIter):
        # [0, 1, 2 .. m-1]
        dataIndex = list(range(m))
        for i in range(m):
            # i和j的不断增大,导致alpha的值不断减少,但是不为0
            alpha = 4/(1.0+j+i)+0.01    # alpha 会随着迭代不断减小,但永远不会减小到0,因为后边还有一个常数项0.0001
            # 随机产生一个 0~len()之间的一个值
            # random.uniform(x, y) 方法将随机生成下一个实数,它在[x,y]范围内,x是这个范围内的最小值,y是这个范围内的最大值。
            randIndex = int(random.uniform(0,len(dataIndex)))
            # sum(dataMatrix[i]*weights)为了求 f(x)的值, f(x)=a1*x1+b2*x2+..+nn*xn
            h = sigmoid(sum(dataMatIn[dataIndex[randIndex]]*weights))
            error = classLabels[dataIndex[randIndex]] - h
            weights = weights + alpha * error * dataMatIn[dataIndex[randIndex]]
            del(dataIndex[randIndex])
    return weights

testLR()
最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 151,511评论 1 330
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 64,495评论 1 273
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 101,595评论 0 225
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 42,558评论 0 190
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 50,715评论 3 270
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 39,672评论 1 192
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 31,112评论 2 291
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 29,837评论 0 181
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 33,417评论 0 228
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 29,928评论 2 232
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 31,316评论 1 242
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 27,773评论 2 234
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 32,253评论 3 220
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 25,827评论 0 8
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 26,440评论 0 180
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 34,523评论 2 249
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 34,583评论 2 249

推荐阅读更多精彩内容