在介绍logistic回归之前，我们需要先大致了解一下回归的概念，假设我们现在有一些数据点，我们用一条线（直线或者曲线都行）对这些点进行拟合，将这些点的分布大致符合这条线的轨迹，这个拟合的过程就称为回归。而我们今天聊的Logistic回归的主要思想就是：根据现有的数据点，对分类边界线建立回归公式（注意是对分类边界线建立回归公式，既拟合出一条最适合分类数据点的线来）。

1.2logistic回归一般过程

logistic回归一般过程就是：(1)收集数据

（2）准备数据

（3）分析数据

（4）训练算法

（5）测试算法

（6）使用算法

1.3sigmoid函数特点

对于logistic回归，我们需要一个函数，接受所有的输入然后预测出类别。例如，在数据只有两个类别的情况下，上述函数应该要输出0或1。《机器学习实战》一书中提到过海维塞德阶跃函数，又称单位阶跃函数，如下图：

这个函数虽然可以很醒目的将0和1分类，但是在从0到1的过程是瞬间跃迁的，这个过程很难处理，我们需要一种更平滑的曲线来过度从0到1的过程，于是，我们便选择了sigmoid函数：

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

在图像中，当x=0时，sigmoid函数值为0.5，随着x的增大，sigmoid的值逼近于1；而随着x的减小，sigmoid的值逼近于0。将坐标轴无限放大，也类似一种阶跃函数。这就实现了0和1的分类，并且0到1的过程中是平滑的。

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

二、初步了解优化算法

在确定选择sigmoid函数作为分类函数之后，现在需要确定最佳回归系数的大小。

2.1梯度上升法

梯度上升的思想是：要找到某个函数的最大值，最好的方法是沿着该函数的梯度方向探寻，将梯度记为 $\bigtriangledown$ ，则f(x,y)的梯度表示为:

$\bigtriangledown 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)必须要在待计算的点上有定义并且可微。

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

在《机器学习实战》书中附上的资源中，下载数据集testSet：

梯度上升算法代码：

def gradAscent(dataMatIn, classLabels):
    dataMatrix = np.mat(dataMatIn)             
    labelMat = np.mat(classLabels).transpose() 
    m, n = np.shape(dataMatrix)
    alpha = 0.001
    maxCycles = 500
    weights = np.ones((n, 1))
    for k in range(maxCycles):              
        h = sigmoid(dataMatrix*weights)     
        error = (labelMat - h)             
        weights = weights + alpha * dataMatrix.transpose()* error 
    return weights

使用matplotlib画出分界线：

def plotBestFit(weights):
    import matplotlib.pyplot as plt
    dataMat, labelMat = loadDataSet()
    dataArr = np.array(dataMat)
    n = np.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 = np.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()

2.3随机梯度上升

在上文中，梯度上升算法在每次更新回归系数时都需要遍历整个数据集，如果样本数量过多，那么该方法的计算复杂度就太高了，于是我们将其改进为随机梯度上升，即以此用一个样本点来更新回归系数。

代码如下：

def stocGradAscent0(dataMatrix, classLabels):
    m, n = np.shape(dataMatrix)
    alpha = 0.01
    weights = np.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

可以看出，使用随机梯度上升之后的分界线反而不准确了，但是这样判断是不准确的，因为判断优化算法的优劣的可靠方法是看它是否收敛，即参数是否到达了稳定值，之后我们还会提供更加优化的方法来改进这个算法。

2.4改进的随机梯度上升算法

def stocGradAscent1(dataMatrix, classLabels, numIter=150):
    m, n = np.shape(dataMatrix)
    weights = np.ones(n)   
    for j in range(numIter):
        dataIndex = list(range(m))
        for i in range(m):
            alpha = 4/(1.0+j+i)+0.0001   
            randIndex = int(np.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

改进的算法通过随机选取样本来更新回归系数，这种方法会减少周期性的波动，此外算法海增加了一个迭代次数作为第三个参数，默认为150，如果给定，那么算法将按照新的参数值进行迭代。

上述是迭代了500次的结果，可以看出，最终分割线还是挺准确的

三、实例：从疝气病症预测病马的死亡率

3.1准备数据

《机器学习实战》书中提供的资源中，本例的数据是缺失的，所有的缺失值必须要用一个实数值来替换，因为我们使用的Numpy数据类型不允许包含缺失值。这里选择实数0来替换所有的缺失值，恰好能适用于Logistic回归。这样做的直觉是我们需要的是一个在更新时不会影响系数的值。另外，由于sigmoid(0)=5，即它对结果的预测不具有任何倾向性，因此上述做法也不会对误差造成任何影响。
同时，如果测试集中一条数据的类别标签已经缺失，那么我们将该类别数据丢弃，因为类别标签与特征不同，很难确定采用某个合适的值来替换。

3.2测试算法

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

def colicTest():
    frTrain = open('D:/machinelearndata/horseColicTraining.txt'); frTest = open('D:/machinelearndata/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(np.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(np.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)))