Logistic回归

我的微信公众号： s406205391; 欢迎大家一起学习，一起进步！！！

假设现在有一些数据点，我们用一条直线对这些点进行拟合（该线称为最佳拟合直线），这个拟合过程就称作回归。利用Logistic回归进行分类的主要思想是：根据现有数据对分类界线建立回归公式，以此进行分类。训练分类器时的做法就是寻找最佳拟合参数，使用的是最优化算法。

Logistic回归的一般过程：

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

基于Logistic回归和Sigmoid函数的分类

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

Logistic回归中，使用了Sigmoid阶跃函数，该函数的特点是：该函数在跳跃点上可以从0瞬间跳跃到1。

image

基于这个函数，我们可以将回归方程，带入这个函数。进而得到一个范围在0~1之间的数字。任何大于0.5的数据被分为1类，小于0.5即被归入0类。这样便可以很好的完成分类。公式如下：

Sigmoid： $h(x)=\frac{1}{1+e^{-z}}$ <font size=2, color='grey'>z表示回归方程</font>

回归方程：z = $w_0x_0 + w_1x_1 + w_2x_2 + ... + w_nx_n$ , 也可以写成 $z=w^Tx$

接下来，我们要做的事情，便是找到w的最佳值，一般有梯度下降法、梯度上升法、随机梯度下降法、随机梯度上升法等。其最终推导公式如下：

$\theta:=\theta + \alpha\cdot[\frac{1}{m}]\cdot x^T\cdot(g(x\cdot\theta)-y)$ <font size=2, color='grey'>g(x)即为Sigmoid函数，α为随机常数，</font>

具体推到过程，可以参考这篇文章。

梯度上升/下降法

梯度上升和梯度下降其实差不多，只需要把上面的公式的加号变成减号就行了。

本文所用的数据集(提取码：ogge)。下面是一个梯度上升算法的具体实现。load_data_set()函数，从文本中提取了一组数据，作为我们的训练的数据集。sigmoid()函数则是Sigmoid的计算。grad_ascent()函数则是梯度上升求参数最佳值的方法。其实主要就是运用了上面的推导公式。

import matplotlib.pyplot as plt
import numpy as np


def load_data_set():
    """读取测试数据集"""
    data_mat = []
    label_mat = []
    for line in open('data\\Ch05\\testSet.txt'):
        line = line.strip().split()
        data_mat.append([1.0, float(line[0]), float(line[1])])
        label_mat.append(int(line[2]))
    return data_mat, label_mat


def sigmoid(in_x):
    """sigmoid阶跃函数"""

    return 1 / (1 + np.exp(-in_x))


def grad_ascent(data_mat_in, class_labels):
    """使用梯度上升法，求得系数"""

    data_matrix = np.mat(data_mat_in)
    label_mat = np.mat(class_labels).transpose()  # 转置矩阵
    m, n = np.shape(data_matrix)
    alpha = 0.001
    max_cycles = 500
    weights = np.ones((n, 1))
    for k in range(max_cycles):
        h = sigmoid(data_matrix * weights)
        error = (label_mat - h)
        weights = weights + alpha * data_matrix.transpose() * error
    return weights


