ML | 5 Logistic回归

文章目录

ML | 5 Logistic回归
- Logistic回归思想
- 基于Logistic回归和Sigmoid函数分类
- 基于最优化方法的最佳回归系数确定
- - 梯度上升法
  - - **梯度上升法例子：**
    - 梯度下降法
  - 训练算法：使用梯度上升找到最佳参数
  - - 画出决策边界
  - 训练算法：随机梯度上升
- 示例：与疝气病症预测病马的死亡率
- - 处理数据中的缺失值
  - 测试算法
- 欢迎关注公众号【三戒纪元】

Logistic回归思想

根据现有数据对分类边界线建立回归公式，以此进行分类。

“回归”源于最佳拟合，表示要找到最佳拟合参数集。

训练分类器时的做法就是寻找最佳拟合参数，使用的是最优化算法。

优点：计算代价不高，易于理解和实现
缺点：容易欠拟合，分类精度可能不高
适用数据类型：数值型和标称型数据

基于Logistic回归和Sigmoid函数分类

我们想要的函数是，能够接受所有的输入，然后预测出类别。

例如，在2个类别的情况下，上述函数输入0或1，这种函数成为海维赛德阶跃函数（Heaviside step function）或者直接称为单位阶跃函数

Sigmoid函数有类似性质
$\sigma(z)=\frac{1}{1+e^{-z}}$

上图为2种尺度下的Sigmoid函数图。上图横坐标为-5到5，曲线变化较为平滑；下图横坐标尺度足够大，Sigmoid函数看起来很像阶跃函数

为了实现Logistic回归分类器，可以在每个特征上都乘以一个回归系数，然后把所有结果值相加，将这个总和带入Sigmoid函数中，进而得到一个范围在0~1之间的数值。任何大于0.5的数据被分入1类，小于0.5的被归为0类。所以 Logistic 回归也可以被看成一种概率估计。

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

Sigmoid 函数的输入记为 z，由下面公式得出：
$w_0x_0 +w_1x_1+w_2x_2+ \ldots +w_nx_n$
采用向量法，记为 $z = w^Tx$

其中, x 是分类器的输入数据，向量 w 就是最佳参数（系数），从而使得分类器尽可能精确

梯度上升法

梯度上升法基于的思想是:要找到某函数的最大值，最好的方法是沿着该函数的梯度方向探寻

如果梯度记为 $\nabla$ ，则函数 f(x,y) 的梯度表示为：
$\nabla f(x,y) = \begin{pmatrix} \frac{\partial f(x,y)}{\partial x} \\ \frac{\partial f(x,y)}{\partial y} \end{pmatrix}$
这个梯度意味着：要沿x的方向移动 $\frac{\partial f(x,y)}{\partial x}$ ，沿y 的方向移动 $\frac{\partial f(x,y)}{\partial y}$ ，其中，函数f(x,y)必须在待计算点上有定义且可微。

梯度上升法例子：

梯度上升算法到达每个点后都会重新估计移动方向。

从P0开始，计算完该点的梯度，函数会根据梯度移动到下一点P1。

在P1点，梯度再次被重新计算，并沿新的梯度方向移动到P2.

如此循环迭代，直到满足停止条件。

迭代的过程中，梯度算子总是保证我们能选取到最佳的移动方向。

记每次移动量为步长，记为 $\alpha$ ，用向量表示梯度上升法的迭代公式如下：
$w:=w+\alpha \nabla_w f(w)$

梯度下降法

梯度下降法与梯度上升法是一样的，只不过公式中的加法需要变成减法：
$w:=w-\alpha \nabla_w f(w)$

梯度上升法用来求函数的最大值，梯度下降法用来求函数的最小值。

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

下图有100个样本点，每个点包含2个数值型特征：X1 和 X2。

在此数据集上，我们将通过使用梯度上升法找到最佳回归系数，也就是拟合出 Logistic 回归模型的最佳参数。

梯度上升伪代码：

每个回归系数初始化为1

重复R次：

计算整个数据集的梯度

使用 $\alpha \times gradient$ 更新回归系数的向量

返回回归系数

编写一个 logRegres.py的回归梯度上升优化算法

from numpy import *

# 加载 100个样本的简单数据集
def loadDataSet():
    dataMat = []; labelMat = []
    fr = open('testSet.txt')
    for line in fr.readlines():
        lineArr = line.strip().split()
        dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])]) # 包含3个特征 X0 X1， X2, X0=1
        labelMat.append(int(lineArr[2]))
    return dataMat,labelMat

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

