05机器学习实战之Logistic 回归

Logistic 回归 概述

Logistic 回归 或者叫逻辑回归 虽然名字有回归，可是它是用来作分类的。其主要思想是: 根据现有数据对分类边界线(Decision Boundary)创建回归公式，以此进行分类。

须知概念

Sigmoid 函数

回归 概念

假设如今有一些数据点，咱们用一条直线对这些点进行拟合（这条直线称为最佳拟合直线），这个拟合的过程就叫作回归。进而能够获得对这些点的拟合直线方程，那么咱们根据这个回归方程，怎么进行分类呢？请看下面。

二值型输出分类函数

咱们想要的函数应该是: 能接受全部的输入而后预测出类别。例如，在两个类的状况下，上述函数输出 0 或 1.或许你以前接触过具备这种性质的函数，该函数称为 海维塞得阶跃函数(Heaviside step function)，或者直接称为 单位阶跃函数。然而，海维塞得阶跃函数的问题在于: 该函数在跳跃点上从 0 瞬间跳跃到 1，这个瞬间跳跃过程有时很难处理。幸亏，另外一个函数也有相似的性质（能够输出 0 或者 1 的性质），且数学上更易处理，这就是 Sigmoid 函数。 Sigmoid 函数具体的计算公式以下:

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

所以，为了实现 Logistic 回归分类器，咱们能够在每一个特征上都乘以一个回归系数（以下公式所示），而后把全部结果值相加，将这个总和代入 Sigmoid 函数中，进而获得一个范围在 0~1 之间的数值。任何大于 0.5 的数据被分入 1 类，小于 0.5 即被纳入 0 类。因此，Logistic 回归也是一种几率估计，好比这里Sigmoid 函数得出的值为0.5，能够理解为给定数据和参数，数据被分入 1 类的几率为0.5。想对Sigmoid 函数有更多了解，能够点开此连接跟此函数互动。

基于最优化方法的回归系数肯定

Sigmoid 函数的输入记为 z ，由下面公式获得:

若是采用向量的写法，上述公式能够写成  ，它表示将这两个数值向量对应元素相乘而后所有加起来即获得 z 值。其中的向量 x 是分类器的输入数据，向量 w 也就是咱们要找到的最佳参数（系数），从而使得分类器尽量地精确。为了寻找该最佳参数，须要用到最优化理论的一些知识。咱们这里使用的是——梯度上升法（Gradient Ascent）。

梯度上升法

梯度上升法的思想

要找到某函数的最大值，最好的方法是沿着该函数的梯度方向探寻。若是梯度记为 ▽ ，则函数 f(x, y) 的梯度由下式表示:

这个梯度意味着要沿 x 的方向移动  ，沿 y 的方向移动  。其中，函数f(x, y) 必需要在待计算的点上有定义而且可微。下图是一个具体的例子。

上图展现的，梯度上升算法到达每一个点后都会从新估计移动的方向。从 P0 开始，计算完该点的梯度，函数就根据梯度移动到下一点 P1。在 P1 点，梯度再次被从新计算，并沿着新的梯度方向移动到 P2 。如此循环迭代，直到知足中止条件。迭代过程当中，梯度算子老是保证咱们能选取到最佳的移动方向。

上图中的梯度上升算法沿梯度方向移动了一步。能够看到，梯度算子老是指向函数值增加最快的方向。这里所说的是移动方向，而未提到移动量的大小。该量值称为步长，记做 α 。用向量来表示的话，梯度上升算法的迭代公式以下:

该公式将一直被迭代执行，直至达到某个中止条件为止，好比迭代次数达到某个指定值或者算法达到某个能够容许的偏差范围。

介绍一下几个相关的概念：

例如：y = w0 + w1x1 + w2x2 + ... + wnxn
梯度：参考上图的例子，二维图像，x方向表明第一个系数，也就是 w1，y方向表明第二个系数也就是 w2，这样的向量就是梯度。
α：上面的梯度算法的迭代公式中的阿尔法，这个表明的是移动步长（step length）。移动步长会影响最终结果的拟合程度，最好的方法就是随着迭代次数更改移动步长。
步长通俗的理解，100米，若是我一步走10米，我须要走10步；若是一步走20米，我只须要走5步。这里的一步走多少米就是步长的意思。
▽f(w)：表明沿着梯度变化的方向。

