ML：Logistic Regression

• 逻辑回归的主要思想：根据现有数据对分类边界线建立回归公式，以此进行分类
• 优点：计算代价不高，易于理解和实现
• 缺点：容易欠拟合，分类精度可能不高
• 适用数据类型：数值型和标称型数据

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

• 逻辑回归需要的函数应该是，能够接受所有的输入，然后预测出类别，对于二分类，输出为 0 或 1，这种函数为单位阶跃函数（海维赛德阶跃函数），该函数在跳跃点上瞬间从 0 跳到 1，这一点在数学上很难处理（不可导）
• Sigmoid 函数具有阶跃函数类似的性质，且是连续的（处处可导），该函数在 0 附近变化很快，在横坐标尺度足够大的情况下，很像阶跃函数，其计算公式如下$$\sigma (z)=\frac{1}{1+e^{-z}}$$
• 在每个特征上都乘以一个回归系数，然后把所有结果值相加，将这个总和代入 Sigmoid 函数中，进而得到一个范围在 0~1 之间的数值，以 0.5 为分界线进行分类
• 接下来的问题就是，如何求解最佳回归系数？

2. 基于最优化方法的最佳回归系数确定

• Sigmoid 函数的输入记为 z，其计算公式为：
  $$z=w_0{x}_0+w_1{x}_1+w_2{x}_2+...+w_n{x}_n$$
  上述公式也可以采用向量的写法，写为 $z={\mathbf{w}}^T\mathbf{x}$，结合 Sigmoid 函数，输出 $\hat{y}$：
  $$\hat{y}=\frac{1}{1+e^{-z}}=\frac{1}{1+e^{-{\mathbf{w}}^T\mathbf{x}}}$$
  $\hat{y}$ 为我们的预测输出，其范围是(0, 1)，假设其为 y 取 1 的概率，即 $p(y=1|\mathbf{x},\mathbf{w})=\hat{y}$，则 $p(y=0|\mathbf{x},\mathbf{w})=1-\hat{y}$，当样本实际类别为 1，即 y = 1 的时候，$\hat{y}$ 越大越好，当样本实际类别为 0，即 y = 0 的时候，$1-\hat{y}$ 越大越好，设成本函数(Cost Function)：
  $$Cost(y,\hat{y})={\hat{y}}^y(1-\hat{y}{)}^{1-y}$$
  y = 1 时，$Cost(y,\hat{y})=\hat{y}$，当 y = 0 时，$Cost(y,\hat{y})=1-\hat{y}$，为了简化计算，将成本函数求自然对数：
  $$Cost(y,\hat{y})=yln\hat{y}+(1-y)ln(1-\hat{y})$$
  上述公式只是对于单个样本而言，对于整个样本集，若样本间相互独立，整个样本集的概率等于各样本概率乘积，乘积取自然对数为累加，得到：
  $$J(\mathbf{w})=\sum _{i=1}^m\left[y^{(i)}ln\hat{y}(\mathbf{w},{\mathbf{x}}^{(i)})+(1-y^{(i)})ln(1-\hat{y}(\mathbf{w},{\mathbf{x}}^{(i)}))\right]$$因此我们的问题转化为求解 $J(\mathbf{w})$ 的最大值
• FYI：这里我们需要求最大值，因此使用梯度上升法，如果是求最小值，如$J_2(\mathbf{w}):=1-J(\mathbf{w})$ ，则使用梯度下降法，原理都是一样的。

2.1 梯度上升法