def gradAscent(dataMatIn, classLabels):
    dataMatrix = mat(dataMatIn)             #convert to NumPy matrix 100x3
    labelMat = mat(classLabels).transpose() #convert to NumPy matrix  100x1
    m,n = shape(dataMatrix) # m: 100   n: 3
    alpha = 0.001 # 向目标移动的步长
    maxCycles = 500 # 迭代次数
    weights = ones((n,1))  # 3x1 [1, 1, 1]
    # 以下是对回归系数的500次迭代
    for k in range(maxCycles):              #heavy on matrix operations
        h = sigmoid(dataMatrix*weights)     #matrix mult 100x1 每个值范围在 0-1
        error = (labelMat - h)              #vector subtraction 100x1
        weights = weights + alpha * dataMatrix.transpose()* error #matrix mult 3x1
    return weights 
 # [[ 4.12414349]
 #  [ 0.48007329]
 #  [-0.6168482 ]]

核心代码是

    # 以下是对回归系数的500次迭代
for k in range(maxCycles):              #heavy on matrix operations
    h = sigmoid(dataMatrix*weights)     #matrix mult 100x1 每个值范围在 0-1
    error = (labelMat - h)              #vector subtraction 100x1
    weights = weights + alpha * dataMatrix.transpose()* error #matrix mult 3x1

计算真实类别与预测类别的差值，按照差值的方向调整回归系数

画出决策边界

def plotBestFit(weights):
    import matplotlib.pyplot as plt
    dataMat,labelMat=loadDataSet()
    dataArr = array(dataMat)
    n = shape(dataArr)[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')
    ax.scatter(xcord2, ycord2, s=30, c='green')
    x = arange(-3.0, 3.0, 0.1)
    y = (-weights[0]-weights[1]*x)/weights[2] # 最佳拟合直线 转换即可得 w[0] + w[1]x + w[2]y = 0
    ax.plot(x, y)
    plt.xlabel('X1'); plt.ylabel('X2');
    plt.show()

需要指出，此处设置了 sigmoid 函数为0，因为0个两个分类的分界处(类别0和类别1），因此设定 $0 = w_0x_0+w_1x_1+w_2x_2$ ，其中 $x 0 = 1$ ，然后解出 X1和X2的关系式

训练算法：随机梯度上升

梯度上升算法在每次更新回归系数时都需要遍历整个数据集，如果数据集特征很多，计算复杂度就太高了。

改进办法是一次仅用1个样本点更新回归系数，该方法称为随机梯度上升算法。

由于可以在新样本到来前对分类器进行增量更新，因此随机梯度上升算法是一个在线学习算法。

随机梯度上升算法伪代码：

所有回归系数初始化为1

对数据集中每个样本

计算该样本的梯度

使用 $\alpha \times gradient$ 更新回归系数的值

返回回归系数数值

具体实现代码为：

def stocGradAscent0(dataMatrix, classLabels):
    m,n = shape(dataMatrix)
    alpha = 0.01
    weights = ones(n)   #initialize to all ones
    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 全是数值，不再是向量；计算过程没有矩阵转换过程，所有变量的数据类型都是NumPy数组

使用该方法分割上述点数据：

    weights = stocGradAscent0(array(dataArr), labelMat)
    print("weight: ", weights) # [ 1.01702007  0.85914348 -0.36579921]
    plotBestFit(weights)

该算法迭代了100次，计算出来的回归系数的效果不如前1个方法迭代500次的效果好。

所以回归系数收敛需要一定的迭代次数。

当随机梯度上升算法迭代200次，可以观察下系数值X0，X1，X2的变化情况

X2经过50次后就达到了稳定值，而X0和X1需要更多的迭代。

有时候样本中存在一些不能正确分类的样本点（样本集非线性可分），导致系数值会有周期性波动，所以每次迭代会引发系数的剧烈改变。

所以算法中最好能够避免来回波动，从而收敛到某个值，另外，收敛速度也要加快。

修改后的随机梯度上升算法：

def stocGradAscent1(dataMatrix, classLabels, numIter=150):
    m,n = shape(dataMatrix)
    weights = ones(n)   # 初始化为1
    for j in range(numIter):
        dataIndex = list(range(m))
        for i in range(m):
            alpha = 4/(1.0+j+i)+0.01    #alpha 每次迭代会下降，但是因为加了常数，所以不会到0
            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]) #  [14.62521125  0.98612686 -1.90030649]
    return weights

这里有2点改进:

$\alpha$ 在每次迭代时都会调整，缓解了上图数据波动会高频波动。

保证 $\alpha \neq0$ 是为了保证多次迭代后新数据仍然具有一定的影响。

如果要处理的问题是动态变化的，那么可适当加大上述常数项，保证新的值获得更大的回归系数。