问：有人会好奇为何有些书籍上说的是梯度降低法（Gradient Decent）?

答： 其实这个两个方法在此状况下本质上是相同的。关键在于代价函数（cost function）或者叫目标函数（objective function）。若是目标函数是损失函数，那就是最小化损失函数来求函数的最小值，就用梯度降低。 若是目标函数是似然函数（Likelihood function），就是要最大化似然函数来求函数的最大值，那就用梯度上升。在逻辑回归中， 损失函数和似然函数无非就是互为正负关系。

只须要在迭代公式中的加法变成减法。所以，对应的公式能够写成

局部最优现象 （Local Optima）

上图表示参数 θ 与偏差函数 J(θ) 的关系图 (这里的偏差函数是损失函数，因此咱们要最小化损失函数)，红色的部分是表示 J(θ) 有着比较高的取值，咱们须要的是，可以让 J(θ) 的值尽可能的低。也就是深蓝色的部分。θ0，θ1 表示 θ 向量的两个维度（此处的θ0，θ1是x0和x1的系数，也对应的是上文w0和w1）。

可能梯度降低的最终点并不是是全局最小点，多是一个局部最小点，如咱们上图中的右边的梯度降低曲线，描述的是最终到达一个局部最小点，这是咱们从新选择了一个初始点获得的。

看来咱们这个算法将会在很大的程度上被初始点的选择影响而陷入局部最小点。

Logistic 回归 原理

Logistic 回归 工做原理

每一个回归系数初始化为 1
重复 R 次:
    计算整个数据集的梯度
    使用 步长 x 梯度 更新回归系数的向量
返回回归系数

Logistic 回归 开发流程

收集数据: 采用任意方法收集数据
准备数据: 因为须要进行距离计算，所以要求数据类型为数值型。另外，结构化数据格式则最佳。
分析数据: 采用任意方法对数据进行分析。
训练算法: 大部分时间将用于训练，训练的目的是为了找到最佳的分类回归系数。
测试算法: 一旦训练步骤完成，分类将会很快。
使用算法: 首先，咱们须要输入一些数据，并将其转换成对应的结构化数值；接着，基于训练好的回归系数就能够对这些数值进行简单的回归计算，断定它们属于哪一个类别；在这以后，咱们就能够在输出的类别上作一些其余分析工做。

Logistic 回归 算法特色

优势: 计算代价不高，易于理解和实现。
缺点: 容易欠拟合，分类精度可能不高。
适用数据类型: 数值型和标称型数据。

附加 方向导数与梯度

Logistic 回归 项目案例

项目概述

在一个简单的数据集上，采用梯度上升法找到 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

绘制在图中，以下图所示:

准备数据: 因为须要进行距离计算，所以要求数据类型为数值型。另外，结构化数据格式则最佳

import numpy as np

data_arr = []
label_arr = []
f = open('D:\\mlInAction\\data\\5.Logistic\\TestSet.txt', 'r')
for line in f.readlines():
    line_arr = line.strip().split()
    # 为了方便计算，咱们将 X0 的值设为 1.0 ，也就是在每一行的开头添加一个 1.0 做为 X0
    data_arr.append([1.0, np.float(line_arr[0]), np.float(line_arr[1])])
    label_arr.append(int(line_arr[2]))

分析数据: 采用任意方法对数据进行分析，此处不须要

训练算法: 使用梯度上升找到最佳参数

定义sigmoid阶跃函数

def sigmoid(x):
    # 这里其实很是有必要解释一下，会出现的错误 RuntimeWarning: overflow encountered in exp
    # 这个错误在学习阶段虽然能够忽略，可是咱们至少应该知道为何
    # 这里是由于咱们输入的有的 x 实在是过小了，好比 -6000之类的，那么计算一个数字 np.exp(6000)这个结果太大了，无法表示，因此就溢出了
    # 若是是计算 np.exp（-6000），这样虽然也会溢出，可是这是下溢，就是表示成零
    # 去网上搜了不少方法，好比 使用bigfloat这个库（我居然没有安装成功，就不尝试了，反正应该是有用的
    return 1.0 / (1 + np.exp(-x))

Logistic 回归梯度上升优化算法

