初探logistic回归

优点：计算代价低，易于实现。
》缺点：容易欠拟合，分类精度可能不高，只能处理二分类问题且必须线性可分。

<!- more-->

1.1 何为logistic回归

逻辑回归的用途十分广泛，主要用来进行预测，也被看做是一种概率估计。这一点可以参考吴军老师的《数学之美》——第28章　逻辑回归和搜索广告，还有刘鹏老师的《计算广告》。

梯度上升法

1.2 构造分类器

为什么不使用海维赛德阶跃函数（单位阶跃函数）：
在跳跃点上瞬间从0跳跃到1，这个过程太激凸...
我们用sigmoid函数，温和不刺激~

SigmoidFunction

这样我们就有了基本的构件模型，下一步就是函数的输入z
我们将特征变量乘以回归系数再全部累加得到Z:

z=w0x0 + w1x1 + w2x2 .... + wnxn

其输出分大于0.5和小于0.5，表示两个类别，也就实现了分类，确定了分类器的函数形式，接下来问题就是求最佳回归系数。（梯度上升法，随机梯度上升法）
主要思想：要找到某函数的最大值，最好的办法是沿着该函数的梯度方向探寻。

梯度是矢量，总是沿着函数值上升最快的方向移动，因此我们沿着梯度方向或者反方向行进时，就能达到一个函数的最大值（梯度上升法）或者最小值（梯度下降法），因此梯度上升算法就是不断更新梯度值，直到梯度不再变化或者变化很小，即函数达到了最大值
梯度算法的迭代公式为（alpha为步长，即每一步移动量）：

而我们常说的梯度下降法就是将其中的加号换为减号。

梯度上升法用来寻找最大值，而梯度下降法用来寻找最小值。

from numpy import *
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
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
    maxCycles = 500
    weights = ones((n, 1))
    for k in range(maxCycles):#注意range用法，此处是从0到maxCycles而不含maxCycles
        h = sigmoid(dataMatrix*weights)
        error = (labelMat - h)
        weights = weights + alpha * dataMatrix.transpose()*error
    return weights


def plotBestFit(weights):
    dataMat, labelMat = loadDataSet()
    dataArr = array(dataMat)
    n = shape(dataArr)[0]
    xcode1 = []; ycode1 = []
    xcode2 = []; ycode2 = []
    for  i in range(n):
        if int(labelMat[i]) == 1:
            xcode1.append(dataArr[i, 1]); ycode1.append(dataArr[i, 2])
        else:
            xcode2.append(dataArr[i, 1]); ycode2.append(dataArr[i, 2])
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.scatter(xcode1, ycode1, s = 30, c = 'red', marker = 's')
    ax.scatter(xcode2, ycode2, 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()


def main():
    dataArr, labelMat =loadDataSet()
    weights = gradAscent(dataArr, labelMat)
    plotBestFit(weights.getA())#getA()将会返回矩阵的array版本



if __name__ == "__main__":
    main()

可以看到划分的结果：

梯度上升法划分结果

那么，该如何提升呢？
由于计算最佳拟合曲线的系数时，梯度上升法都要遍历整个数据集，当数据量十分巨大时，很难想象性能是多么差，于是我们引入一种在线学习方法，说白了就是增量式的更新系数。
我们叫随机梯度上升法。

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

结果如下图所示：

随机梯度上升法

我们会发现，数据的分类比较糟糕，并没有达到理想的效果，原因是它的迭代次数远没有梯度上升法的500次多，需要对其进一步修改。

def stocGradAscent1(dataMatrix, classLabels, numIter=150):
    m,n = shape(dataMatrix)
    weights = ones(n)   #initialize to all ones
    for j in range(numIter):
        dataIndex = range(m)
        for i in range(m):
            alpha = 4/(1.0+j+i)+0.0001    #apha decreases with iteration, does not 
            randIndex = int(random.uniform(0,len(dataIndex)))#go to 0 because of the constant
            h = sigmoid(sum(dataMatrix[randIndex]*weights))
            error = classLabels[randIndex] - h
            weights = weights + alpha * error * dataMatrix[randIndex]
            del(dataIndex[randIndex])
    return weights

经过修改之后再运行可得结果如下图所示：

改进的随机梯度上升法

最后，我们想要预测一下病马的死亡率，根据给出的数据文件horseColicTraining.txt和horseColicTest.txt，在看过数据之后会发现这是一个脏数据，需要对数据进行预处理即数据清洗，将缺失数据项的数据进行补充，缺失标签的进行剔除。

def classifyVector(inX, weights):#计算sigmoid的值，输入为特征向量和回归系数
    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, 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(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():#调用colicTest()10次并求平均值
    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))

def main():
    # dataArr, labelMat =loadDataSet()
    # weights = stocGradAscent1(array(dataArr), labelMat)
    # plotBestFit(weights)
    multiTest()
if __name__ == "__main__":
    main()

---

逻辑回归的目的在于寻找一开始我们找到的sigmoid阶跃函数的最佳拟合系数，而寻找的方法使用了最优化方法（梯度上升法/梯度下降法）。逻辑回归主要适用于对数据进行分类，而后在分类的基础上进行进一步分析。
（完）