机器学习实战——Logistic回归

Logistic回归，虽然名为回归，但实际是一个分类算法。利用Logistic回归进行分类的主要思想是：根据现有的一些数据对分类边界线建立回归公式，以此进行分类。这一算法的“回归”源于最佳拟合，表示要找到最佳拟合参数集。

Logistic回归的优点是：计算代价不高，易于理解和实现

缺点是：容易欠拟合，分类精度可能不高。

sigmoid函数

为了实现分类的目的，我们需要的拟合函数应该是可以接受所有的输入，同时能够根据输入到输出的映射预测出类别。

海维赛德阶跃函数（又称单位阶跃函数）就是这样的一个函数，但是该函数不连续，在跳跃点直接从0跳转至1，在很多情况下不好处理。因此，这一算法使用另一种函数，sigmoid函数，也称S型生长曲线。这一函数具备许多优秀的性质，可以帮助我们很好地拟合数据。Sigmoid函数的具体计算公式如下：

$\sigma$(z) = $\frac{1}{1+exp(-z)}$

这一函数连续可微，可以将输入映射到0，1之间，求导运算简洁,即：

$\frac{d\sigma (z)}{dz}$ = $\sigma$(z) × (1-$\sigma$(z))

下图是sigmoid函数的图像：

使用最优化方法确定回归系数

如上文所述，记sigmoid函数的输入为z，

z = w₀x₀ + w₁x₁ + w₂x₂ + $\cdots$ + w_nx_n

将系数w(coefficients)和输入均表示成列向量的形式，那么上式可以写成

$z = w^Tx $

最优化方法

梯度上升法

类似于梯度下降法，只是取梯度正方向和梯度负方向存在差异。

梯度上升法的基本思想是：要找到某个函数的最大值，那么最好是沿着梯度方向进行寻找。如果将梯度记为$▽$,则函数f(x,y)的梯度由下式表示：

$▽$f(x,y) = ($\frac{\partial f(x,y)}{\partial x}$,$\frac{\partial f(x,y)}{\partial y}$)

梯度上升法要求所求函数在计算的点上可微，sigmoid函数很好地满足了这一要求。

梯度上升法到达每个点上都会重新估计移动的方向，从起始点开始，计算该点的梯度，沿着该梯度方向进行移动，移动的距离为步长，也称学习率(learning rate)。以$\alpha$表示步长，梯度上升法的迭代公式可以表示为：

w = w + $\alpha$ $▽$_wf(w)

下图是梯度上升法的一个可视化例子：

梯度上升法的实现

def gradAscent(dataMatIn,classLabels):
    dataMatrix = mat(dataMatIn)		#为便于计算，转换为numpy矩阵
    labelMat = mat(classLabels).transpose()		
    m,n = shape(dataMatrix)
    alpha = 0.001		#步长
    maxCycle = 200		#最大迭代次数
    weights = ones((n,1))	#初始化权重系数
    for i in range(maxCycle):
        h = sigmoid(dataMatrix*weights)		
        error = (labelMat - h)
        weights = weights + alpha * dataMatrix.transpose() * error
    return weights

该函数输入的两个参数分别是训练的数据集，和类别标签。为了便于计算，将二者均转换为numpy矩阵，之所以对类别标签进行转置操作，是因为在训练时，我使用了《机器学习实战》一书提供的数据集，导入该数据集后，dataMatIn为一个100×3的矩阵，classLabels为1×100的行向量。为便于矩阵运算，所以进行转置。

在循环体中，系数的迭代公式与上文中介绍的并不一致。原书中省略了一个简单的数学推导。

数学推导如下：

w的迭代公式为：w = w + $\alpha$ $▽$_wf(w)

从而得到上边代码中的迭代公式。

画出决策边界

def plotBestFit(weights):
    #weights = weights.getA()
    weights = array(weights)
    dataMat,labelMat = loadDataSet()
    x1cords = []
    y1cords = []
    x2cords = []
    y2cords = []
    dataArr = array(dataMat)
    n = shape(dataArr)[0]
    for i in range(n):
        if labelMat[i] == 1:
            x1cords.append(dataArr[i,1])
            y1cords.append(dataArr[i,2])
        else:
            x2cords.append(dataArr[i,1])
            y2cords.append(dataArr[i,2])
        
    fig = plt.figure()
    ax = plt.subplot(111)
    ax.scatter(x1cords,y1cords,s=30,color='r',marker='s')
    ax.scatter(x2cords,y2cords,s=30,color='b',marker='o')
    x = arange(-3.0,3.0,0.1)
    y = (-weights[0]-weights[1]*x)/weights[2]
    ax.plot(x,y)
    plt.xlabel('x-axis')
    plt.ylabel('y-axis')
    plt.show()

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

结果如图所示：

随机梯度上升法

下面先给出第一个版本的随机梯度上升法的实现：

def stocGradAscent0(dataMatrix,classLabels):
    m,n = shape(dataMatrix)
    alpha = 0.01
    weights = ones(n)   #返回array类型，1×n大小
    for i in range(m):
        h = sigmoid(sum(dataMatrix[i]*weights))     #list和array相乘作内积之后求和
        error = classLabels[i] - h
        weights = weights + alpha * error * dataMatrix[i]
    return weights

运行该算法三个系数的迭代情况如下：

