分类和Logistic回归

      在上一章节中，介绍了简单的线性回归，给出一系列的离散点，利用回归模型找到一条最佳的拟合直线，其中在求解最佳拟合直线的过程中利用到了批梯度下降算法和随机梯度下降算法以及最小二乘法。可以看到，回归是一个连续的模型，那么怎样将这样一个连续的模型用在分类问题上呢，这就是这一节中将要介绍的对数回归模型。

     对数回归本质上就是线性回归，只不过在线性回归的基础上加上了一个函数而已，这一点和单层感知器和相似，在单层感知器的输入上加上一个激励函数，就会得到不同的输出。对数回归也是一样，只需要在我们上一节中的回归模型的基础上加上一个函数，就可以将上面的线性回归模型应用到分类问题上。

1.Sigmoid函数

      Sigmoid函数也叫做逻辑函数，形式如下所示：

                                                     

     函数对应的图像如下所示：

                                       

      可以看到，对于不同的输入z，函数的输出总是在0到1之间，当z的值为0是，函数值是0.5。正是因为函数的这样的性质我们可以将线性回归中的输出值作为z作用大函数上，最终的结果就是一个在0到1之间的数字，以0作为分界线，当线性回归的结果大于0的时候，分类结果就是1，当线性回归的结果小于0的时候，分类的结果就是0。有了这种思想，我们就可以将回归函数重新定义如下：

                                            

2 模型定义

       对数回归可以用来分类0/1问题，也就是预测结果属于0或者1的二分类问题，我们假设问题满足伯努利分布，也就是如下的形式：

                                        

      其中是我们将要求取的参数，将上述的式子用一个式子表示就是如下的形式：

                                        

     假设所有的样本都是独立分布的，利用最大似然估计

                                        

     对式子两边同时取对数，得到

                                

    现在式子中只有一个是未知的，现在的目标是要将式子最大化，采用梯度上升的做法

                                

     对单一的样本的似然函数进行求导得：

                

    可以看到，这里的结果和线性回归中的结果是一样的，所以采用梯度上升的算法对进行更新的过程为：
                   

3 实验过程

    本次实验中采用的数据集如下：

                 

 3.1 批梯度上升算法实现过程           

#批梯度上升算法def gradAscent(X,Y):    X=np.array(X)    m,n=np.shape(X)    Y=np.array(Y).reshape(m,1) #转化为列向量的形式    alpha = 0.01               #定义步长    iteration=400    thea = np.ones((n, 1)) #定义初始的点    for k in range(iteration):        H=sigmoid(X.dot(thea))        error=Y-H        thea=thea+alpha*X.T.dot(error)    return thea

   实验结果：

         

 3.2 随机梯度上升算法      

#随机梯度算法def stocgradAscent(X,Y):    X=np.array(X)    m,n=np.shape(X)    Y=np.array(Y).reshape(m,1)    thea=np.ones((n,1))    alpha = 0.01  # 定义步长    for i in range(m):        H=sigmoid(X[i].dot(thea))        error=Y[i]-H        thea=thea+alpha*error*X[i].reshape(n,1)    return thea

    实验结果如下：

      

    可以看到，随机梯度上升算法的结果在本例中并不比批梯度上升算法的结果好，但是只从这样的结果中直接的去比较算法的好坏是不公平的，在批梯度上升中，算法在整个的数据集上经过了 400次的迭代过程。而判断一个算法好坏的可靠方法是看它是否收敛，也就是是在经过有限步骤的迭代之后，参数是否到达了一个稳定值。在《机器学习实战》一书中，针对这个例子给出了系数迭代过程中变化趋势图：

          

   可以看到X的系数在达到稳定的时候经过的次数是不一样的，在大的波动停止后，还会有一些小的周期性波动，产生这种现象的原因是存在着一些不能正确分类的样本点，在每次迭代的时候会引发系数的剧烈改变，我们针对随机梯度上升算法做出如下的该进：

  1.在每次迭代的过程中，改变步长alpha的值，这样会减小波动性。

  2.随机的选择样本点进行梯度上升来更新系数值，这样可以避免周期性的波动，在每次迭代完成之后从列表中删除这
  个已经选择的值。

该进之后的随机梯度上升算法如下：

#改进的随机梯度下降上升算法def stocgradAscent1(X,Y):    X=np.array(X)    m,n=np.shape(X)    Y=np.array(Y).reshape(m,1)    iteration = 400    thea=np.ones((n,1))    for j in range(iteration):        dataindex=[]        for k in range(m):            dataindex.append(k)        for i in range(m):            alpha=4/(1.0+i+j)+0.01            randinx=int(random.uniform(0,len(dataindex)))            h=sigmoid(X[i].dot(thea))            error=Y[i]-h            thea=thea+alpha*error*X[i].reshape(n,1)            dataindex.remove(dataindex[randinx])    return thea

  该进之后的效果如下：

          

  可以看到，在改进之后，随机梯度上升算法一样可以取得很好的效果。

  文中实验所有的代码如下：