def plot_best_fit(weights):
    data_mat, label_mat = load_data_set()
    data_array = np.array(data_mat)
    n = np.shape(data_array)[0]
    xcord1 = []
    ycord1 = []
    xcord2 = []
    ycord2 = []
    for i in range(n):
        if int(label_mat[i]) == 1:
            xcord1.append(data_array[i, 1])
            ycord1.append(data_array[i, 2])
        else:
            xcord2.append(data_array[i, 1])
            ycord2.append(data_array[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('x1')
        plt.ylabel('x2')
        plt.show()

if __name__ == '__main__':
    d, l = load_data_set()
    w = grad_ascent(d, l)
    plot_best_fit(w)

下图是梯度上升法，迭代500次的分类结果。从结果上看，分类的还是比较好的。

image

随机梯度上升法

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

def sto_grad_ascent(data_matrix, class_labels):
    """随机梯度上升法"""

    m, n = np.shape(data_matrix)
    alpha = 0.001
    weights = np.ones(n)
    for i in range(m):
        h = sigmoid(data_matrix[i] * weights)
        error = (class_labels[i] - h)
        weights = weights + alpha * data_matrix[i] * error
    return weights

if __name__ == '__main__':
    d, l = load_data_set()
    w = sto_grad_ascent(np.array(d), l)
    plot_best_fit(w)

下图是随机梯度上升法，拟合的结果。不如梯度上升法完美，但是也需要考虑，后者是迭代了500次才得到的。一个判断优化算法优劣的可靠方法是看它是否收敛，也就是说参数是否达到了稳定值，是否还会不断地变化。

image

改进的随机梯度上升算法

对随机梯度算法进行一些改进，得到下面的代码。其改进之处在于。一方面，alpha在每次迭代的时候都会调整，这会缓解随机梯度上升算法的数据波动。另外，虽然alpha会随着迭代次数不断减小，但永远不会减小到0，这是因为公式中存在常数项，必须保证在多次迭代后，该常数项仍然具有一定的影响。另一个改进之处在于，这里通过随机选取样本来更新回归系数。这种方法将减少周期性波动。

def stoc_grad_ascent(data_matrix, class_labels, num_iter=500):
    """改进的随机梯度上升算法"""

    m, n = np.shape(data_matrix)
    weights = np.ones(n)
    for j in range(num_iter):
        data_index = list(range(m))
        for i in range(m):
            alpha = 4 / (1 + j + i) + 0.01
            rand_index = int(np.random.uniform(0, len(data_index)))
            h = sigmoid(sum(data_matrix[rand_index] * weights))
            error = class_labels[rand_index] - h
            weights = weights + alpha * error * data_matrix[rand_index]
            del data_index[rand_index]
    return weights

改进之后的拟合结果，要明显优于改进之前得结果。

在这里插入图片描述

示例：从疝气病症预测病马病死率

接下来将使用Logistic回归来预测患有疝病的马的存活问题。这里的数据包含368个样本和28个特征。我们在日常数据处理中，常会遇到数据集包含缺失值的情况，那么在Logistic回归算法中，缺失值是否可以用0来代替呢？

答案是可以的。我们观察回归系数的更新公式：

weights = weights + alpha * error * dataMatrix[randIndex]

如果dataMatrix的某个特征对应值为0，那么该特征的系数将不做更新，即：

weights = weights

另外，由于sigmoid(0)=0.5, 即它对结果的预测不具有任何倾向性。基于上述原因，将缺失值用0代替既可以保留现有数据，也不需要对优化的算法进行修改。

下面的代码中，classify_vector()构建了一个Logistic分类器，只需要将训练你的到的回归系数带入进去，就可以对样本进行分类了。colic_test()函数则是示例的具体代码

def classify_vector(in_x, weithts):
    prob = sigmoid(sum(in_x * weithts))
    if prob > 0.5:
        return 1.0
    else:
        return 0.0


def colic_test():
    fr_train = open('data\\Ch05\\horseColicTraining.txt')
    fr_test = open('data\\Ch05\\horseColicTest.txt')

    training_set = []
    training_labels = []
    for line in fr_train.readlines():
        line = line.strip().split('\t')
        line_arr = list(map(float, line[0: 21]))
        training_set.append(line_arr)
        training_labels.append(float(line[21]))

    train_weights = stoc_grad_ascent(np.array(training_set), training_labels, 500)
    error_count = 0
    num_test_vec = 0
    for line in fr_test.readlines():
        num_test_vec += 1
        line = line.strip().split('\t')
        line_arr = list(map(float, line[0: 21]))
        if int(classify_vector(np.array(line_arr), train_weights)) != int(line[21]):
            error_count += 1

    error_rate = error_count / num_test_vec
    print(f"the error rate of this test is: {error_rate}")
    return error_rate

if __name__ == '__main__':
    colic_test()