def grad_ascent(data_arr, class_labels):
    """
    梯度上升法，其实就是由于使用了极大似然估计，这个你们有必要去看推导，只看代码感受不太够
    :param data_arr: 传入的就是一个普通的数组，固然你传入一个二维的ndarray也行
    :param class_labels: class_labels 是类别标签，它是一个 1*100 的行向量。
                    为了便于矩阵计算，须要将该行向量转换为列向量，作法是将原向量转置，再将它赋值给label_mat
    :return: 
    """
    # 注意一下，我把原来 data_mat_in 改为data_arr,由于传进来的是一个数组，用这个比较不容易搞混
    # turn the data_arr to numpy matrix
    data_mat = np.mat(data_arr)
    # 变成矩阵以后进行转置
    label_mat = np.mat(class_labels).transpose()
    # m->数据量，样本数 n->特征数
    m, n = np.shape(data_mat)
    # 学习率，learning rate
    alpha = 0.001
    # 最大迭代次数，伪装迭代这么屡次就能收敛2333
    max_cycles = 500
    # 生成一个长度和特征数相同的矩阵，此处n为3 -> [[1],[1],[1]]
    # weights 表明回归系数， 此处的 ones((n,1)) 建立一个长度和特征数相同的矩阵，其中的数所有都是 1
    weights = np.ones((n, 1))
    for k in range(max_cycles):
        # 这里是点乘  m x 3 dot 3 x 1
        h = sigmoid(data_mat * weights)
        error = label_mat - h
        # 这里比较建议看一下推导，为何这么作能够，这里已是求导以后的
        weights = weights + alpha * data_mat.transpose() * error
    return weights

你们看到这儿可能会有一些疑惑，就是，咱们在迭代中更新咱们的回归系数，后边的部分是怎么计算出来的？为何会是 alpha * dataMatrix.transpose() * error ?由于这就是咱们所求的梯度，也就是对 f(w) 对 w 求一阶导数。具体推导以下:

 

可参考http://blog.csdn.net/achuo/article/details/51160101

画出数据集和 Logistic 回归最佳拟合直线的函数

import matplotlib.pyplot as plt

def plot_best_fit(data_mat, label_mat, weights):
    """
    可视化
    :param weights: 
    :return: 
    """
    data_arr = np.array(data_mat)
    n = np.shape(data_mat)[0]
    x_cord1 = []
    y_cord1 = []
    x_cord2 = []
    y_cord2 = []
    for i in range(n):
        if int(label_mat[i]) == 1:
            x_cord1.append(data_arr[i, 1])
            y_cord1.append(data_arr[i, 2])
        else:
            x_cord2.append(data_arr[i, 1])
            y_cord2.append(data_arr[i, 2])
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.scatter(x_cord1, y_cord1, s=30, color='k', marker='^')
    ax.scatter(x_cord2, y_cord2, s=30, color='red', marker='s')
    x = np.arange(-3.0, 3.0, 0.1)
    # print(x)
    y = (-weights[0] - weights[1] * x) / weights[2]
    # type(y)
    y = np.ravel(y)  # y原来是一个二维，须要转化为1维
    """
    y的由来，卧槽，是否是没看懂？
    首先理论上是这个样子的。
    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   
    """
    ax.plot(x, y)
    plt.xlabel('x1')
    plt.ylabel('y1')
    plt.show()

注意

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

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

如下是随机梯度上升算法的实现代码:全部回归系数初始化为 1
对数据集中每一个样本
    计算该样本的梯度
    使用 alpha * gradient 更新回归系数值
返回回归系数值

def stoc_grad_ascent0(data_mat, class_labels):
    """
    随机梯度上升，只使用一个样本点来更新回归系数
    :param data_mat: 输入数据的数据特征（除去最后一列）,ndarray
    :param class_labels: 输入数据的类别标签（最后一列数据）
    :return: 获得的最佳回归系数
    """
    m, n = np.shape(data_mat)
    alpha = 0.01
    weights = np.ones(n)
    for i in range(m):
        # sum(data_mat[i]*weights)为了求 f(x)的值， f(x)=a1*x1+b2*x2+..+nn*xn,
        # 此处求出的 h 是一个具体的数值，而不是一个矩阵
        h = sigmoid(sum(data_mat[i] * weights))
        error = class_labels[i] - h
        # 仍是和上面同样，这个先去看推导，再写程序
        weights = weights + alpha * error * data_mat[i]
    return weights

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

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