• 基本思想：要找到某函数的最大值，最好的方法是沿着该函数的梯度方向探寻
• 对于函数 $f(x,y)$，其梯度：
  $$\nabla f(x,y)=\left(\begin{array}{c}\frac{\partial f(x,y)}{\partial x}\\ \frac{\partial f(x,y)}{\partial y}\end{array}\right)$$
  该向量的方向即为移动方向，因为梯度算子总是只想函数值增长最快的方向。我们使用迭代的方式来探寻函数最大值，每次都沿着梯度方向移动，而移动量的大小称为步长，记作 $\alpha$，对于$J(\mathbf{w})$，迭代公式如下：
  $$\mathbf{w}:=\mathbf{w}+\alpha {\nabla }_{\mathbf{w}}J(\mathbf{w})=\mathbf{w}+\alpha \frac{\partial J(\mathbf{w})}{\partial \mathbf{w}}$$
  对于 $w_j$:
  $$w_j:=w_j+\alpha \frac{\partial J(\mathbf{w})}{\partial w_j}$$
  因此我们需要求解 $J(\mathbf{w})$ 的偏导数，由前面知道：
  $$J(\mathbf{w})=\sum _{i=1}^m\left[y^{(i)}ln{\hat{y}}^{(i)}+(1-y^{(i)})ln(1-{\hat{y}}^{(i)})\right]$$
  而
  $${\hat{y}}^{(i)}=\frac{1}{1+e^{-z^{(i)}}}$$
  $$z^{(i)}={\mathbf{w}}^T{\mathbf{x}}^{(i)}$$
  对于$J(\mathbf{w})$ 的偏导数：
  $$\frac{\partial J(\mathbf{w})}{\partial w_j}=\frac{\partial J(\mathbf{w})}{\partial {\hat{y}}^{(i)}}\frac{\partial {\hat{y}}^{(i)}}{\partial z^{(i)}}\frac{\partial z^{(i)}}{\partial w_j}$$
  其中：
  $$\frac{\partial J(\mathbf{w})}{\partial {\hat{y}}^{(i)}}=\sum _{i=1}^m\left[\frac{y^{(i)}}{{\hat{y}}^{(i)}}+\frac{y^{(i)}-1}{1-{\hat{y}}^{(i)}}\right]$$
  $$\frac{\partial {\hat{y}}^{(i)}}{\partial z^{(i)}}=\frac{\partial }{\partial z^{(i)}}\left(\frac{1}{1-e^{-z^{(i)}}}\right)=\frac{-e^{-z^{(i)}}(-1)}{(1-e^{-z^{(i)}}{)}^2}=\frac{1}{1-e^{-z^{(i)}}}\frac{-e^{-z^{(i)}}}{1-e^{-z^{(i)}}}=\frac{1}{1-e^{-z^{(i)}}}\frac{1}{1-e^{z^{(i)}}}={\hat{y}}^{(i)}(1-{\hat{y}}^{(i)})$$
  $$\frac{\partial z^{(i)}}{\partial w_j}=\frac{\partial }{\partial w_j}(w_0{x}_0^{(i)}+w_1{x}_1^{(i)}+w_2{x}_2^{(i)}+...+w_n{x}_n^{(i)})=x_j^{(i)}$$
  三式相乘：
  $$\frac{\partial J(\mathbf{w})}{\partial w_j}=\sum _{i=1}^m\left[\frac{y^{(i)}}{{\hat{y}}^{(i)}}+\frac{y^{(i)}-1}{1-{\hat{y}}^{(i)}}\right]{\hat{y}}^{(i)}(1-{\hat{y}}^{(i)})x_j=\sum _{i=1}^n\left[(y^{(i)}-{\hat{y}}^{(i)})x_j^{(i)}\right]$$
  由此得到迭代公式：
  $$w_j:=w_j+\alpha \sum _{i=1}^m\left[(y^{(i)}-{\hat{y}}^{(i)})x_j^{(i)}\right]$$

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

'''Logistic 回归梯度上升优化算法'''

'''载入数据集，返回特征矩阵、标签列表'''
def loadDataSet():
    dataMat = []
    labelMat = []
    with open('Ch05/testSet.txt') as f:
        for line in f.readlines():
            lineArr = line.strip().split()
            dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])    # x0 = 1
            labelMat.append(int(lineArr[2]))
    return dataMat, labelMat

'''Sigmoid 函数'''
def sigmoid(inX):
    return 1 / (1 + np.exp(-inX))

