机器学习：逻辑回归

虽然名字里带回归，但实际上是一种分类方法，主要用于两分类问题，即只有两种分类

优点：计算代价不高，易于理解和实现
缺点：容易欠拟合，分类精度可能不高

原理

线性回归函数
  (small z = f(X) = XW)
  其中
    X 是特征值
    W 是回归系数
  X 和 W 都是向量，可展开为
    (small z = XW = X_{0}W_{0} + X_{1}W_{1} + ... + X_{n}W_{n})
  线性方程其实应该是
    (small z = XW + b)
  为此这里固定
    (small X_{0}=1)
    (small W_{0}=b)
  其他 X 值才是用户输入，这样变成两个向量相乘方便计算

逻辑回归函数 (Sigmoid 函数)
  
  (large y=g(z)=frac{1}{1+e^{-z}})
  
  该函数模拟阶跃函数 (在某个跳跃点从 0 瞬间跳到 1，跳跃点两边的值固定为 0 和 1)
  可以得出
    (small y=left{egin{matrix}0.12&z=-2\0.5&z=0\0.88&z=2end{matrix} ight.)
  且满足
    (small g(z) + g(-z) = 1)
  
  在 z 轴比较长的情况下看起来就像跳跃点为 0 的阶跃函数

  

分类器
  结合线性回归函数和逻辑回归函数
  
    (large y=g(z)=frac{1}{1+e^{-z}}=frac{1}{1+e^{-XW}})
  
  对特征向量 X 进行计算，得出 0~1 之间的值
  大于 0.5 属于分类 1，小于 0.5 属于分类 0
  所以逻辑回归也可以被看成是一种概率估计

训练分类器
  通过最优化算法（通常是梯度下降算法），寻找最佳回归系数 W 的值

梯度下降算法

要找到某函数的最小值，最好的方法是沿着该函数的梯度方向向下探寻
就像等高线，选择最高点的最快方法就是不断沿着梯度走
  
    
  
逻辑回归函数 (small g(z)) 的结果实际是一个 (small (0,1)) 区间的概率估计
为此损失函数采取逻辑损失函数
  
  (large L = -(y)log(g(z)) - (1-y)log(1-g(z)))
  
当 (small y = 0) 时，(small L = -log(1-g(z)))，(small g(z))范围是(small (0,1))，对应 L 范围是 (small (0,infty))
当 (small y = 1) 时，(small L = -log(g(z)))，(small g(z))范围是(small (0, 1))，对应 L 范围是 (small (infty, 0))

梯度下降算法的目的就是找最佳的 (small w) 参数，使得所有样本的 (small L) 值的总和最小
  
  (large L(w) = - sum_{i=0}^{n}( (y)log(g(z)) + (1-y)log(1-g(z)) ))
  
梯度算法的迭代公式如下
  (large w := w - alpha abla L(w))
其中
  (small abla L(w)) 是函数 (small L) 在 (small w) 处的梯度，代表使 (small L) 值变化率最大的方向
  (alpha) 是步长，即沿梯度方向变化的大小，必须取一个很小的值

该公式一直被迭代执行，即(small w) 沿着梯度方向不断减少，使(small L)值以最快的速率不断下降
直至达到某个停止条件，如迭代次数达到阈值，或(small L)值达到某个可以允许的误差范围
进一步了解梯度，参考 导数、偏导数、方向导数、梯度、梯度下降
  
已知对下列函数求导
  ( ormalsize g(z) = frac{1}{1+e^{-z}})
  
  ( ormalsize f(x) = log(x))
可得
  (large frac{partial g(z)}{partial z}=g(z)(1-g(z)))
  
  (large frac{partial f(x)}{partial x}=frac{1}{x})
  
梯度在每个轴上的分量是函数在该轴的偏导数
  
  (large abla L(w) = (frac{partial L(w)}{partial w_{0}}, ... , frac{partial L(w)}{partial w_{j}}, ... , frac{partial L(w)}{partial w_{m}}))
  
结合以上内容求解 (small L(w)) 的偏导数（i 代表样本，共 n 个，j 代表特征，共 m 个）
 
  (large frac{partial L(w)}{partial w_{j}}=-sum_{i=0}^{n}(frac{partial((y)log(g(z)))}{partial w_{j}}+frac{partial((1-y)log(1-g(z)))}{partial w_{j}}))
  
      (large = -sum_{i=0}^{n}(yfrac{1}{g(z)}frac{partial g(z)}{partial w_{j}}-(1-y)frac{1}{1-g(z)}frac{partial g(z)}{partial w_{j}}))
  
      (large = -sum_{i=0}^{n}((frac{y}{g(z)}-frac{(1-y)}{1-g(z)})frac{partial g(z)}{partial w_{j}}))
  
      (large = -sum_{i=0}^{n}((frac{y}{g(z)}-frac{(1-y)}{1-g(z)})g(z)(1-g(z))frac{partial z}{partial w_{j}}))
  
      (large = -sum_{i=0}^{n}((y(1-g(z))-(1-y)g(z))frac{partial xw}{partial w_{j}}))
  
      (large = -sum_{i=0}^{n}((y-g(z))frac{partial xw}{partial w_{j} }))
  
      (large = -sum_{i=0}^{n}((y-g(z))x_{j}))
  
      (large = -sum_{i=0}^{n}(e^{(i)}x_{j}^{(i)}))
  
      (large = -egin{bmatrix}x_{j}^{(0)}&...&x_{j}^{(i)}&...&x_{j}^{(n)}end{bmatrix}egin{bmatrix}e^{0}\...\e^{i}\...\e^{n} end{bmatrix})
 
将所有偏导组成梯度
 
   (large abla L(w) = (frac{partial L(w)}{partial w_{0}}, ... , frac{partial L(w)}{partial w_{j}}, ... , frac{partial L(w)}{partial w_{m}}))
  
      (large = -egin{bmatrix}x_{0}^{(0)}&...&x_{0}^{(i)}&...&x_{0}^{(n)}\x_{j}^{(0)}&...&x_{j}^{(i)}&...&x_{j}^{(n)}\x_{m}^{(0)}&...&x_{m}^{(i)}&...&x_{m}^{(n)}end{bmatrix}egin{bmatrix}e^{0}\...\e^{i}\...\e^{n} end{bmatrix})
  
      (large = -X^{T}E)
 
进一步求得梯度下降公式
 
  (large W = W - alpha abla L(W))
  
     (large = W + alpha X^{T}E)
  
     (large = W + alpha X^{T}(Y-g(XW)))
 
其中
  样本集有 n 条数据，m 个特征
  X 是 (small (n，m)) 矩阵
  W 是 (small (m，1)) 矩阵
  Y 是 (small (n，1)) 矩阵
  (small alpha > 0)

随机梯度下降（SGD - Stochastic Gradient Descent）

可以看到，每一次迭代都要对所有样本计算，计算量太大
可以改为每一次迭代都随机取一部分样本计算就可以

代码

# coding=utf-8
import numpy as np
import random


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


# 梯度下降算法
# 根据样本（X,Y) 算法最佳的 W
def gradDescent(sampleData, classLabels):
    """
    sampleData  - 样本特征，(n,m) 的二维数组，n 是样本数，m 是特征数
                  每行的第一个值 X0 固定为 1，从 X1 开始才是真正的特征值，目的是简化向量的计算
                   y = x1*w1 + ... + xm*wm + b
                     = x1*w1 + ... + xm*wm + x0w0
                     = XW

    classLabels - 样本标签，(1,n) 的一维数组
    """

    # 转为 NumPy 矩阵
    dataMatrix = np.mat(sampleData)

    # 将 (1,n) 转为 (n,1) 方便后面的矩阵计算
    labelMatrix = np.mat(classLabels).transpose()

    # n 个样本，m 个特征值
    n, m = np.shape(dataMatrix)

    # 梯度下降的步长
    alpha = 0.001

    # 最大迭代次数
    maxCycles = 500

    # 初始化 W 为 (m,1) 数组, 默认值为 1
    weights = np.ones((m, 1))

    # 迭代
    for k in range(maxCycles):
        # (n,m) 矩阵乘以 (m,1) 矩阵，得到 (n,1) 矩阵，再通过逻辑回归函数得到样本的 Y
        h = sigmoid(dataMatrix*weights)

        # 两个 (n,1) 矩阵，得到每个样本的误差
        error = (labelMatrix - h)

        # w = w + a*(X^T)*(Y-g(XW))
        #   = w + a*(X^T)*E
        weights = weights + alpha * dataMatrix.transpose() * error

    return weights


# 分类
def classify(data, weights):
    dataMatrix = np.mat(data)
    resultMatrix = sigmoid(dataMatrix * weights) > 0.5
    
    """
    for result in resultMatrix:
        print(result.item())
    """
    return resultMatrix


# 随机梯度
def stocGradDescent0(sampleData, classLabels):
    dataMatrix = np.mat(sampleData)
    labelMatrix = np.mat(classLabels).transpose()

    n, m = np.shape(dataMatrix)
    alpha = 0.01
    weights = np.ones((m, 1))

    for i in range(n):
        # 每次迭代只取一个样本，迭代次数为样本个数
        h = sigmoid(dataMatrix[i]*weights)
        error = labelMatrix[i] - h
        weights = weights + alpha * dataMatrix[i].transpose() * error

    return weights


# 改进的随机梯度
def stocGradDescent1(sampleData, classLabels, numIter=150):
    dataMatrix = np.mat(sampleData)
    labelMatrix = np.mat(classLabels).transpose()

    n, m = np.shape(dataMatrix)
    weights = np.ones((m, 1))

    # 自己选择迭代次数
    for j in range(numIter):
        dataIndex = range(n)

        # 每次迭代又迭代了每一个样本
        for i in range(n):
            # 每次迭代都改变步长，0.0001 用于防止出现 0 的情况
            alpha = 4/(1.0+j+i)+0.0001

            # 随机选择一个样本
            randIndex = int(random.uniform(0, len(dataIndex)))

            h = sigmoid(dataMatrix[randIndex]*weights)
            error = labelMatrix[randIndex] - h

            # 计算新的 W
            weights = weights + alpha * dataMatrix[randIndex].transpose() * error

            # 删除该样本下标
            del(dataIndex[randIndex])

    return weights