针对这个问题，咱们改进了以前的随机梯度上升算法，以下:

def stoc_grad_ascent1(data_mat, class_labels, num_iter=150):
    """
    改进版的随机梯度上升，使用随机的一个样原本更新回归系数
    :param data_mat: 输入数据的数据特征（除去最后一列）,ndarray
    :param class_labels: 输入数据的类别标签（最后一列数据
    :param num_iter: 迭代次数
    :return: 获得的最佳回归系数
    """
    m, n = np.shape(data_mat)
    weights = np.ones(n)
    for j in range(num_iter):
        # 这里必需要用list，否则后面的del无法使用
        data_index = list(range(m))
        for i in range(m):
            # i和j的不断增大，致使alpha的值不断减小，可是不为0
            alpha = 4 / (1.0 + j + i) + 0.01
            # 随机产生一个 0～len()之间的一个值
            # random.uniform(x, y) 方法将随机生成下一个实数，它在[x,y]范围内,x是这个范围内的最小值，y是这个范围内的最大值。
            rand_index = int(np.random.uniform(0, len(data_index)))
            h = sigmoid(np.sum(data_mat[data_index[rand_index]] * weights))
            error = class_labels[data_index[rand_index]] - h
            weights = weights + alpha * error * data_mat[data_index[rand_index]]
            del(data_index[rand_index])
    return weights

In [54]:

import numpy as np

data_arr = []
label_arr = []
f = open('D:\\mlInAction\\data\\5.Logistic\\TestSet.txt', 'r')
for line in f.readlines():
    line_arr = line.strip().split()
    # 为了方便计算，咱们将 X0 的值设为 1.0 ，也就是在每一行的开头添加一个 1.0 做为 X0
    data_arr.append([1.0, np.float(line_arr[0]), np.float(line_arr[1])])
    label_arr.append(int(line_arr[2]))

In [55]:

data_arr

Out[55]:

[[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]]

In [56]:

label_arr

Out[56]:

[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]

In [57]:

def sigmoid(x):
    # 这里其实很是有必要解释一下，会出现的错误 RuntimeWarning: overflow encountered in exp
    # 这个错误在学习阶段虽然能够忽略，可是咱们至少应该知道为何
    # 这里是由于咱们输入的有的 x 实在是过小了，好比 -6000之类的，那么计算一个数字 np.exp(6000)这个结果太大了，无法表示，因此就溢出了
    # 若是是计算 np.exp（-6000），这样虽然也会溢出，可是这是下溢，就是表示成零
    # 去网上搜了不少方法，好比 使用bigfloat这个库（我居然没有安装成功，就不尝试了，反正应该是有用的
    return 1.0 / (1 + np.exp(-x))

In [58]:

def grad_ascent(data_arr, class_labels):
    """
    梯度上升法，其实就是由于使用了极大似然估计，这个你们有必要去看推导，只看代码感受不太够
    :param data_arr: 传入的就是一个普通的数组，固然你传入一个二维的ndarray也行
    :param class_labels: class_labels 是类别标签，它是一个 1*100 的行向量。
                    为了便于矩阵计算，须要将该行向量转换为列向量，作法是将原向量转置，再将它赋值给label_mat
    :return: 
    """
    # 注意一下，我把原来 data_mat_in 改为data_arr,由于传进来的是一个数组，用这个比较不容易搞混
    # turn the data_arr to numpy matrix
    data_mat = np.mat(data_arr)
    # 变成矩阵以后进行转置
    label_mat = np.mat(class_labels).transpose()
    # m->数据量，样本数 n->特征数
    m, n = np.shape(data_mat)
    # 学习率，learning rate
    alpha = 0.001
    # 最大迭代次数，伪装迭代这么屡次就能收敛2333
    max_cycles = 500
    # 生成一个长度和特征数相同的矩阵，此处n为3 -> [[1],[1],[1]]
    # weights 表明回归系数， 此处的 ones((n,1)) 建立一个长度和特征数相同的矩阵，其中的数所有都是 1
    weights = np.ones((n, 1))
    for k in range(max_cycles):
        # 这里是点乘  m x 3 dot 3 x 1
        h = sigmoid(data_mat * weights)
        error = label_mat - h
        # 这里比较建议看一下推导，为何这么作能够，这里已是求导以后的
        weights = weights + alpha * data_mat.transpose() * error
    return weights

