Logistic回归 from机器学习实战第五章

回归：对于数据点，利用直线对其拟合的过程

Sigmoid函数：from度娘


特点：在x=0处像阶跃函数
输入 z =w0x0+w1x1+…+wnxn

优化方法：梯度上升（下降）法：

沿梯度方向找极值
原理：由于我们的任务是求得经验损失函数的最小值，所以实际上是一个“下坡”的过程。在每一个点上，我们希望往下走一步（假设一步为固定值0.5米），使得下降的高度最大，那么我们就要选择坡度变化率最大的方向往下走，这个方向就是经验损失函数在这一点梯度的反方向。每走一步，我们都要重新计算函数在当前点的梯度，然后选择梯度的反方向作为走下去的方向。随着每一步迭代，梯度不断地减小，到最后减小为零。这就是为什么叫“梯度下降法”。

核心：每走一步，迭代一次，计算当前的梯度，选择梯度最大（最小）的方向走

图上的γ 是一个系数，我们这里用0.001，迭代次数设置为500次
训练算法和画图的方法代码如下

import numpy as np
import matplotlib.pyplot as plt
import random


def loadDataSet():
    """
    测试用例的加载,设有x1,x2,补上x0的常数系数为1方便计算
    :return:
    """
    dataMat=[];
    labelMat=[];
    fr =open("testSet.txt")
    for line in fr.readlines():
        lineArr=line.strip().split()
        dataMat.append([1,float(lineArr[0]),float(lineArr[1])])
        labelMat.append(int(lineArr[-1]))
    return dataMat,labelMat


def sigmoid(inX):
    return 1.0/(1+np.exp(-inX))


def gradAscent(dataMatIn,classLabels):
    """
    梯度上升法
    :param dataMatIn: 数据集
    :param classLabels: 标签集
    :return:
    """
    #先转化为 Numpy的矩阵
    dataMat =np.mat(dataMatIn)
    m,n=np.shape(dataMat)
    labelMat=np.mat(classLabels).transpose()

    #系数/循环次数
    k =0.001
    maxcycle =500
    weights =np.ones((n,1))
    for i in range(maxcycle):
        h = sigmoid(dataMat*weights)
        error = (labelMat-h)
        weights =weights+k*dataMat.transpose()*error
    #numpy.matrix.transpose 矩阵转置

    return weights


