机器学习实战四：Logisitic回归

前言

这一章开始就学习了最优化算法，由于机器学习实战对于代码比较侧重，所以对于数学原理需要一定的基础。之前没有基础，到网上去找了一些算法中涉及到的数学原理进行学习，大概理解了整个流程，对于一些具体的数学推导还是存在着疑问，所以后边不单单是学习这本书，还要找一些数学原理比较多的书籍进行理解。

理论

回归是一种极易理解的模型，相当于y=f(x)，通过自变量x也就是我们输入的数据点，来获取因变量y，即预测分类。最简单的是线性回归，逻辑回归模型实在线性回归的基础上，套用了一个逻辑函数，这个逻辑函数就是

1

上面这个函数为一个单位阶跃函数，由于逻辑回归多用于二分类，所以我们使用这个函数可以很好的将数据点分为0或者是1，来进行分类。也就是说对于待分类的数据，我们将数据进行处理，然后带入这个逻辑函数后，只会返回两种结果，结果为0或者结果为1，我们就基于这两种结果对数据进行分类。而在代入这个函数之前，还需要对数据进行处理，这个处理过程就是我们的回归过程。

逻辑回归和线性回归在原理上是一致的，在回归时，整个步骤如下：

（1）找到一个合适的预测函数，一般表示为h函数即hypothesis，该函数是我们需要找的分类函数，他的输入就是我们的待分类数据点，输出就是我们对其预测的结果。这个过程非常关键我们需要对数据有一定得了解或分析，知道或者猜测函数的大概形式，比如线性函数还是非线性函数。

逻辑回归的h函数有两部分组成，一个是我们上边说的逻辑函数，另一个即我们进行数据处理的过程；线性边界函数，线性边界的形式如下

2

所以综上所述，我们整个h函数就是

3

我们要让输入的数据值经过我们的h函数就分为0或者1。那么后边的步骤我们就需要找到最佳的回归系数。

（2）可以说我们现在的工作就是要找到一条最佳的直线来对分类结果不同的点进行最佳的划分，这个过程需要不断地修正，不断地与结果进行对比。构造一个Cost函数（损失函数），该函数表示预测的输出（h）与训练数据类别（y）之间的偏差，可以是二者之间的差（h-y）或者是其他的形式。综合考虑所有训练数据的“损失”，将Cost求和或者求平均，记为J(θ)函数，表示所有训练数据预测值与实际类别的偏差。当J(θ)最小时，就表示数据预测值与实际类别的偏差值就最小，使用梯度下降算法求J(θ)的最小值。J的值就是H-y，首先我们来看梯度下降算法公式：

5

这里w:=表示值是覆盖原来的w，而不是等于，相当于一个迭代的过程。我们这里构造J(θ)：

7

使用梯度下降算法来求得J(θ)的最小值

8

所以最后梯度下降算法的公式经过推导过后为

9

w为要找的最佳线性系数，α为步长，通过不断地迭代训练数据使J(θ)走到最小值。

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
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])])
        labelMat.append(int(lineArr[2]))
    return dataMat,labelMat
def sigmoid(inX):
    return 1.0/(1+exp(-inX))
def gradAscent(dataMatIn, classLabels):
    dataMatrix= mat(dataMatIn)
    labelMat=mat(classLabels).transpose()
    m,n=shape(dataMatIn)
    alpha=0.001
    maxCycles=500
    weights=ones((n,1))
    for k in range(maxCycles):
        h=sigmoid(dataMatrix*weights)
        error=(labelMat-h)
        weights=weights+alpha*dataMatrix.transpose()*error
    return weights

loadDataSet，主要功能是打开testSet.txt，并逐行读取。每行前两个值为X1和X2，第三个值为对应的类别标签4

gradAscent就是梯度上升算法，该函数有两个参数：dataMatIn和classLabels，transpose为矩阵转置函数，将1×100的行向量变为100×1的列向量。weights权重也就是上边公式中的w初始化为1，maxCycles是迭代次数初始化为500，变量h是一个列向量，sigmoid函数实际上输入的是一个矩阵的乘法，即所有训练数据与权重weights的乘法，然后将这个数据和训练数据的分类结果进行比较，找到差值，这个差值就是我们的J(θ)即h-y的值。然后对这个w向量进行修正，也就是我们的梯度上升算法。

