机器学习实战之Logistic回归

逻辑回归实际上是分类算法，称之为回归是因为要选择最佳参数拟合出分类边界。

1.sigmoid函数

当输入z=0时，输出为1；随着x值增大，函数值逼近1；随着z减小，函数值逼近零。当从较长的尺度来看，sigmoid函数很像一个节约函数。

利用sigmoid函数作为分类器，我们可以将每个特征乘上一个回归系数，然后求和作为输入z，输出大于0.5的归为1类，输出不大于0.5的归为0类。

即令 ,

x为特征向量,为待求回归参数。

2.梯度下降法

梯度下降法是一种求解函数最小值的优化算法。即每次迭代都沿着梯度下降的方向。

其算法为：

                                            

其中为学习率，需要选择恰当。

带入求梯度后，Logistic回归梯度下降法求参数的算法为：

记特征矩阵为m*n的矩阵X，则上式变为

 其中和y为向量。

3.代码

（1）数据导入

from numpy import *
#数据导入，输出为两个list，分别为特征矩阵和标签
def dataLoad(filename):
    dataMat=[]
    labelMat=[]
    fr=open(filename)
    for line in fr.readlines():
        lineArr=line.strip().split()
        dataMat.append([1.0,float(lineArr[0]),float(lineArr[1])])
        labelMat.append(float(lineArr[2]))
    return dataMat,labelMat

（2）梯度下降实现

#定义sigmoid函数
def sigmoid(inX):
    return (1.0/(1.0+exp(-inX)))

def gradAscent(datamatIn,labelclss):
    dataMatrix=mat(datamatIn)
    labelMatrix=mat(labelclss).transpose()
    m,n=shape(dataMatrix)
    alpha=0.001#学习率
    cycles=500#循环次数
    weights=ones((n,1))#参数初始化
    for i in range(cycles):
        error=labelMatrix-sigmoid(dataMatrix*weights)
        weights=weights+alpha*dataMatrix.transpose()*error #梯度下降算法
    return weights.tolist()

（3）可视化

def dataPlot():
    import matplotlib.pyplot as plt
    dataMat,labelMat=dataLoad('D:\python\code\machinelearninginaction\Ch05\\testSet.txt')
    weights=gradAscent(dataMat,labelMat)
    m,n=shape(dataMat)
    x1=[];y1=[]
    x2=[];y2=[]
    for i in range(m):
        if labelMat[i]==1:
            x1.append(dataMat[i][1])
            y1.append(dataMat[i][2])
        if labelMat[i]==0:
            x2.append(dataMat[i][1])
            y2.append(dataMat[i][2])
    fig=plt.figure()
    ax=fig.add_subplot(111)
    ax.scatter(x1,y1,s=30,c='red',marker='s')
    ax.scatter(x2,y2,s=30,c='green')
    x = arange(-3.0, 3.0, 0.1)
    y = (-weights[0][0] - weights[1][0] * x) / weights[2][0]
    ax.plot(x, y)
    plt.xlabel('X1')
    plt.ylabel('X2')
    plt.show()

运行结果如图：

可以看到，分类效果还不错。但是以上算法中，每一次参数的更新都对整个数据集进行了遍历，这种算法称为批处理，当数据集较大时计算量太大，可以采用一种随机梯度法。

（4）随机梯度下降

随机梯度下降每次只采用一个样本点来更新参数

def stoGradAscent0(datamatIn,labelclass):
    m,n=shape(datamatIn)
    alpha=0.01
    weights=ones(n)
    for i in range(m):
        h=sigmoid(sum(datamatIn[i]*weights))
        error=labelclass[i]-h
        weights=weights+[error*alpha*a for a in datamatIn[i]] #随机梯度下降
    return weights

随机梯度下降中每次添加一个样本特征，与之前不同的是，这里的error和alpha都是数值。

运行结果：

效果不如第一次的好，这是因为随机梯度下降运算次数太少，我们可以做些改进

（5）改进随机梯度下降法

def stoGradAscent1(datamatIn,labelclass,numberIter=100):
    m,n=shape(datamatIn)
    weights=ones(n)
    for i in range(numberIter):
        dataIndex=list(range(m))
        for j in range(m):
            alpha=4/(1+i+j)+0.0001 #设置步长随着迭代次数增加而减小
            randIndex=int(random.uniform(0,len(dataIndex)))#随机挑选一个样本
            h=sigmoid(sum(datamatIn[randIndex]*weights))
            error=labelclass[randIndex]-h
            weights=weights+[error*alpha*a for a in datamatIn[randIndex]]
            del(dataIndex[randIndex])#删除已经使用过的样本
    return weights

该程序在原有基础上做了三点优化：

• 1.alpha每次迭代的时候都会调整
• 2.随机选取样本来更新回归系数，这种做法会减少周期性波动
• 3.增加了迭代参数的设置。默认为100次，如果给定，将会按照新的参数值进行迭代

运行结果：

（6）实例 预测疝气病马死亡率

def classifyVector(inX,weights):
    prob=sigmoid(sum(inX*weights))
    if prob>0.5:
        return 1
    else:
        return 0

def colicTest():
    frTrain=open('D:\python\code\machinelearninginaction\Ch05\horseColicTraining.txt')
    frTest=open('D:\python\code\machinelearninginaction\Ch05\horseColicTest.txt')
    trainSet=[]
    trainLabel=[]
    for line in frTrain.readlines():
        currLine=line.strip().split('\t')
        lineArr=[float(each) for each in currLine]
        trainSet.append(lineArr[:21])
        trainLabel.append(lineArr[21])
    trainWeights=stoGradAscent1(trainSet,trainLabel,numberIter=500)
    errorCount=0
    numTest=0
    for line in frTest.readlines():
        numTest+=1
        currLine=line.strip().split('\t')
        lineArr = [float(each) for each in currLine]
        testFeat=lineArr[:21]
        predicLabel=lineArr[21]
        if classifyVector(testFeat,trainWeights)!=predicLabel:
            errorCount+=1
    errorRate=float(errorCount/numTest)
    print("the error rate of this test is: %f" % errorRate)
    return errorRate


def multiTest():
    numTests = 20; errorSum=0.0
    for k in range(numTests):
        errorSum += colicTest()
    print ("after %d iterations the average error rate is: %f" % (numTests, errorSum/float(numTests)))

运行结果：

20次平均预测错误率为35%，因为原始数据有很多缺失，这已经是一个不错的结果了。