logistic（分享自其他博主）

作者：小村长  出处：http://blog.csdn.net/lu597203933 欢迎转载或分享，但请务必声明文章出处。 （新浪微博：小村长zack, 欢迎交流！）

1：简单概念描述

假设现在有一些数据点，我们用一条直线对这些点进行拟合(该线称为最佳拟合直线)，这个拟合过程就称为回归。训练分类器就是为了寻找最佳拟合参数，使用的是最优化算法。

这就是简单的线性回归问题，可以通过最小二乘法求解其参数，最小二乘法和最大似然估计见：http://blog.csdn.net/lu597203933/article/details/45032607。 但是当有一类情况如判断邮件是否为垃圾邮件或者判断患者癌细胞为恶性的还是良性的，这就属于分类问题了，是线性回归所无法解决的。这里以线性回归为基础，讲解logistic回归用于解决此类分类问题。

基于sigmoid函数分类：logistic回归想要的函数能够接受所有的输入然后预测出类别。这个函数就是sigmoid函数，它也像一个阶跃函数。其公式如下：

其中： z = w₀x₀+w₁x₁+….+w_nx_n,w为参数， x为特征

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

亦即公式表示为：

g(z)曲线为：

此时就可以对标签y进行分类了：

其中θ^Tx=0 即θ0+θ1*x1+θ2*x2=0 称为决策边界即boundarydecision。

Cost function:

线性回归的cost function依据最小二乘法是最小化观察值和估计值的差平方和。即：

但是对于logistic回归，我们的cost fucntion不能最小化观察值和估计值的差平法和，因为这样我们会发现J(θ)为非凸函数，此时就存在很多局部极值点，就无法用梯度迭代得到最终的参数(来源于AndrewNg video)。因此我们这里重新定义一种cost function

通过以上两个函数的函数曲线，我们会发现当y=1，而估计值h=1或者当y=0，而估计值h=0,即预测准确了，此时的cost就为0，，但是当预测错误了cost就会无穷大，很明显满足cost function的定义。

可以将上面的分组函数写在一起：

这样得到总体的损失函数J(θ)为：

       梯度上升法：基于的思想是要找到某函数的最大值，最好的方法是沿着该函数的梯度方向探寻。

该公式将一直被迭代执行，直到达到某个停止条件为止，比如迭代次数达到某个指定值或者算法达到某个可以允许的误差范围。

这样我们依据上面的J(θ)就可以得到梯度上升的公式：

当然上图中少了个求和符号。这样就得到

当然对于随机化的梯度迭代每次只使用一个样本进行参数更新，就为：

这也是下面代码中公式的来源。

例如：data=[1,2,3;4,5,6;7,8,9;10,11,12]为4个样本点，3个特征的数据集，，此时标签为[1,0,0,0]，

那么用梯度上升  表达的就是当j=0时，就是第一列[1,4,7,10]与标签差的乘积。。做解释如下：

for k in range(maxCycles):
        h = sigmoid(dataMatrix*weights)
        error = (labelMat - h)
        weights = weights + alpha * dataMatrix.transpose() * error //此段代码为《机器学习实战代码》，主要用来进行梯度上升过程

       

labelMat         dataMatrix                         error

如上图所示，labelMat表示数据集矩阵（4*3），说明共有四条样例，3列表明有三个特征，也就是说有三个系数需要进行学习。从而可知weight向量为（3*1）形式。假设当前利用上图中的dataMatrix与当前的weight向量计算后得到的值带入sigmoid函数后，得到了一个h向量（该h向量为利用当前系数计算得出的分类值，由于讨论0，1分类情况，所以该h向量的值应该是在0，1附近的值。h向量很显然是4*1的，与labelMat格式一致）。在得到h向量后，与真实的labelMat做想减操作，得到偏差值，也就是入股error向量。下面进入到代码最重要的部分，第四行代码。将dataMatrix转置。如下图：

 