'''梯度上升算法，返回权重值'''
def gradAscent(dataMatIn, classLabels):
    dataMatrix = np.mat(dataMatIn)    # shape(dataMatrix) = (m, n)
    labelMat = np.mat(classLabels).T    # shape(labelMat) = (m, 1)
    m, n = np.shape(dataMatrix)
    alpha = 0.001
    maxCycles = 500
    weights = np.ones((n, 1))    # shape(weights) = (n, 1)
    for k in range(maxCycles):
        h = sigmoid(dataMatrix * weights)   # shape(h) = (m. n) * (n, 1) = (m, 1)
        error = (labelMat - h)
        weights = weights + alpha * dataMatrix.T * error
    return weights

[In]: dataArr, labelMat = loadDataSet()
[In]: gradAscent(dataArr, labelMat)
[Out]: matrix([[ 4.12414349],
        	   [ 0.48007329],
        	   [-0.6168482 ]])

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

def plotBestFit(weights):
    import matplotlib.pyplot as plt
    dataMat, labelMat = loadDataSet()
    dataArr = np.array(dataMat)
    n = dataArr.shape[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', label='1')
    ax.scatter(xcord2, ycord2, s=30, c='green', label='0')
    x = np.arange(-3., 3., 0.1)
    y = (-weights[0] - weights[1] * x) / weights[2]
    ax.plot(x, y.T)
    plt.xlabel('X1')
    plt.ylabel('X2')
    plt.legend()
    plt.show()

2.4 训练算法：随机梯度上升

• 上述梯度上升算法有个问题，它需要在每次更新回归系数时都遍历整个数据集，这一点在数据集很大的情况下，计算复杂度太高，对应的改进方法为每次只用一个样本来更新回归系数，称为随机梯度上升算法
• 由于随机梯度上升算法可以在新样本到来时对分类器进行增量式更新，因此该算法是一个在线学习算法。相应地，一次处理所有数据被称作批处理

'''训练算法：随机梯度上升'''
def stocGradAscent0(dataMatrix, classLabels):
    m, n = np.shape(dataMatrix)
    alpha = 0.01
    weights = np.ones(n)
    for i in range(m):
        h = sigmoid(sum(dataMatrix[i] * weights))
        error = classLabels[i] - h
        weights = weights + alpha * error * dataMatrix[i]
    return weights

• 随机梯度上升算法中，h 和 error 都是数值
• 从结果可以看到，上面的随机梯度上升算法结果与普通梯度上升相比还是有差距。然而判断一个算法优劣的可靠方法是看它是否收敛，也就是参数是否达到稳定值。上面的随机梯度上升算法在整个数据集上运行200次，三个回归系数的变化情况如下图，可以看到收敛状况也不是很好，而且还存在周期性波动（主要是因为一些不能正确分类的样本点导致每次迭代时系数都会发生改变）
  
• 对于上述问题，可以进行两处改进，一是对于 alpha 的取值，可以每次迭代都进行调整，从而减缓数据波动。alpha 随着迭代次数不断减小，我们设置一个常数项，使其不会减小到 0。每次迭代 alpha 减小 1/(i+j)，随着迭代次数的增加，调整值也在减小。第二处改进是每次随机选择样本，可以减少周期性波动

'''改进的随机梯度上升算法'''

def stocGradAscent1(dataMatrix, classLabels, numIter = 40):
    m, n = np.shape(dataMatrix)
    weights = np.ones(n)
    for j in range(numIter):
        for i in range(m):
            alpha = 4 / (1.0 + j + i) + 0.01
            randIndex = int(np.random.uniform(0, m))
            h = sigmoid(sum(dataMatrix[randIndex] * weights))
            error = classLabels[randIndex] - h
            weights = weights + alpha * error * dataMatrix[randIndex]
            np.delete(dataMatrix, randIndex)
    return weights

• 我们迭代了 40 次的结果如上，可见其效果和普通梯度上升算法相近，且只迭代了 40 次，而后者需要迭代 500 次。从回归系数的变化来看，改进后的算法可以收敛得更快，而且没有了周期性波动