logistic回归

logistic回归

一、概念

logistic回归是统计学习中的经典分类方法，属于对数线性模型，所以也被称为对数几率回归。这里要注意，虽然带有回归的字眼，但是该模型是一种分类算法，逻辑斯谛回归是一种线性分类器，针对的是线性可分问题。利用logistic回归进行分类的主要思想是：根据现有的数据对分类边界线建立回归公式，以此进行分类。这里的“回归”一词源于最佳拟合，表示要找到最佳拟合参数集，因此，logistic训练分类器时的做法就是寻找最佳拟合参数，使用的是最优化方法。

二、Logistic回归函数：

事件的几率，是指该事件发生的概率与该事件不发生的概率的比值。如果事件发生的概率是p，那么该事件的几率是p/(1-p)。取该事件发生几率的对数，定义为该事件的对数几率或logit函数：

 

事件发生的概率p的取值范围为[0,1]，对于这样的输入，计算出来的几率只能是非负的，而通过取对数，便可以将输出转换到整个实数范围内，这样，我们就可以将对数几率记为输入特征值的线性表达式 ：

   其中，p(y =1|x)是条件概率分布，表示当输入为x时，实例被分为1类的概率，依据此概率我们能得到事件发生的对数几率。但是，我们的初衷是做分类器，简单点说就是通过输入特征来判定该实例属于哪一类别或者属于某一类别的概率。所以我们取logit函数的反函数，令w^T*x的线性组合为输入，p为输出，经如下推导：

 公式1就是logistic函数。大家应该对Φ(x)很熟悉，是一个sigmoid函数，类似于阶跃函数的S型生长曲线。

 

上图给出了sigmoid函数在不同坐标尺度下的两条曲线图。当x为0时，sigmoid函数值为0.5。随着x的增大，对应的sigmoid函数的值将逼近于1；而随着x的减小，sigmoid函数的值将逼近于0。而第二幅图中我们能看到在横坐标的刻度足够大是，在x=0处sigmoid函数看起来很像阶跃函数。

     那么对于公式1，我们可以这样解释：为了实现logistic回归分类器，我们可以在每个特征上都乘以一个回归系数，然后把所有的结果值相加，将这个总和带入sigmoid函数中。进而得到一个范围在0-1之间的数值。最后设定一个阈值，在大于阈值时判定为1，否则判定为0。以上便是逻辑斯谛回归算法是思想，公式就是分类器的函数形式。

三、最佳回归系数的确定

我们已经确定了logisti分类函数，有了分类函数，我们输入特征向量就可以得出实例属于某个类别的概率。但这里有个问题，权重w（回归系数）我们是不确定的。正如我们想的那样，我们需要求得最佳的回归系数，从而使得分类器尽可能的精确。

       如何才能获得最佳的回归系数呢？这里就要用到最优化方法。在很多分类器中，都会将预测值与实际值的误差的平方和作为损失函数（代价函数），通过梯度下降算法求得函数的最小值来确定最佳系数。前面我们提到过某件事情发生的概率为p，在logistic回归中所定义的损失函数就是定义一个似然函数做概率的连乘，数值越大越好，也就是某个样本属于其真实标记样本的概率越大越好。如，一个样本的特征x所对应的标记为1,通过逻辑斯蒂回归模型之后，会给出该样本的标记为1和为-1的概率分别是多少，我们当然希望模型给出该样本属于1的概率越大越好。既然是求最大值，那我们用到的最优化算法就是梯度上升，其实也就是与梯度下降相反而已。

1、我们需要先定义一个最大似然函数L，假定数据集中的每个样本都是独立的，其计算公式如下：

 L(w)就是我们以上说的对于损失函数的最原始的定义，但我们还要进一步处理，那就是取对数，首先，在似然函数值非常小的时候，可能出现数值溢出的情况，使用对数函数降低了这种情况发生的可能性。其次，我们可以将各因子的连乘转换为和的形式，利用微积分中的方法，通过加法转换技巧可以更容易地对函数求导。

2、取似然函数的对数：

 这里似然函数是取最大值的，我们可以直接将其确定为损失函数然后使用梯度上升算法求最优的回归系数。但是既然是损失函数，应该还是取最小值好一点，且我们最耳熟能详的方法也是梯度下降，很多书中讲到这里也都是用的梯度下降算法，所以我将以上公式取反来使用梯度下降算法来求最小值。

 为了更好地理解这一最终损失函数，先了解一下如何计算单个样本实例的成本：

 通过上述公式发现：如果y=0，则第一项为0；如果y=1，则第二项等于0；

下图解释了在不同Φ(z)值时对单一示例样本进行分类的代价：

 

 

从上图中可以看到，若y=1，在Φ(z)值越大的情况下，红色的曲线越趋近于0。Φ(z)越大，那么样本被正确划分到1类的概率就越大，这时损失就越小。相对的，Φ(z)越小，意味着样本越有可能属于0类，但是我们知道y=1，在样本属于1类，所以在分类错误的情况下，代价将趋近于无穷。绿线代表着y=0是情况，当然是跟红线相反，接下来我们就要使用梯度下降求最小值了。首先，计算对数似然函数对j个权重的偏导：

 在进入下一步之前，先计算一下sigmoid函数的偏导：

 将（2）式带入（1）式中，

 我们的目标是求得使损失函数最小化的权重w，所以按梯度下降的方向不断的更新权重：

 由于我们是同时更新所有的权重的，因此可以将更新规则记为：

 其中，Δ w 定义为：

 所以综上，我们将梯度下降的更新规则定义为：

 四、梯度上升算法

梯度上升法基于的思想是：要找到某函数的最大值，最好的方法是沿着该函数的梯度方向探寻。如果梯度记为∇，则函数f(x,y)的梯度由下式表示：

 

梯度上升算法到达每个点后都会重新估计移动的方向。从P0开始，计算完该点的梯度，函数就根据梯度移动到下一点P1。在P1点，梯度再次被重新计算，并沿新的梯度方向移动到P2.如此循环迭代，知道满足通知条件。迭代的过程中，梯度算子总是保证我们能选取到最佳的移动方向

       上图中的梯度上升算法沿梯度方向移动了一步。可以看到，梯度算子总是指向函数值增长最快的方向。这里所说的是移动方向，而未提到移动量的大小。该量值称为步长，记作α \alphaα。用向量来表示的话，梯度上升算法的迭代公式如下：

 

 

 五、 测试程序

from numpy import *


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)  # convert to NumPy matrix
    labelMat = mat(classLabels).transpose()  # convert to NumPy matrix
    m, n = shape(dataMatrix)
    alpha = 0.001
    maxCycles = 500
    weights = ones((n, 1))
    for k in range(maxCycles):  # heavy on matrix operations
        h = sigmoid(dataMatrix * weights)  # matrix mult
        error = (labelMat - h)  # vector subtraction
        weights = weights + alpha * dataMatrix.transpose() * error  # matrix mult
    return weights


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()


if __name__ =='__main__':
    dataArr,labelMat=loadDataSet()
    weights=gradAscent(dataArr,labelMat)
    plotBestFit(weights.getA())

运行结果：