转置后得到3*4的矩阵，将该矩阵与4*1格式的error相乘，得到3*1的向量，该向量可以理解为就是梯度下降的方向，再乘以alpha（步长），再与当前的weight向量做相加操作，得到新的weight值。（也就是学习出来的新的系数向量。）也就是博文作者总结的：在一次迭代过程中，对于系数向量中的第j个系数的变化量，就是由未转置的datamatrix的第j列与标签差的乘积和。再回头看给出的这个式子:发现就是刚才提及的意思，计算第j个系数时，依次让第i个样例的第j个系数乘以真实标签与第i个样例计算得出标签的差值做乘积，再sum以下。（这个式子里没有向量的运算在其中）

 

为什么要采用上面的函数作为cost function？

Andrew Ng给的解释是因为最小估计值和观察值的差平方和为非凸函数，通过函数曲线观察得到上面的cost function满足条件。

这里给出另外一种解释——最大似然估计

我们知道h_θ(x)≥0.5<后面简用h>,此时y=1， 小于0.5，y=0. 那么我们就用h作为y=1发生的概率，那么当y=0时，h＜0.5，此时不能用h作为y=0的概率，<因为最大似然的思想使已有的数据发生的概率最 大化，小于0.5太小了>，我们可以用1-h作为y=0的概率，这样就可以作为y=0的概率了，，然后只需要最大化联合概率密度函数就可以了。

这样联合概率密度函数就可以写成：

再转换成对数似然函数，就和上面给出的似然函数一致了。

代码如下：

****************************梯度上升***********************************************

def gradAscent(dataMatIn,classLabels):
    dataMatrix=mat(dataMatIn)#利用numpy带有的矩阵，使运算简单
    labelMat=mat(classLabels).transpose()#labelMat转置称为列向量
    m,n=shape(dataMatrix)
    alpha=0.001
    maxCycles=500
    weights=ones((n,1))#weights的长度与特征长度相同
    for k in range(maxCycles):
        h=sigmoid(dataMatrix*weights)#利用当前weight与数据集计算出值，代入sigmoid函数中，计算得出预测的分类
        error=labelMat-h#计算预测分类与真实分类的差值(如果写成h-labelMat，在计算sigmoid时发生overflow)
        weights=weights+alpha*dataMatrix.transpose()*error#weight中的第i个值（系数），就是数据集中的第i列（共有n个数）与偏差值（共有n个数）之对应乘积之和
    return weights#得出最终的系数向量

***************************随机梯度上升***********************************************

def stroGradAscent0(dataMatrix,classLabels):#dataMatrix是array
    m,n=shape(dataMatrix)
    alpha=0.01
    weights=ones(n)#行向量,array
    for i in range(m):
        h=sigmoid(sum(dataMatrix[i]*weights))#array之间的计算,matrix[i]仍然是matrix，所以dataMatrix需要时array类型
        error=classLabels[i]-h
        weights=weights+alpha*error*dataMatrix[i]#每次都只是用一个样本进行更新（为什么不用转置，参看上文中的推导公式）
    return weights

***************************随机梯度上升的优化**************************

def stocGradAscent1(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.0001    #alpha会随着迭代次数的不断增加而不断减小，但是由于有常数项，所以永远不会减到0。保证新数据永远有影响
            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

********************************************************************

现在可以利用已有的函数，进行实际应用:从疝气症预测病马的死亡率：代码很简单，不做解释了

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(' ')
        lineArr =[]
        for i in range(21):
            lineArr.append(float(currLine[i]))
        trainingSet.append(lineArr)
        trainingLabels.append(float(currLine[21]))
    trainWeights = stocGradAscent1(array(trainingSet), trainingLabels, 1000)
    errorCount = 0; numTestVec = 0.0
    for line in frTest.readlines():
        numTestVec += 1.0
        currLine = line.strip().split(' ')
        lineArr =[]
        for i in range(21):
            lineArr.append(float(currLine[i]))
        if int(classifyVector(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))