In [69]:

weights = grad_ascent(data_arr, label_arr)
weights

Out[69]:

matrix([[ 4.12414349],
        [ 0.48007329],
        [-0.6168482 ]])

In [76]:

import matplotlib.pyplot as plt


def plot_best_fit(data_mat, label_mat, weights):
    """
    可视化
    :param weights: 
    :return: 
    """
    data_arr = np.array(data_mat)
    n = np.shape(data_mat)[0]
    x_cord1 = []
    y_cord1 = []
    x_cord2 = []
    y_cord2 = []
    for i in range(n):
        if int(label_mat[i]) == 1:
            x_cord1.append(data_arr[i, 1])
            y_cord1.append(data_arr[i, 2])
        else:
            x_cord2.append(data_arr[i, 1])
            y_cord2.append(data_arr[i, 2])
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.scatter(x_cord1, y_cord1, s=30, color='k', marker='^')
    ax.scatter(x_cord2, y_cord2, s=30, color='red', marker='s')
    x = np.arange(-3.0, 3.0, 0.1)
    # print(x)
    y = (-weights[0] - weights[1] * x) / weights[2]
    # type(y)
    y = np.ravel(y)  # y原来是一个二维，须要转化为1维
    """
    y的由来，卧槽，是否是没看懂？
    首先理论上是这个样子的。
    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   
    """
    ax.plot(x, y)
    plt.xlabel('x1')
    plt.ylabel('y1')
    plt.show()

In [77]:

plot_best_fit(data_arr, label_arr, weights)

 

In [72]:

def stoc_grad_ascent0(data_mat, class_labels):
    """
    随机梯度上升，只使用一个样本点来更新回归系数
    :param data_mat: 输入数据的数据特征（除去最后一列）,ndarray
    :param class_labels: 输入数据的类别标签（最后一列数据）
    :return: 获得的最佳回归系数
    """
    m, n = np.shape(data_mat)
    alpha = 0.01
    weights = np.ones(n)
    for i in range(m):
        # sum(data_mat[i]*weights)为了求 f(x)的值， f(x)=a1*x1+b2*x2+..+nn*xn,
        # 此处求出的 h 是一个具体的数值，而不是一个矩阵
        h = sigmoid(sum(data_mat[i] * weights))
        error = class_labels[i] - h
        # 仍是和上面同样，这个先去看推导，再写程序
        weights = weights + alpha * error * data_mat[i]
    return weights

In [78]:

def stoc_grad_ascent1(data_mat, class_labels, num_iter=150):
    """
    改进版的随机梯度上升，使用随机的一个样原本更新回归系数
    :param data_mat: 输入数据的数据特征（除去最后一列）,ndarray
    :param class_labels: 输入数据的类别标签（最后一列数据
    :param num_iter: 迭代次数
    :return: 获得的最佳回归系数
    """
    m, n = np.shape(data_mat)
    weights = np.ones(n)
    for j in range(num_iter):
        # 这里必需要用list，否则后面的del无法使用
        data_index = list(range(m))
        for i in range(m):
            # i和j的不断增大，致使alpha的值不断减小，可是不为0
            alpha = 4 / (1.0 + j + i) + 0.01
            # 随机产生一个 0～len()之间的一个值
            # random.uniform(x, y) 方法将随机生成下一个实数，它在[x,y]范围内,x是这个范围内的最小值，y是这个范围内的最大值。
            rand_index = int(np.random.uniform(0, len(data_index)))
            h = sigmoid(np.sum(data_mat[data_index[rand_index]] * weights))
            error = class_labels[data_index[rand_index]] - h
            weights = weights + alpha * error * data_mat[data_index[rand_index]]
            del(data_index[rand_index])
    return weights

In [81]:

weights1 = stoc_grad_ascent1(np.array(data_arr), np.array(label_arr))
weights1

Out[81]:

array([13.97204276,  1.24615972, -1.9102347 ])

In [82]:

plot_best_fit(data_arr, label_arr, weights1)