机器学习实战-logistics回归

一.引言

假设我们现有一些数据点，我们用一条直线对这些点进行拟合，这个拟合的过程就称作回归。利用logistic回归进行分类的主要思想是：根据现有数据对分类边界线建立回归公式，以此进行分类。

     sigmoid=1/(1+exp(-inX)),该函数的曲线类似于一个s型，在x=0处，函数值为0.5.

　　于是，为了实现logistic分类器，我们可以在每个特征上都乘以一个回归系数，然后所有的相乘结果进行累加，将这个总结作为输入，输入到sigmoid函数中，从而得到一个大小为0~1之间的值，当该值大于0.5归类为1，否则归类为0，这样就完成了二分类的任务。所以logistic回归可以看成是一种概率估计。

二.基于最优化方法的最佳回归系数确定

　　每个特征乘一个回归因数，然后把所有结果值相加，放入sigmoid函数，进而得到0-1间数值，任何大于0.5分为1类，反之0类。先采用梯度上升法。

2.1 梯度上升法

　　梯度上升法(等同于我们熟知的梯度下降法，前者是寻找最大值，后者寻找最小值)，它的基本思想是：要找到某函数的最大值，最好的方法就是沿着该函数的梯度方向搜寻。如果函数为f，梯度记为D，a为步长，那么梯度上升法的迭代公式为：w：w+a*Dwf(w)。该公式停止的条件是迭代次数达到某个指定值或者算法达到某个允许的误差范围。我们熟知的梯度下降法迭代公式为：w：w-a*Dwf(w)

使用梯度上升法寻找最佳参数

　　假设有100个样本点，每个样本有两个特征：x1和x2.此外为方便考虑，我们额外添加一个x0=1，将线性函数z=wTx+b转为z=wTx(此时向量w和x的维度均价1).那么梯度上升法的伪代码如下：

　　初始化每个回归系数为1

　　重复R次：

　　　　计算整个数据集梯度

　　　　使用alpha*gradient更新回归系数的向量

　　返回回归系数

代码：

def sigmoid(inX):
    return 1.0/(1+exp(-inX))
#这样写有时候出现下面警告RuntimeWarning: overflow encountered in exp
改成下面的就好了

def sigmoid(inX):
    if inX >= 0:
        return 1.0 / (1 + exp(-inX))
    else:
        return exp(inX) / (1 + exp(inX))

#根据梯度上升法，求函数最大值。梯度下降法求函数最小值
#w = w+af(w)(a是步长，f(w)是偏导)
def gradAscent(dataMatIn,classLabels):
    dataMatrix = mat(dataMatIn)
    labelMat = mat(classLabels).transpose()
    m,n = shape(dataMatrix)
    alpha = 0.001
    maxCycles = 500
    weights = ones((n,1))
    #print(type(weights))
    for k in range(maxCycles):
        h = sigmoid(dataMatrix * weights) #此处开始是matrix和array乘积后来是两个matrix乘积。是矩阵相乘（100×3）*（3×1）=100×1
        error = labelMat-h
        weights = weights + alpha*dataMatrix.transpose()*error #这个地方是推导出来的
    return weights

其中weights = weights + alpha*dataMatrix.transpose()*error推导过程如下：

           写的很好 https://blog.csdn.net/cxjoker/article/details/83002197

2.2随机梯度上升

梯度上升法每次更新回归系数都需要遍历整个数据集，当样本数量较小时，该方法尚可，但是当样本数据集非常大且特征非常多时，那么梯度上升法的计算复杂度就会特别高。一种改进的方法是一次仅用一个样本点来更新回归系数，即随机梯度上升法。由于可以在新样本到来时对分类器进行增量式更新，因此随机梯度上升法是一个在线学习算法。随机梯度上升法可以写成如下伪代码：

所有回归系数初始化为1

对数据集每个样本

　　计算该样本的梯度

　　使用alpha*gradient更新回归系数值

 返回回归系数值

实现代码如下

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
        print(error)
        weights = weights + (alpha * error)* dataMatrix[i]
    return mat(weights)

2.3 随机梯度上升改进算法
 

   我们知道，评判一个优化算法的优劣的可靠方法是看其是否收敛，也就是说参数的值是否达到稳定值。此外，当参数值接近稳定时，仍然可能会出现一些小的周期性的波动。这种情况发生的原因是样本集中存在一些不能正确分类的样本点(数据集并非线性可分)，所以这些点在每次迭代时会引发系数的剧烈改变，造成周期性的波动。显然我们希望算法能够避免来回波动，从而收敛到某个值，并且收敛速度也要足够快。
上面的算法是固定的步长，固定的步长，不稳定，会产生震荡，所以下面的算法采用不固定的步长。
在距离目标值远的时候，步长大，距离目标值近的时候，步长小。

def stocGradAscent1(dataMatrix,classLabels,numIter=150):
    m,n = shape(dataMatrix)
    weights = ones(n)
    for j in range(numIter):
        dataIndex = list(range(m))
        for i in range(m):
            # 每次更新参数时设置动态的步长，且为保证多次迭代后对新数据仍然具有一定影响
            # 添加了固定步长0.1
            alpha = 4/(1.0+j+i)+0.01
            # 随机获取样本
            randIndex = random.randint(0,len(dataIndex))
            h = sigmoid(sum(dataMatrix[randIndex]*weights))
            error = classLabels[randIndex] - h
            weights = weights + (alpha * error)* dataMatrix[randIndex]
            del(dataIndex[randIndex])
    return mat(weights)

1 alpha在每次迭代更新是都会调整，这会缓解数据波动或者高频运动。此外，alpha还有一个常数项，目的是为了保证在多次迭代后仍然对新数据具有一定的影响，如果要处理的问题是动态变化的，可以适当加大该常数项，从而确保新的值获得更大的回归系数。

2 第二个改进的地方是选择随机的样本对参数进行更新，由于增加了随机性，这就防止参数值发生周期性的波动。

3 并且采用上述改进算法后，收敛速度更快

三.画图

#下一步分析数据，画出数据集和logistic回归最佳拟合直线
def plotBestFit(wei):
    import matplotlib.pyplot as plt
    weights = wei.getA() #将自身矩阵变量转化为ndarray类型的变量
    dataMat,labelMat = loadfile()
    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) #生成一个array，起-3.0，终3.0（不包含），步长0.1
    # 设置sigmoid函数为0，拟合曲线为0 = w0*x0+w1*x1+w2*x2, 故x2 = (-w0*x0-w1*x1)/w2, x0为1,x1为x，x2为y，故有：
    y = (-weights[0]-weights[1]*x)/weights[2] #由x求出y
    ax.plot(x,y)
    plt.xlabel('X1');plt.ylabel('X2');
    plt.show()

四.logistic回归实例：从疝气病症预测病马死亡率

from logis1 import *
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).reshape(21,1)
    #print(type(trainWeights))
    errorCount = 0;numTestVec = 0.0
    for line in frTest.readlines():
        numTestVec += 1
        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(errorrate)
    return errorrate
def multiTest():
    numTests = 10;errorSum = 0.0
    for k in range(numTests):
        errorSum += colicTest()
    print(errorSum/10.0)

if __name__=="__main__":
    multiTest()