机器学习实战第五章 Logistic回归

Logistic回归

书上讲解的顺序其实并不是很好，下面是打乱顺序的讲解。

要解决的问题

逻辑回归是一种用于解决分类问题的统计学习方法。它主要用于预测二分类问题，即将输入数据分为两个不同的类别。

逻辑回归的目标是根据输入特征的线性组合，通过一个特定的函数（一般是sigmoid函数）将输入映射到一个概率值，该概率值表示样本属于某个类别的可能性。Sigmoid函数（也称为逻辑函数或Logistic函数）的形式为$$\sigma (z)=\frac{1}{1+e^{-z}}$$

• 这里插一句，为什么通常使用sigmoid函数呢？因为它的原函数是正态分布的，比较符合常见的事物的分布

一些更进一步的理解

如果说要提取重点，那么有以下这么几个点：

1. 实际上我们要做的事情还是计算特征的权重，然后根据权重去做下一步的事情
2. sigmoid函数给出的结果是，我们根据输入的特征预测的事物属于某一类（二分类）的概率
3. 根据sigmoid函数，我们能够更好的去计算损失（就是预测值与真实值之前的差距）并做进一步的优化

• 这里再插一句 上述第三点中提到了去做进一步优化，这一步通常使用的是梯度上升/下降算法。简而言之就是寻找能够更好的权重。下面我们将通过实验来介绍这一点。

实例

下面将采用书中附赠的代码来进行进一步的讲解，所有值得关注的点均已在代码中注释：

from numpy import *

# 加载数据集
def loadDataSet():
    dataMat = []
    labelMat = []
    fr = open('testSet.txt')
    for line in fr.readlines():
        lineArr = line.strip().split()
        # 在前面加了个1.0，是为了方便计算，同时由于是二维数据，即有两个自变量，所以有三个系数，1.0的意思是这里有一个点
        dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])
        labelMat.append(int(lineArr[2]))
    return dataMat, labelMat


# 定义sigmoid函数，也就是激活函数
def sigmoid(inX):
    return 1.0 / (1 + exp(-inX))


# 梯度上升算法
def gradAscent(dataMatIn, classLabels):
    # （以下两⾏）转换为NumPy矩阵数据类型
    dataMatrix = mat(dataMatIn)
    labelMat = mat(classLabels).transpose()
    
    # dataMatrix的⾏（m）和列（n）数
    m, n = shape(dataMatrix)
    # 变化率，也可以称为学习率、步长等，主要用于控制梯度上升速度
    alpha = 0.001
    # 最大迭代次数
    maxCycles = 500
    # 创建与列数相同的⾮零向量，初始化特征的权重值
    weights = ones((n, 1))
    # 开始迭代
    for k in range(maxCycles):
        # 计算当前选择样本的sigmoid值，也即样本属于1类的概率
        h = sigmoid(dataMatrix * weights)
        # 计算真实类别与预测类别的差值，这里的差值是指概率上的差值
        error = (labelMat - h)
        # 按照差值的⽅向调整回归系数，这里使用了梯度上升算法，就是找上升最快的方向
        # 最终的目的是为了找到最佳参数weights，使得概率最大。迭代这么多次是为了使得概率最大（error趋向于0），此时weights就是最佳参数，同时也会稳定
        weights = weights + alpha * dataMatrix.transpose() * error
    return weights

dataArr, labelMat = loadDataSet()
weights = gradAscent(dataArr, labelMat)
print(weights)

运行结果如下图所示

画图的代码此处不再给出

可以看到拟合的非常好。然而为什么我们将二维坐标转化为了三维坐标呢？假如说我们只用二维坐标，那么我们求解权重的过程就是求解
$$ax+by=0$$
不难发现，这个解出来的直线是必过原点的，然而这样的效果并不好，因此我们需要增加一个额外的变量来增加自由度。这样一来就变成了
$$ax+by+c=0$$
这样求出的解就能够很好的划分两类点了，这就是我们代码中将变量升维的原因。

随机梯度上升

上述代码+理论其实已经非常完整了，但是依旧存在些许缺陷，例如现在我们是100个样本就要转化为100*3的矩阵，假如后续有了更多的样本，那么计算量几乎是无法接受的。因此可以将代码实现改为每次只选一个样本更新回归系数。这就是随机梯度上升算法。

• 补充：像这种可以随时更新的称为在线学习，一次处理所有数据而无法随时更新的被称作批处理

随机梯度上升代码及调用如下：

# 随机梯度上升算法
def stocGradAscent0(dataMatrix, classLabels):
    m, n = shape(dataMatrix)
    alpha = 0.01
    weights = ones(n)
    # 主要区别在于，每次只用一个样本更新，之前是用所有样本更新500次
    for i in range(m):
        h = sigmoid(sum(dataMatrix[i] * weights))
        error = classLabels[i] - h
        weights = weights + alpha * error * dataMatrix[i]
    return weights

dataArr, labelMat = loadDataSet()
weights = stocGradAscent0(array(dataArr), labelMat)
print(weights)
plotBestFit(weights)