def plotBestFit(weights):
    """
    画图
    :param weights: w
    :return: 图像
    """
    dataMatIn,labelIn =loadDataSet()
    dataMat=np.mat(dataMatIn)
    n =np.shape(dataMat)[0]
    #这是数据个数(列)

    #分类，根据标签存在不同的数据集里
    xcord1=[]
    xcord2=[]
    ycord1=[]
    ycord2=[]

    for i in range(n):
        if(int(labelIn[i])==1):
            xcord1.append(dataMat[i,1])
            ycord1.append(dataMat[i,2])
        else:
            xcord2.append(dataMat[i,1])
            ycord2.append(dataMat[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()

需要注意的是，这里先把矩阵通过np.mat(）方法转化成了numpy矩阵，这样才可以使用部分numpy的函数
结果如下：

500次迭代，足以获得稳定的（正确的）结果了

方法2：随机梯度法

梯度上升（下降）法有一个问题，就是每次迭代（更新系数）都要遍历整个数据集，计算起来太费劲了，如果样本、特征都很多，计算复杂度就很高

随机梯度法是每次只用一个样本点来更新系数
代码如下：

def gradAscent2(dataMatIn,classLabels):
    """
    随机梯度上升算法
    :param dataMatIn: 数据集
    :param classLabels: 分类集
    :return:
    """
    # 先转化为 Numpy的矩阵
    dataMat=np.array(dataMatIn)
    m, n = np.shape(dataMat)

    k = 0.01
    weights =np.ones(n)

    for j in range(m):
        h=sigmoid(sum(dataMat[j]*weights))
        error =classLabels[j]-h
        weights =weights +k*error*dataMat[j]

    return weights

跑了一下，结果惨绝人寰

改进随机梯度上升法

def gradAscent3(dataMatIn,classLabels,numIter=150):
    """
    改进随机梯度上升算法
    :param dataMatIn: 数据集
    :param classLabels: 分类集
    :param numIter:迭代次数
    :return:
    """
    # 先转化为 Numpy的矩阵
    dataMat=np.array(dataMatIn)
    m, n = np.shape(dataMat)


    weights =np.ones(n)

    for i in range(numIter):
        dataindex =list(range(m))
        #先设置成最后
        for j in range(m):
            k =4/(1.0+i+j) +0.01
            #设置一个渐变，k值逐渐减少
            randindex = int(random.uniform(0,len(dataindex)))

            h=sigmoid(sum(dataMat[randindex]*weights))
            error =classLabels[randindex]-h
            weights =weights +k*error*dataMat[randindex]
            del(dataindex[randindex])

    return weights

解读如下：
1.k的值不再是一个常量，是一个随着迭代次数增加而变小的量（且永远不会减到0）/且避免了k的值严格下降
2.通过先建立一份List index清单，随机选择拿来做比较的样本点，避免了如上一段代码一样，每次都是从第一个数据到最后一个数据，避免了周期性波动。
注意：书中代码本实是

dataindex =range(m)

而在python3中，会报错

TypeError: 'range' object doesn't support item deletion

查阅资料发现，这是由于在python2里面返回的是一个数组对象
而python3中 返回的是一个range对象
改为

dataindex =list(range(m))

就解决了问题
每一组迭代，都选用[0,dataindex)里面选一个随机数,拿该点来处理数据，之后，把这个点删掉，下次再剩下的里面选，这样每轮迭代都是每个数据用过一次

得到的拟合曲线看起来好了一点

实战案例 疝气病症预测死马率

问题1：丢失的数据怎么处理？
方案有1>0
2>-1(异常值)
3>所有样品的均值
4>相似样本的均值
5>换算法

我最初觉得方案4最好，但是计算量会比较大
书中给的是：将空白值用0代替
这是由于weights的计算过程中
w = w0x0+w1x1+…+wnxn
w =w+kerrordataMatrix[randinx]
如果某个特征值为0，
w =w 不会对数据造成影响
而且这样的方案比较简单

给定的数据，缺失值已经用0补上，因此并没有实际的工作量
代码如下：

def colicTest():
    '''
    测试函数
    :return:
    '''
    #导入训练集
    fr =open("D:\\pyt_example\\01\\venv\\horseColicTraining.txt")
    trainMat=[]
    labelMat=[]

    for line in fr.readlines():
        curline =line.strip().split('\t')
        lineArr =[]
        for i in range(21):
            lineArr.append(float(curline[i]))

        trainMat.append(lineArr)
        labelMat.append(float(curline[-1]))


    weights =gradAscent3(trainMat, labelMat, 500)


    '''然后处理测试集'''
    errorcount=0.0
    numTest =0.0
    frtest =open("D:\\pyt_example\\01\\venv\\horseColicTest.txt")

    for line in frtest.readlines():
        numTest+=1.0
        curline =line.strip().split('\t')
        lineArr =[]
        for i in range(21):
            lineArr.append(float(curline[i]))

        testans = classifyvector(lineArr ,weights)

        if int(testans) !=int(curline[-1]):
            errorcount+=1.0


    errotrate = errorcount/numTest*100

    print("errotrate is %.2f%%" %errotrate)
    return errotrate

运行结果的错误率还是很高的：

>>> reload(Logistic)
<module 'Logistic' from 'D:\\pyt_example\\01\\venv\\Logistic.py'>
>>> Logistic.colicTest()
errotrate is 34.33%
34.32835820895522
>>>

书中最后设置了一个函数，调用该方法10次计算了平均值。
在数据集和测试集一致的情况下，得到的误差在30%-40%之间
平均误差35%
我觉得这还是迭代的次数过少不够稳定吧

PS:改为应用梯度上升法，得到的误差约28%，要稍微好一点

书中讲第七章，会继续使用这些数据，因此下一篇更新可能会先学习第七章