机器学习第二周-逻辑回归

机器学习第二周-逻辑回归

本周学习内容

• 分类问题：朴素贝叶斯
• 分类问题：逻辑回归

小插曲：在写第一版的时候，没有及时保存，电脑发生故障死机，导致本版丢失了大部分学习笔记，本版为第二版。所幸大部分代码保存在jupyter，得以保存。以后一定要定时保存啊。

第五章：逻辑回归

Sigmoid函数

Logistic回归虽然叫做回归，但实际上是一种分类算法，其思想是根据现有数据对分类边界线建立回归公式，以此进行分类。二值输出型分类函数，我们之前熟悉的有单位阶跃函数。但是因为单位阶跃函数存在0到1的瞬间跳跃，跳跃处导数为无穷大，因此，在后续分析中不利于我们处理。我们可以选择另外一个值域在0到1之间的函数:Sigmoid函数。图像如下图，x为0时，其函数值为0；x值越大越接近1；x值越小越接近0

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

Sigmoid函数的输入记为z，

梯度上升法

这个梯度意味着要沿 x 的方向移动 ，沿 y 的方向移动 。

这个梯度上升法和之前接触到的梯度下降法是一致的。梯度上升法是找函数局部的最大值，而梯度下降法是用来寻找局部最小值。每一次都向梯度方向移动一定的步长，我们在公式中用α表示步长。

从公式中我们可以看出来，梯度上升法的迭代公式是用的加法，梯度下降法用的是减号。梯度上升法和梯度下降法都有局部最优值的问题，这就说明可能并不能找到全局最优，而是陷入到一个部分最优值中去

Logistic回归原理

每个回归系数初始化为 1
重复 R 次:
    计算整个数据集的梯度
    使用 步长 x 梯度 更新回归系数的向量
返回回归系数

代码实战-简单数据集使用逻辑回归分类

解析数据

from numpy import*
# 导入文件数据
def loadDataSet(file_name):
    # dataMat为原始数据， labelMat为原始数据的标签
    dataMat = []     #原始数据的特征
    labelMat = []     #原始数据的标签
    fr = open(file_name)   #打开要解析的文件路径
    for line in fr.readlines():
        lineArr = line.strip().split()
        # 为了方便计算，我们将 X0 的值设为 1.0 ，也就是在每一行的开头添加一个 1.0 作为 X0
        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.0 / (1 + exp(-inX))

逻辑回归梯度上升算法

def gradAscent(dataMatIn,classLabels):
    dataMatrix=mat(dataMatIn)  #转换为numpy矩阵
    labelMat=mat(classLabels).transpose()   #transpose（）将行向量转换为列向量，转置
    m,n=shape(dataMatrix)   #m是行数（样本数），n是列数（特征数）
    alpha=0.001    #学习率，本书也把他叫做步长
    maxCycles=500   #最大的迭代次数
    weights=ones((n,1))   #生成一个初始全为1的w权重，长度自然等于特征的个数
    for k in range(maxCycles):
        h=sigmoid(dataMatrix*weights)   #这里计算出的h的每一行就是每一个样本通过sigmoid算出的输出值
        error=(labelMat-h)     #error的每一行代表每一个样本输出值和真实值的误差
        weights=weights+alpha*dataMatrix.transpose()*error   #这是一个递归公式经过一系列推到之后得出的公式
    return weights

对于倒数第二行代码，可参考http://blog.csdn.net/achuo/article/details/51160101

数据展示

#数据展示
def plotBestFit(weights):
    '''
    dataArr:样本数据的特征
    labelMat:样本数据的类别标签，即目标变量
    weights:回归系数
    '''
    dataMat,labelMat=loadDataSet('testSet.txt')
    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)
    # w0+w1*x+w2*y=0 => y = (-w0-w1*x)/w2   
    y = (-weights[0]-weights[1]*x)/weights[2]
    ax.plot(x, y)
    plt.xlabel('X'); plt.ylabel('Y')
    plt.show()

分类结果输出显示

其中参数集为

[[ 4.12414349]
 [ 0.48007329]
 [-0.6168482 ]]

随机梯度上升法

对于梯度上升算法，在每一次更新回归系数的时候都需要遍历整个数据集，这样在大的数据集上是需要大量的计算的。我们可以使用一个退而求其次的方法，随机梯度上升法，虽然迭代次数可能会更多，但是计算时间可以明显降低。随机梯度上升法的思想就是每一次仅用一个样本点来更新回归系数。