如果看样本使用次数的话，当前的代码每个样本只用了1次，而之前的代码则是用了500次，那么效果差是显而易见的：

有很多都分类错了地方，如果对代码进行一些小小的改动：

# 随机梯度上升算法
def stocGradAscent0(dataMatrix, classLabels):
    m, n = shape(dataMatrix)
    alpha = 0.01
    # ❶（以下三⾏）创建与列数相同的⾮零向量
    weights = ones(n)
    for j in range(500):
        for i in range(m):
            # ❷（以下两⾏）计算当前选择样本的sigmoid值
            h = sigmoid(sum(dataMatrix[i] * weights))
            error = classLabels[i] - h
            # ❸（以下两⾏）更新回归系数
            weights = weights + alpha * error * dataMatrix[i]
    return weights

那么效果就可以达到跟之前几乎同样优秀的水准：

此外值得一提的是，逻辑回归需要进行大量计算才能逐渐收敛。如下图：
而且即使收敛了，也会有可能在小范围内震荡：

最后依旧震荡的原因是，不论我们怎么分类，总会有一些点无法正确划分（这是数据与边界线类型的选择共同造成的结果），此时可以对随机梯度上升算法做进一步优化来加快收敛速度并减少波动。

改进随机梯度上升

为了加快收敛，并减少波动，我们可以对其步长做改动：

# 改进的随机梯度上升算法
def stocGradAscent1(dataMatrix, classLabels, numIter=150):
    m, n = shape(dataMatrix)
    weights = ones(n)
    for j in range(numIter):
        dataIndex = list(range(m))
        for i in range(m):
            # 随着迭代进行，逐渐减小步长
            alpha = 4 / (1.0 + j + i) + 0.01
            # 随机选取更新，能够减少周期性波动
            randIndex = int(random.uniform(0, len(dataIndex)))
            h = sigmoid(sum(dataMatrix[randIndex] * weights))
            error = classLabels[randIndex] - h
            weights = weights + alpha * error * dataMatrix[randIndex]
            del (dataIndex[randIndex])
    return weights

总的来说，这个方法的主要思想就是：

1. 控制步长
2. 随机选取数据来改善周期性波动
   其中我觉得第一点相对而言更重要一些，对比图如下：
   
   可以看到二者的区别其实并不太大，只是后者周期性波动的问题没有得到解决，这应该是因为我们没有将控制步长变换的最小值控制到更小的范围，如果可以的话应该能够进一步缓解周期性波动的问题。

从疝气病症预测病马的死亡率

除了数据集需要预处理，此处并没有什么特殊的地方。又因为我们直接使用了随书附赠的干净数据集，因此直接给出代码如下：

# ❶（以下两⾏）打开训练集和测试集
def classifyVector(inX, weights):
    prob = sigmoid(sum(inX * weights))
    if prob > 0.5:
        return 1.0
    else:
        return 0.0

def colicTest():
    frTrain = open('horseColicTraining.txt')
    frTest = open('horseColicTest.txt')
    # ❷（以下两⾏）训练集和测试集的列表
    trainingSet = [];
    trainingLabels = []
    # ❸（以下两⾏）训练集的逐行读取
    for line in frTrain.readlines():
        currLine = line.strip().split('\t')
        # ❹（以下两⾏）训练集的特征值
        lineArr = []
        for i in range(21):
            lineArr.append(float(currLine[i]))
        # ❺（以下两⾏）训练集的标签值
        trainingSet.append(lineArr)
        trainingLabels.append(float(currLine[21]))
    # ❻（以下两⾏）计算回归系数向量
    trainWeights = stocGradAscent1(array(trainingSet), trainingLabels, 500)
    # ❼（以下两⾏）测试集的逐行读取
    errorCount = 0;
    numTestVec = 0.0
    for line in frTest.readlines():
        numTestVec += 1.0
        currLine = line.strip().split('\t')
        # ❽（以下两⾏）测试集的特征值
        lineArr = []
        for i in range(21):
            lineArr.append(float(currLine[i]))
        # ❾（以下两⾏）测试集的标签值
        if int(classifyVector(array(lineArr), trainWeights)) != int(
                currLine[21]):
            errorCount += 1
    # ❿（以下⼏⾏）打印错误率
    errorRate = (float(errorCount) / numTestVec)
    print("the error rate of this test is: %f" % errorRate)
    return errorRate

# ❶（以下两⾏）调⽤函数10次并求平均值
def multiTest():
    numTests = 10;
    errorSum = 0.0
    for k in range(numTests):
        errorSum += colicTest()
    print("after %d iterations the average error rate is: %f" %
          (numTests, errorSum / float(numTests)))

multiTest()

小结

开始逐渐学到了梯度下降、收敛等概念。本章中并没有直接给出梯度下降的定义，因为梯度就是变化速度最快的方向，因此算出来前后结果的差值之后直接用就行了。至于收敛可以通过控制步长来实现，此处的数据集比较简单，没有太多可以说的地方。