#对数回归import numpy as npimport matplotlib.pyplot as pltimport random#读取文本数据，三维样本def ReadFile(filename):    X=[]    Y=[]    lines=open(filename,encoding='utf-8').readlines()    for i in range(len(lines)):        line=lines[i].strip().split()        X.append([1.0,float(line[0]),float(line[1])])        Y.append(int(line[2]))    return X,Y#定义sigmoid函数def sigmoid(x):    return 1.0/(1+np.exp(-x))#批梯度上升算法def gradAscent(X,Y):    X=np.array(X)    m,n=np.shape(X)    Y=np.array(Y).reshape(m,1) #转化为列向量的形式    alpha = 0.01               #定义步长    iteration=400    thea = np.ones((n, 1)) #定义初始的点    for k in range(iteration):        H=sigmoid(X.dot(thea))        error=Y-H        thea=thea+alpha*X.T.dot(error)    return thea#随机梯度算法def stocgradAscent(X,Y):    X=np.array(X)    m,n=np.shape(X)    Y=np.array(Y).reshape(m,1)    thea=np.ones((n,1))    alpha = 0.01  # 定义步长    for i in range(m):        H=sigmoid(X[i].dot(thea))        error=Y[i]-H        thea=thea+alpha*error*X[i].reshape(n,1)    return thea#改进的随机梯度下降上升算法def stocgradAscent1(X,Y):    X=np.array(X)    m,n=np.shape(X)    Y=np.array(Y).reshape(m,1)    iteration = 400    thea=np.ones((n,1))    for j in range(iteration):        dataindex=[]        for k in range(m):            dataindex.append(k)        for i in range(m):            alpha=4/(1.0+i+j)+0.01            randinx=int(random.uniform(0,len(dataindex)))            h=sigmoid(X[i].dot(thea))            error=Y[i]-h            thea=thea+alpha*error*X[i].reshape(n,1)            dataindex.remove(dataindex[randinx])    return thea#作出分类曲线def plot(X,Y,thea):    X = np.array(X)    m, n = np.shape(X)    Y = np.array(Y).reshape(m, 1)    xcord1=[];ycord1=[]    xcord2=[];ycord2=[]    for i in range(m):        if Y[i]==1:            xcord1.append(X[i][1])            ycord1.append(X[i][2])        if Y[i]==0:            xcord2.append(X[i][1])            ycord2.append(X[i][2])    fig=plt.figure()    ax=fig.add_subplot(111)    type1=ax.scatter(xcord1, ycord1, s=30, c='red', marker='s',label='1')    type2=ax.scatter(xcord2, ycord2, s=30, c='green',marker='s',label='0')    plt.legend(loc='upper right')    x=np.arange(-4,4)    y=((-thea[0] - thea[1] * x) / thea[2])    ax.plot(x, y)    plt.xlabel("X1")    plt.ylabel("X2")    plt.show()if __name__=="__main__":    X,Y=ReadFile('ex2data1.txt')    #thea=stocgradAscent1(X,Y)    thea=gradAscent(X,Y)    #thea=stocgradAscent(X,Y)    print(thea)    plot(X,Y,thea)

4 利用Logistic解决多分类问题

     Logistic主要应用在二分类问题中，我们上述的实验就是在一个二分类的问题上进行的，Logistic同样也可以应用在多分类问题中。在多分类问题中，我们可以将多分类问题看作是多个二分类问题。在本实验中，采用传统的手写字识别的数据集，利用Logistic回归实现手写字识别的思想大致如下：这是一个包含10个类别的分类问题，将整个识别问题看作为多个二分类问题。例如在识别数字0的时候，我们可以将整个问题看作为0和非0处理，将属于类别0的数据项看作是1，将不属于数字0的数据项看作是0，对于其他的类别做同样的处理。在二分类问题中，求得的是一个向量的形式，而在多分类问题中，因为将问题看作是多个分类问题，所以训练的结果是多个分类器，最终得到的是一个矩阵，矩阵的每一行代表一个分类器。

     实验的结果如下：

       

    本实验中所有代码如下：

#利用logistic回归实现多分类问题import numpy as npfrom os import listdir#定义数据读取函数def LoadData(dir):    filelist=listdir(dir)    m=len(filelist)    X=np.zeros((m,1024))    Y=np.zeros((m,1))    for i in range(m):        arr=np.zeros((1,1024))        filename=filelist[i]        fr=open('%s/%s' %(dir,filename))        for j in range(32):            line=fr.readline()            for k in range(32):                arr[0, 32 * j + k] = int(line[k])        X[i,:]=arr        name=filename.split('.')[0]        label=name.split('_')[0]        Y[i]=int(label)    return X,Y#定义Sigmoid函数def sigmoid(x):    return .5 * (1 + np.tanh(.5 * x))    #return 1.0/(1+np.exp(-x))#利用批梯度上升算法为每个类别训练一个分类器def gradAscent(X,Y):    m,n=np.shape(X)    num_labels=10                         #类别的数目    all_thea=np.ones((num_labels,n))      #定义将要返回的系数矩阵    alpha = 0.01                          #定义步长    iteration=400    for i in range(num_labels):        index = []        thea = np.ones((n, 1))           #定义初始的点        for j in range(m):            if Y[j]==i:                index.append(1)            else:                index.append(0)        y=np.array(index).reshape(m,1)        for k in range(iteration):            H=sigmoid(X.dot(thea))            error=y-H            thea=thea+alpha*X.T.dot(error)        all_thea[i,:]=thea.reshape(1,n)    return all_thea#定义预测函数def predict(X,thea):    m,n=np.shape(X)    num_labels = 10  # 类别的数目    lables=[]    r=sigmoid(X.dot(thea.T))    for i in range(m):        lables.append(np.argmax(r[i]))    return lablesif __name__=="__main__":    TrainX,TrainY=LoadData("trainingDigits")    TestX,TestY=LoadData("testDigits")    all_thea=gradAscent(TrainX,TrainY)    labels=predict(TestX,all_thea)    y=np.array(labels).reshape(np.shape(TestX)[0],1)    k=0    for i in range(np.shape(TestX)[0]):        print('%d的识别结果是:%d'%(TestY[i],y[i]))        if y[i]==TestY[i]:            k=k+1    print('true rate is:','%.4f' %(k/(np.shape(TestY)[0])))

参考文献

   Ng机器学习讲义

   http://www.cnblogs.com/netuml/p/5727745.html

   http://blog.csdn.net/stdcoutzyx/article/details/9113681