在线学习:可以增量学习。

批处理：一次处理所有数据。

梯度上升法伪代码

所有回归系数初始化为 1
对数据集中每个样本
    计算该样本的梯度
    使用 alpha x gradient 更新回归系数值
返回回归系数值

python代码

def stocGradAscent0(dataMatrix, classLabels):
    m,n = shape(dataMatrix) #m为样本数，n为特征数
    alpha = 0.01 #步长
    # n*1的矩阵
    # 函数ones创建一个全1的数组
    weights = 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

从上面的例子我们可以看到，在随机梯度上升法中，我们更新回归系数的时候不再将所有样本全部带入到里面算了，而是仅选举其中一个对回归系数进行更新。

评价一个方法好会的因素不能只看结果，我们需要看看其收敛性。也就是是否能达到稳定状态，不会产生来回的波动。

改进版随机梯度上升法

def stocGradAscent1(dataMatrix, classLabels, numIter=150):
    m,n = shape(dataMatrix)
    weights = ones(n)   # 同样的初始化回归系数
    # numIter=150代表迭代次数
    for j in range(numIter):
        dataIndex = list(range(m))
        for i in range(m):
            alpha = 4/(1.0+j+i)+0.0001    # alpha 随着迭代次数的不断增多，而不断减少，但不能为0.其中j是迭代次数， i是样本点的下标
            # random.uniform(x, y) 方法将随机生成下一个实数，它在[x,y]范围内,x是这个范围内的最小值，y是这个范围内的最大值。
            randIndex = int(random.uniform(0,len(dataIndex)))
            # 这个地方类似于第三章，从样本中随机抽一个，进行计算之后再将其删掉
            h = sigmoid(sum(dataMatrix[dataIndex[randIndex]]*weights))
            error = classLabels[dataIndex[randIndex]] - h
            weights = weights + alpha * error * dataMatrix[dataIndex[randIndex]]
            del(dataIndex[randIndex]) #删掉已经使用过的
    return weights

改进之处：随着梯度迭代次数的不断增加，其步长也不断减少。在降低alpha的函数中，alpha每次减少1/(j+i) ，其中j是迭代次数， i是样本点的下标①。这样当j<<max(i)时，alpha就不是严格下降的。

为了减少周期性波动，我们将之前的算法改进为每次随机从样本集中抽取一个，使用完将其删除掉，然后进行多次的迭代。

以下是使用改进版随机梯度上升法获得的结果。

代码实战-从疝气病症预测病马的死亡率

使用Logistic回归来预测患有疝病的马的存活问题。这里的数据①包含368个样本和28 个特征。

开发流程

(1) 收集数据：给定数据文件。  
(2) 准备数据：用Python解析文本文件并填充缺失值。  (3) 分析数据：可视化并观察数据。 
(4) 训练算法：使用优化算法，找到最佳的系数。 
(5) 测试算法：为了量化回归的效果，需要观察错误率。根据错误率决定是否回退到训练 阶段，通过改变迭代的次数和步长等参数来得到更好的回归系数。  
(6) 使用算法：实现一个简单的命令行程序来收集马的症状并输出预测结果并非难事，这 可以做为留给读者的一道习题

在这个数据集上有30%的值是缺失的，通过这个例子，可以学会如何处理数据缺失的问题。

准备数据

我们可以采用以下的一些方法，处理数据缺失值

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

在处理本数据集，我们需要考虑两点：

第一，对于缺失的特征值我们用0表示，原因是，如果使用0值，则该样本并不会影响到回归系数，也就是说不会对回归系数产生更新。另外，sigmoid（0）=0.5，他对预测结果没有任何的倾向性。

第二，如果类标签丢失，我们将该样本丢弃。

Logistic分类函数

书中提供的数据集是已经根据上面的方法处理好了的

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运行结果

the error rate of this test is: 0.253731
the error rate of this test is: 0.328358
the error rate of this test is: 0.313433
the error rate of this test is: 0.298507
the error rate of this test is: 0.298507
the error rate of this test is: 0.328358
the error rate of this test is: 0.283582
the error rate of this test is: 0.343284
the error rate of this test is: 0.447761
the error rate of this test is: 0.283582
after 10 iterations the average error rate is: 0.317910

总结

通过对logistic回归模型学习，我学会了基本的logistic回归模型的思想，梯度上升法，sigmoid函数。通过两个例子，学习了随机梯度上升法，以及改进版的随机梯度上升法。并且了解了处理异常数据的常用方法。