该算法经过大量迭代才能趋于稳定，另外因为数据集线性不可分等原因，之后仍然会出现小的波动。

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

分类结果：

结果该算法错分了约三分之一的样本

所以对该算法进行改进：对步长进行动态的调整，同时随机选取样本进行系数的更新，来减少周期性波动。

def stocGradAscent1(dataMatrix,classLabels,numIters = 150):
    m,n = shape(dataMatrix)
    weights = ones(n)
    for j in range(numIters):
        #dataIndex = range(m)  #python3range返回range对象
        dataIndex = list(range(m))
        for i in range(m):
            alpha = 4.0/(1+i+j) + 0.01
            randIndex = int(random.uniform(0,len(dataIndex)-1))
            h = sigmoid(sum(dataMatrix[randIndex] * weights))
            error = classLabels[randIndex] - h
            weights = weights + alpha * error * dataMatrix[randIndex]
            del(dataIndex[randIndex])  
    return weights

ataArr,labelMat = loadDataSet()
weights = stocGradAscent1(array(dataArr),labelMat)
plotBestFit(weights)

运行该算法三个系数的迭代情况为：

分类结果为：

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

本次示例使用uci中的数据。由于这一数据集存在30%的缺失值，所以分类效果不会特别理想。

数据集中缺失值的处理

1.直接丢弃该样本

2.重新采集数据

3.使用一个全局常量填充缺失值

4.使用属性的平均值填充缺失值

5.使用与给定元组同一类的所有样本的平均值

6.使用最可能的值填充缺失值：如使用基于Bayesian公式或判定树之类的算法

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.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.0		#计算错误数量
    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()
    errorRate = float(errorSum)/numTests
    print("After %d iteration, the error rate of this test is %f" % (numTests,errorRate))

the error rate of this test is: 0.313433
the error rate of this test is: 0.358209
the error rate of this test is: 0.298507
the error rate of this test is: 0.373134
the error rate of this test is: 0.417910
the error rate of this test is: 0.298507
the error rate of this test is: 0.402985
the error rate of this test is: 0.358209
the error rate of this test is: 0.343284
the error rate of this test is: 0.417910
After 10 iteration, the error rate of this test is 0.358209

完整代码：

from numpy import *
import operator
import matplotlib.pyplot as plt

#导入数据
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

#logistic函数
def sigmoid(inX):
    return 1.0/(1+exp(-inX))


#梯度上升
def gradAscent(dataMatIn,classLabels):
    dataMatrix = mat(dataMatIn)
    labelMat = mat(classLabels).transpose()
    m,n = shape(dataMatrix)
    alpha = 0.001
    maxCycle = 200
    weights = ones((n,1))
    for i in range(maxCycle):
        h = sigmoid(dataMatrix*weights)
        error = (labelMat - h)
        weights = weights + alpha * dataMatrix.transpose() * error
    return weights

'''test:
dataArr,labelMat = loadDataSet()
gradAscent(dataArr,labelMat)
'''

def plotBestFit(weights):
    #weights = weights.getA()
    weights = array(weights)
    dataMat,labelMat = loadDataSet()
    x1cords = []
    y1cords = []
    x2cords = []
    y2cords = []
    dataArr = array(dataMat)
    n = shape(dataArr)[0]
    for i in range(n):
        if labelMat[i] == 1:
            x1cords.append(dataArr[i,1])
            y1cords.append(dataArr[i,2])
        else:
            x2cords.append(dataArr[i,1])
            y2cords.append(dataArr[i,2])
        
    fig = plt.figure()
    ax = plt.subplot(111)
    ax.scatter(x1cords,y1cords,s=30,color='r',marker='s')
    ax.scatter(x2cords,y2cords,s=30,color='b',marker='o')
    x = arange(-3.0,3.0,0.1)
    y = (-weights[0]-weights[1]*x)/weights[2]
    ax.plot(x,y)
    plt.xlabel('x-axis')
    plt.ylabel('y-axis')
    plt.show()
    
'''
dataArr,labelMat = loadDataSet()
weights = gradAscent(dataArr,labelMat)
plotBestFit(weights)
'''
    
def stocGradAscent0(dataMatrix,classLabels):
    m,n = shape(dataMatrix)
    alpha = 0.01
    weights = ones(n)   #返回array类型，1×n大小
    for i in range(m):
        h = sigmoid(sum(dataMatrix[i]*weights))     #list和array相乘作内积
        error = classLabels[i] - h
        weights = weights + alpha * error * dataMatrix[i]
    return weights

def stocGradAscent1(dataMatrix,classLabels,numIters = 150):
    m,n = shape(dataMatrix)
    weights = ones(n)
    for j in range(numIters):
        #dataIndex = range(m)  #python3range返回range对象
        dataIndex = list(range(m))
        for i in range(m):
            alpha = 4.0/(1+i+j) + 0.01
            randIndex = int(random.uniform(0,len(dataIndex)-1))
            h = sigmoid(sum(dataMatrix[randIndex] * weights))
            error = classLabels[randIndex] - h
            weights = weights + alpha * error * dataMatrix[randIndex]
            del(dataIndex[randIndex])  
    return weights

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.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.0
    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()
    errorRate = float(errorSum)/numTests
    print("After %d iteration, the error rate of this test is %f" % (numTests,errorRate))