在降低 $\alpha$ 的函数中， $\alpha$ 每次减少 $\frac1{1+i+j}$ ，其中j 是迭代次数，i是样本点的下标。当 $\max(i)$ 时， $\alpha$ 就不是严格下降的。

通过随机选取样本更新回归系数：这种方法将减少周期性波动。每次随机从列表中选出1个值，然后从列表中删除该值。

随着迭代进行，各系数的变化情况如下：

最终拟合的回归直线为：

示例：与疝气病症预测病马的死亡率

该数据集包含368个样本和28个特征，数据集中有30%的值是缺失的。

数据集内容如下：

2.000000	1.000000	38.500000	66.000000	28.000000	3.000000	3.000000	0.000000	2.000000	5.000000	4.000000	4.000000	0.000000	0.000000	0.000000	3.000000	5.000000	45.000000	8.400000	0.000000	0.000000	0.000000
1.000000	1.000000	39.200000	88.000000	20.000000	0.000000	0.000000	4.000000	1.000000	3.000000	4.000000	2.000000	0.000000	0.000000	0.000000	4.000000	2.000000	50.000000	85.000000	2.000000	2.000000	0.000000
2.000000	1.000000	38.300000	40.000000	24.000000	1.000000	1.000000	3.000000	1.000000	3.000000	3.000000	1.000000	0.000000	0.000000	0.000000	1.000000	1.000000	33.000000	6.700000	0.000000	0.000000	1.000000
1.000000	9.000000	39.100000	164.000000	84.000000	4.000000	1.000000	6.000000	2.000000	2.000000	4.000000	4.000000	1.000000	2.000000	5.000000	3.000000	0.000000	48.000000	7.200000	3.000000	5.300000	0.000000
2.000000	1.000000	37.300000	104.000000	35.000000	0.000000	0.000000	6.000000	2.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	74.000000	7.400000	0.000000	0.000000	0.000000
2.000000	1.000000	0.000000	0.000000	0.000000	2.000000	1.000000	3.000000	1.000000	2.000000	3.000000	2.000000	2.000000	1.000000	0.000000	3.000000	3.000000	0.000000	0.000000	0.000000	0.000000	1.000000
1.000000	1.000000	37.900000	48.000000	16.000000	1.000000	1.000000	1.000000	1.000000	3.000000	3.000000	3.000000	1.000000	1.000000	0.000000	3.000000	5.000000	37.000000	7.000000	0.000000	0.000000	1.000000
1.000000	1.000000	0.000000	60.000000	0.000000	3.000000	0.000000	0.000000	1.000000	0.000000	4.000000	2.000000	2.000000	1.000000	0.000000	3.000000	4.000000	44.000000	8.300000	0.000000	0.000000	0.000000
2.000000	1.000000	0.000000	80.000000	36.000000	3.000000	4.000000	3.000000	1.000000	4.000000	4.000000	4.000000	2.000000	1.000000	0.000000	3.000000	5.000000	38.000000	6.200000	0.000000	0.000000	0.000000
2.000000	9.000000	38.300000	90.000000	0.000000	1.000000	0.000000	1.000000	1.000000	5.000000	3.000000	1.000000	2.000000	1.000000	0.000000	3.000000	0.000000	40.000000	6.200000	1.000000	2.200000	1.000000

首要问题是解决缺失的数据。

处理数据中的缺失值

可选做法：

使用可用特征的均值填补缺失值；
使用特殊值填补缺失值，如-1；
忽略有缺失值的样本；
使用相似样本的均值添补缺失值；
使用另外的机器学习算法预测缺失值。

测试算法

#  以回归系数和特征向量作为输入来计算对应的 Sigmoid值
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])) # 标签
    # 迭代1000次求回归系数
    trainWeights = stocGradAscent1(array(trainingSet), trainingLabels, 1000)
    # 读取测试数据集
    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


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)))

classifyVector 以回归系数和特征向量作为输入来计算对应的 Sigmoid值。如果 Sigmoid 值大于0.5函数返回1，否则返回0。
colicTest() 具有完全独立的功能，多次运行得到的结果可能稍有不同，这是因为其中有随机的成分在里面。
multiTest() 调用 colicTest() 函数10次并求取结果的平均值。

实际运行结果：

the error rate of this test is: 0.313433
the error rate of this test is: 0.447761
the error rate of this test is: 0.388060
the error rate of this test is: 0.417910
the error rate of this test is: 0.238806
the error rate of this test is: 0.283582
the error rate of this test is: 0.402985
the error rate of this test is: 0.417910
the error rate of this test is: 0.417910
the error rate of this test is: 0.417910
after 10 iterations the average error rate is: 0.374627