Logistic回归

项目代码

Logistic回归

假设有一堆数据点，用一条直线对这些点进行拟合（该条直线成为最佳拟合直线），这个拟合过程就称为回归。
利用Logistic回归进行分类的思想：根据现有数据对分类边界线建立回归公式，以此进行分类。这里的回归源于最佳拟合，表示要找到最佳拟合参数集。

利用Logistic回归如何分类

假设在二分类情况下，我们希望有某个函数在接受所有输入之后可以预测出类别（即可以输出0或1，表示不同类别）。这种函数称为海维赛德阶跃函数（或者称为单位阶跃函数）。然而海维赛德阶跃函数在跳跃点上从0瞬间跳跃到1，这个瞬间跳跃过程有时很难处理。因此需要一个更好的函数，这就是这节用到的sigmoid函数，其数学表达式为：
$sigmoid(z)=\frac{1}{1+e^{-z}}$
下图是sigmoid函数在不同坐标尺度下的曲线图：

当$z=0$时，$sigmoid(z)=\frac{1}{2}$；随着$z$的增大，$sigmoid(z)$函数值逼近1；随着$z$的减小，$sigmoid(z)$函数值逼近0。在坐标尺度比较大的情况下，$sigmoid(z)$函数很像一个阶跃函数。
因此，为实现Logistic回归分类器，我们可以在每个特征上都乘一个回归系数，然后把所有结果值相加后带入$sigmoid(z)$函数，任何大于0.5的数据被分为类别1，小于0.5的被归入类别0。
Python代码实现$sigmoid函数$:

def sigmoid(inX):
    return 1 / (1 + np.exp(-inX))

确定了分类器的函数形式之后，怎么确定最佳回归系数？？？

利用优化算法确定最佳回归系数

$sigmoid$函数的输入记为$z$，则：
$z=w_0{x}_0+w_1{x}_1+w_2{x}_2+...+w_n{x}_n+$
采用向量形式可以写为：$z=w^{T}x$，其中向量$w$就是我们需要找的最佳回归系数。
那么怎么求出最佳回归系数$w$呢？？？

梯度上升法求$w$

梯度上升法的思想：要找到某函数的最大值，最好的方法是沿着该函数的梯度方向探寻。如果梯度记为$\nabla$，则函数$f(x,y)$的梯度可以表示为：
$\nabla f(x,y)={⟮}_{\frac{\delta f(x,y)}{\delta y}}^{\frac{\delta f(x,y)}{\delta x}}⟯$，
这个公式意味着沿$x$的方向移动$\frac{\delta f(x,y)}{\delta x}$，沿$y$的方向移动$\frac{\delta f(x,y)}{\delta y}$，函数$f(x,y)$要在待计算的点上有定义并且可微。如图（图来自机器学习实战一书）：

从图中可以看出梯度算子总是指向函数值增长最快的方向，移动量的大小称为步长$\alpha$。则梯度的迭代公式可以表示为：
$w:=w$$+$$\alpha$$\nabla$${}_w$$f(w)$
该公式在达到某个停止条件之前一直迭代执行，例如迭代次数达到某个指定值或算法达到某个可以允许的误差范围。

图片来自机器学习实战一书。
梯度上升法的伪代码如下：

'''
1. 初始化回归系数w
2. for N次数：
3.     计算整个数据集的梯度gradient
4.     使用w = w + alpha*gradient更新回归系数
5.     返回回归系数
'''

python代码实现:

def grad_ascent(x_train, y_train, learning_rate=0.001, epochs=500):
    x_mat = np.matrix(x_train)
    y_mat = np.matrix(y_train).transpose()
    rows, cols = x_mat.shape
    args = np.ones((cols + 1, 1))
    add_1 = np.ones((rows,))
    x_mat = np.column_stack((x_mat, add_1))
    for k in range(epochs):
        y_pre = sigmoid(x_mat * args)
        loss = y_mat - y_pre
        args = args + learning_rate * x_mat.transpose() * loss
    w = args[:-1]
    b = np.array(args[-1])
    return w, b

画出决策边界的代码：

def plot_best_fit(x_train, y_train, weight, bias):
    import matplotlib.pyplot as plt
    rows = x_train.shape[0]
    x0 = []; x1 = []; y0 = []; y1 = []
    for i in range(rows):
        if y_train[i] == 0:
            x0.append(x_train[i, 0])
            y0.append(x_train[i, 1])
        else:
            x1.append(x_train[i, 0])
            y1.append(x_train[i, 1])
    ax = plt.subplot(111)
    ax.scatter(x0, y0, s=30, c='r')
    ax.scatter(x1, y1, s=30, c='g')
    x = np.arange(-3, 3, 0.1)
    y = (-bias - weight[0] * x) / weight[1]
    ax.plot(x, y.transpose())
    plt.xlabel('X1')
    plt.ylabel('Y1')
    plt.show()
    ```
加载数据集的代码：

```python
def read_data(data_path):
    data = np.loadtxt(data_path)
    return data[:, :-1], data[:, -1]

结果：

随机梯度上升法

梯度上升法在每次更新回归系数时都需要遍历整个数据集，该方法处理100个左右的数据还可以，但如果有数亿的样本和上万的特征，该方法的计算复杂度太高。因此将该算法改为以此仅用一个样本点更新回归系数，该方法称为随机梯度上升算法，伪代码和实现代码及结果如图(图来自机器学习实战一书)：



该方法的结果并非最佳拟合，这是因为上一个实现是在整个数据集上迭代500次才得到。判断一个优化算法优劣的可靠方法是看它是否收敛，也就是说参数是否达到稳定值。因此需要对随机梯度算法进行改进，代码如下：

def stoc_grad_ascent(x_train, y_train, learning_rate=0.01, epochs=500):
    rows, cols = x_train.shape
    args = np.ones((cols + 1, ))
    add_1 = np.ones((rows,))
    x_train = np.column_stack((x_train, add_1))
    for j in range(epochs):
        data_index = np.arange(rows)
        for i in range(rows):
            learning_rate = 4/(j+i+1)+0.01
            # learning_rate = learning_rate*0.95
            rand_index = int(np.random.uniform(0,len(data_index)))
            y_pre = sigmoid(sum(x_train[rand_index] * args))
            loss = y_train[rand_index] - y_pre
            args = args + learning_rate * x_train[rand_index] * loss
            np.delete(data_index, rand_index)
    w = args[:-1]
    b = args[-1]
    return w, b

改进之处：1) 每次迭代调整$\alpha$值；2) 通过随机取样更新回归系数。
其结果如图：