整个过程下来，迭代了500次，为了找到w最佳参数，每次都需要将训练集的所有数据和w向量相乘，然后使用sigmoid函数输出0或1表示训练数据的分类。由于需要不断迭代来找到最佳的w向量，将输出h和真正的分类数据label相减，然后使用梯度上升算法。修正出一组新的weights，继续进行迭代。

下面我们通过代码画出决策边界，更为直观地确定我们的w参数是否正确

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
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]
    ax.plot(x, y)
    plt.xlabel('X1'); plt.ylabel('X2');
    plt.show()

figure_1

这条最佳拟合直线虽然效果很好但是却需要大量计算，主要计算量大的一步为h=sigmoid(dataMatrix*weights)

由于是矩阵乘法，所以其实每次运算都需要三百次乘法，这个计算量太消耗时间，所以后边对其进行改进，随机梯度上升算法：

1
2
3
4
5
6
7
8
9
10
11
12
13
def stocGradAscent(dataMatrix, classLabels, numIter=150):
    m,n=shape(dataMatrix)
    weights=ones(n)
    for j in range(numIter):
        dataIndex=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

这里对h进行，h不是一个向量而是一个数值，也就是说我们随机抽取一组训练集数据让他对我们的w向量进行修正，而不是对于整个矩阵进行h运算然后与label的列向量进行求差然后进行修正。并且对alpha进行改进，alpha在每次迭代时都会调整，这会缓解数据波动或者高频波动，因为存在常数项，所以尽管alpha每次都会减小但不会减少为0。最终效果：

figure_2

实践

使用Logistic回归来预测患有疝病的马的存活问题。这里的数据包括368个样本和28个特征。疝病是描述马肠胃痛的术语。然而，这种病不一定源于马的胃肠问题，有的指标比较主观，有的指标难以测量，例如马的疼痛级别。数据还存在一个问题，数据集中有30%的值是缺失的。

对于缺失值直接补为0，然后保存为两个文件：horseColicTest.txt，horseColicTraining.txt。然后利用分类器来对马生死问题进行预测，也就是每匹马都患有疝病，我们通过测量的马的数据来预测怎样特征的马在患有疝病后能够存活。

把测试集的每个特征向量乘以最优化得来的回归系数，再将该乘积结果求和，最后输入到Sigmoid函数中。如果对应的Sigmoid值大于0.5就预测类别标签为1，否则为0.

1
2
3
4
def classifyVector(inX, weights):
    prob= sigmoid(sum(inX*weights))
    if prob>0.5: return 1.0
    else: return 0.0

classifyVector函数，它以回归系数和特征向量作为输入来计算对应的Sigmoid值，如果Sigmoid值大于0.5函数返回1，否则返回0.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
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]))
    trainingweights=stocGradAscent(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),trainingweights))!=int(currLine[21]):
            errorCount+=1
    errorRate = (float(errorCount)/numTestVec)
    print "the error rate of this test is: %f" %errorRate
    return errorRate

colicTest函数用于打开测试集和训练集，并对数据进行格式化处理的函数。该函数首先导入训练集，同前面一样，数据的最后一列仍然是类别标签。数据导入之后使用函数stocGradAscent来计算回归系数向量。这里进行了一千次迭代。在系数计算完成后，导入测试集并计算分类错误率。因为有随机的成分在里边，所以多次运行后结果会稍有不同。

1
2
3
4
5
def multiTest():
    numTest=10; errorSum=0.0
    for k in range(numTest):
        errorSum+=colicTest()
    print "after %d iterations the average error rate is: %f"  %(numTest,errorSum/float(numTest))

multiTest函数调用colicTest函数10次，并求结果的平均值，由于存在30%的数据损失，所以最后输出结果会有37%，但是效果已经非常好了。

6

Comments

⬆︎TOP