<学习笔记> Logistic回归

2017-12-13

ps:感谢《机器学习实战》指导

这次的内容和最优化有关系，提到最优化，我们其实见过很多了，我们的原则就是用最少的付出换取
最丰厚的汇报。比如如何设计发动机使油耗最少而功率最大；如何规划路线使我们用最短的时间从始发点到
目的点；如何用最少的背包装最多的东西….
接下来我们需要稍微理解一下回归这个词。
假设现在我们有一些数据点，我们用一条直线对我们这些数据点进行拟合，这个过程就是回归。
我们的这条线可以写成一个函数，这个过程当中，我们通过不断地调整我们函数中的各种参数，使得我们的
数据点们平均最大化的接近我们的这条直线。
Logistic回归的思想：
根据现有数据对分类边界线建立回归公式，以此进行分类。
优点：
计算代价不高，易理解易实现；
缺点：
容易欠拟合，分类的精度可能不高；
适用类型：数值型和标称型的数据类型。
我们想要得到的函数是能够接受所有的输入然后做出分类。
在两个类的情况下，我们很容易想到一个函数，赫维赛德跃迁函数：

连续形式：
H(x)={0,if x<0;1,x>0;1/2,x=0}
离散形式：
Hn={0,n<0;1,n>=0;n为整数}

函数可以接受所有输出然后输出0或1，这样我们可以测出类别。
但是这个函数在0和1跃迁的这个瞬间过程不容易处理，所以我们选择另一个函数,Sigmoid函数：

1	σ(z)=1/(1+exp(-z))

根据这个函数的图像，我们可以做出以下分析：
当z为0时，我们的函数值为0.5；
随着z的增大，函数值趋近于1，随着z的减小，函数值趋近于0；如果我们从一个很大的尺度看这个函数的图像，
他其实很像一个跃迁函数。
为了实现我们的分类器，我们可以在每个特征上都乘以一个回归系数，然后把所有得到的结果求和，代入到
我们的Sigmoid函数，得到一个范围在0，1之间的数。
然后将大于0.5的归于“1类”,将小于0.5的归于“0类”.
确定分类函数形式后，我们下一个要解决的问题是确定最佳回归系数。
Sigmoid函数的输入记为z，可由下列公式得出：

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

将上述公式写成向量的形式：z=w^Tx,即将w转置乘x。
向量x为我们的输入，向量w即是我们要找的最佳参数。
为了寻找最佳参数，我们要用到最优化理论的一些知识。
1.梯度上升法
主要思想：如果想要找到某个函数的最大值，那就沿着该函数的梯度方向探寻。
在高等数学中，函数f(x,y)的梯度看作一个向量，记作：

1	grad(f)=(∂f/∂x,∂f/∂y)

函数f必须在待计算的点上有定义并且可微。
梯度上升算法每到达一个点的时候都会重新估计下次的移动方向。
梯度算子总是指向函数增长最快的方向。我们将每次在方向上的移动量称为步长，记作α，如果用向量的方式表示，
我们的的梯度算法的迭代公式为：

1	w:=w+αgrad(f(w))

这个公式在达到终止条件前会一直迭代。
以上就是我们进行分类器构造的理论工具。
接下来我们用python来实现我们的算法。

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(intX):
    from math import exp
    return 1.0/(1+exp(-intX))
def grandascent(datamatin,classlabels):
    from numpy import *
    datamatrix=mat(datamatin)
    labelmatrix=mat(classlabels).transpose()
    #转换为numpy矩阵的数据类型；
    m,n=shape(datamatrix)
    alpha=0.001
    maxcycles=500
    weights=ones((n,1))
    for k in range(maxcycles):
        h=sigmoid(datamatrix*weights)
        error=labelmatrix-h
        weights=weights+alpha*datamatrix.transpose()*error
    return weights
def plotbestfit(wei):
    import matplotlib.pyplot as plt
    from numpy import *
    weights=wei.getA()
    datamat,labelmat=loaddataset()
    dataarr=array(datamat)
    n=shape(dataarr)[0]
    xcords1=[]
    ycords1=[]
    xcords2=[]
    ycords2=[]
    for i in range(n):
        if int(labelmat[i])==1:
            xcords1.append(dataarr[i,1])
            ycords1.append(dataarr[i,2])
        else:
            xcords2.append(dataarr[i,1])
            ycords2.append(dataarr[i,2])
    fig=plt.figure()
    ax=fig.add_subplot(111)
    ax.scatter(xcords1,ycords1,s=30,c='red',marker='s')
    ax.scatter(xcords2,ycords2,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()

随机梯度上升算法
我们之前的梯度上升算法在每次更新回归系数的时候都要遍历整个数据集，在数据集很大的时候
这种做法就会变得笨拙。我们通过下面这种方式改进我们的算法。
每次只使用一个数据点对回归系数进行更新，该方法称为随机梯度上升算法。
实现随机梯度上升算法的python实现：

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.01
            #调整alpha值来缓和波动
            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

接下来我们将理论运用到实践当中。
实例：从疝气病症预测病马的死亡率
研究过程：
（1）.收集数据：这里使用我们预备好的数据文件
（2）.准备数据：用python解析文本，填充空缺
（3）.分析数据：将我们的数据结果可视化并进行观察总结
（4）.训练算法：使用优化算法，找到最佳的数据
（5）.测试算法：观察错误率，通过改变迭代次数和步长来优化我们的参数，在过程当中根据错误率来决定我们是否要继续进行训练
（6）.使用算法：实现简单的真正的从症状预测结果
分析过程要面对的一些问题
1.数据缺失：鉴于我们所得到的数据也许有时会很来之不易，如果重新进行收集可能会给我们带来很大的损失，
所以我们要考虑一些其他的方法：
（1）.使用特征值的均值填补
（2）.使用我们自设定的特殊值
（3）.干脆忽略掉存在缺失的小部分样本
（4）.使用相似的样本抽取特征填补
（5）.使用其他的算法预测并填补
结合我们这里的问题，我们决定选择实数0来填补，这样做适用于logistic回归。
如若发现类别标签缺失，直接舍弃。
使用logistic方法进行分类：将测试集上的每个特征向量乘最优化得到的回归系数，再将该结果求和，然后作为sigmoid函数的输入即可。
接下来继续完善我们的分类器：

def classifyvector(intx,weights):
    prob=sigmoid(sum(intx*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)
    errcount=0
    numtestvec=0.0
    for line in frtest.readline():
        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]):
            errcount+=1
    errrate=(float(errcount/numtestvec))
    print("the errorrate is%d",errrate)
    return errrate
def multitest():
    numtests=10
    errorsum=0.0
    for k in range(numtests):
        errorsum+=colictest()
    print("after 10 iterations the average error rate is %f",float(errorsum/numtests))