自我代码提升之逻辑回归

(点击上方数据取经团，可快速关注)

作者：数据取经团-王嘉琪

      对于数据挖掘或是机器学习等相关算法的学习，用代码将之实现是从理论到实践的有效途径，而目前在数据科学应用领域，很多成熟的工具（诸如python的sklearn模块或是各种各样的R包），为我们的学习和工作提供了便利。但这并不意味着自己动手去实现算法失去了意义，在学习阶段，自我去实现算法的过程不仅仅是对理论的一种提升，更是对自己coding实践能力的一种历练。

  在数据建模过程中，就有机器学习的监督学习而言，通常会有一个 目标变量（因变量）和若干个可能会对目标变量产生影响的变量（自变量），选用一个模型，结合训练数据集（XY已知）进行训练拟合，再用模型去进行预测。在模型的选择中，Y的分布类型可能会对模型的性质产生影响，通常情况下，因变量Y为正态分布的数值型变量时，采用的是回归预测方式，如线性回归等；因变量Y为贝努利分布的类别变量时，所采用的模型应该是一个分类器，而Logistic回归就是其中之一。

Logistic回归介绍

  二分类问题的Logistic模型的基本形式如下：

  因变量Y可能属于0或者1两个类别，n为自变量的数量，输出的结果为Y属于1类的概率，概率值P只能介于0和1之间，当P越接近于1，则说明Y越可能属于1类，否则，Y就越可能属于0类。

  就逻辑曲线的形状特征而言，越远离曲线的中心点（自然对数的指数结果为0，P等于0.5），那么自然对数指数每发生单位变动，P的变动幅度则会趋缓；当越接近中心点时，那么自然对数指数每发生单位变动，P的变动幅度则会变得陡峭。

  作为比较简单的几类分类器算法之一，逻辑回归具有一定的自身特点。首先，逻辑回归的形势比较简洁，训练结果的系数值α可以清楚地知道对应的自变量对因变量Y的影响，而输出的P可以作为分类的概率值，容易理解其含义；其次，逻辑回归属于结构简单的机器学习模型，相比于复杂的模型，训练效率非常高，具有一定的“性价比”；再次，逻辑回归可以同时适用于类别型和数值型的自变量。然而，也有一些缺陷，如模型训练会导致自变量多重共线性而影响泛化能力，非线性变化对概率值的区分度的问题等。

  对于逻辑回归的系数α（包含了截距项和自变量系数）的求解，首先用最大似然法估计模型参数，将每个样本个体的Xj和α写成向量的形式以方便计算和理解。得到了似然函数L(α)，而最后的目标则是使得似然函数最大化。

  对似然函数取对数得到对数似然函数，从损失函数的角度去理解，在L前加负号亦可以称之为对数损失函数（区别于回归方法的平方损失函数），并且将似然函数写成矩阵运算形式。

  根据最大似然估计的思想，我们需要使得对数似然函数最大化，可以理解成求解出最优化（max(L)）的组合α，所以可以采用梯度上升的方法来求解。梯度上升的思想是：在求解最大化问题的时候，函数输入点的取值每次沿着当前的函数正梯度方向，以某个距离值进行变化，经过若干次迭代以后达到最大值的点。

  所以对α的求解，应该先确定向量α的初始值（可随机确定，期望为0），然后再每次迭代过程中，计算出当前的对数似然函数L的梯度，再用梯度乘以一个距离s，即为本次迭代α的变化

PYTHON的逻辑回归简单实现

  在了解逻辑回归的求解过程推导之后，可以采用PYTHON去进行实现，主要是基于矩阵处理的numpy模块，对于测试采用的数据集，这里采用疝病数据预测马匹的死亡率的相关数据集，包含了训练数据集合测试数据集，前者用来训练模型，后者用以测试。

  首先给出读取数据集的函数，并且将数据集读取，划分自变量矩阵和因变量向量：

#数据集读取函数
def get_data(file_name,sep='\t'):
    data_open = open(file_name)
    lines = data_open.readlines()
    data = []    
    for line in lines:
        line = line.replace('\n','')
        line = line.split(sep)
        a=[]        
        for value in line:
            value = float(value)
            a.append(value)
        data.append(a)    
    return np.array(data)
#训练集读取
horse_train = get_data(u'C:\\Users\\intern\\Desktop\\example\
\机器学习算法实现最终版\\horseColicTraining.txt')
horse_train_x = horse_train[:,range(horse_train.shape[1]-1)]
horse_train_y = horse_train[:,horse_train.shape[1]-1]
#测试集读取
horse_test = get_data(u'C:\\Users\\intern\\Desktop\\example\
\机器学习算法实现最终版\\horseColicTest.txt')
horse_test_x = horse_test[:,range(horse_test.shape[1]-1)]
horse_test_y = horse_test[:,horse_test.shape[1]-1]

  训练数据集包含了21个维度的自变量，样本数为299，测试数据集包含了67条记录。在导入数据之后，开始对逻辑回归模型的函数进行编写：

#梯度上升求解LR#L= (Yln(1/(1+exp(-XK)))+(1-Y)*ln(1-1/(1+exp(-XK)))) ==> d(L)/d(K) = X.T•(Y-1/(1+exp(-XK)))#K_n+1 = K_n + d(L)/d(K_n) * s
class Logistic_Regression():
    def log_fun(self,x): #sigmoid函数
        return 1/(1+np.exp(-x))        
    def fit_gd(self,x,y,iter=1000,s=0.01,e_break=0.001): #逻辑回归参数求解
        x = np.matrix(np.hstack((np.ones((x.shape[0],1)),x)))
        length = x.shape[1]
        K_init = np.matrix(np.random.rand(length)-0.5).T #初始化系数
        y = np.matrix(y).T
        i = 1
        while i <= iter:#梯度上升迭代
            g = x.T.dot(np.array((y - self.log_fun(x.dot(K_init)))))
            K_init = K_init + g*s
            i += 1
        self.a = K_init
        return self.a        
    def predict_prob(self,x): #预测函数（概率）
        return np.array(self.log_fun(x.dot(self.a[1:,:])+self.a[0,0]).T)[0,:]    
    def predict_type(self,x,thre_var=0.5): #预测函数（类别）
        type_pre = self.predict_prob(x)        
        for i in range(type_pre.shape[0]):            
            if type_pre[i] > thre_var:
                type_pre[i] = 1
            else:
                type_pre[i] = 0
        return type_pre        
    def predict_accuracy(self,x,y,thre_var=0.5): #预测准确度函数
        y_pre = self.predict_type(x,thre_var)
        accuracy = sum(y_pre==y)/float(y.shape[0])        
        return accuracy

  这里将模型定义为一个类，类中包含了四个函数，sigmoid函数（用于其它函数的调用），梯度求解训练函数（主体，输入数据集因变量和自变量来训练模型，输出模型系数），概率预测函数和类别预测函数。然后，我们开始用疝病数据集去训练函数，设置的迭代次数为10000次，上升步长为0.001，并且将训练完成的模型对测试集进行类别预测。

logistic_model = Logistic_Regression() #实例化模型对象
logistic_model.fit_gd(horse_train_x,horse_train_y,iter=10000,s=0.001) #模型训练
horse_pred = logistic_model.predict_type(horse_test_x) #预测
logistic_model.predict_accuracy(horse_test_x,horse_test_y) #准确度

       输出预测准确度，为0.73。

利用Sklearn机器学习模块训练逻辑回归模型

  Python的Sklearn模块提供了强大的机器学习建模的功能，包含了当前主流的机器学习各类模型，下列代码是对该模块进行逻辑回归建模的简单应用，有关该模块在调参，测试等深入应用请读者自行研究。

#采用sklearn模块的逻辑回归示例
from sklearn.linear_model import LogisticRegression #加载模块
lr = LogisticRegression(C=1000, random_state=0) #实例化
lr.fit(horse_train_x,horse_train_y) #模型训练
sum(lr.predict(horse_train_x) == horse_train_y)/float(len(horse_train_y)) #分类准确度

拓展1：牛顿法求解逻辑回归

  牛顿法也是比较常用的优化求解算法之一，其基本思想是，用目标函数在某一点的二阶泰勒展开式作为其替代函数，为了求出极值点（满足优化目标），令目标函数一阶导为0，对应地，其替代函数（泰勒展开式）的一阶导为0，对替代函数求解方程可以得到新的点，改点比原先的那个点更加接近优化目标，所以通过这样迭代若干次找到最优点。

  和梯度下降方法相比，牛顿法的二阶收敛性质使得其在实现优化的速度上要远远高于梯度下降的方法。但是在牛顿法的各个步骤中，对二阶导海塞矩阵H的求逆操作在高阶情况下会严重影响效率，所以通常采用一个正定矩阵G来对H进行替代，即拟牛顿法

  本文对牛顿法和DFP拟牛顿法求解逻辑回归给出了PYTHON的实现代码，对这些方法的具体论述可以参见李航的《统计学习方法》。在代码运行之后我们可以看出，牛顿法只需要较少的迭代次数，就使得模型达到了较好的预测效果。

  对牛顿法的实现参见下图：

#附：拟牛顿法求解LogisticRegression
#相比于梯度下降对模型的求解（一阶收敛），牛顿法，拟牛顿法具有更快的下降收敛速度（二阶收敛）
#牛顿法是利用目标函数的二阶泰勒展开式作为其替代函数，
#并且得到对应的梯度，按照极值点梯度为0的思想，求出极值点xi
#用新的xi会带入二阶泰勒展开式的梯度中，重复迭代，直到停止

#因为牛顿法涉及二阶展开的海塞矩阵及其求逆运算，该步骤较为复杂
#通过对海塞矩阵H的逆矩阵的估计来改进和简化的方法为拟牛顿法
#DFP作为常见的拟牛顿法之一的算法，采用一个满足条件的正定矩阵G来替代H的逆
#每次迭代通过两个附加项P和Q来更新G
class logistic_NTmodel():
    def log_fun(self,x): #sigmoid函数
        return 1/(1+np.exp(-x))
    def pos_def_ma(self,n): #随机产生一个正定矩阵（用于拟牛顿法）
        import numpy as np
        a = np.diag(np.random.rand(n))
        c = np.array(np.random.rand(n*n)).reshape(n,n)       
        a = c.T.dot(a).dot(c)
        return(a)    
    def grad_los(self,w,x,y): #logistic损失函数的一阶导数（梯度）
        grad = -x.T.dot(np.array(y - self.log_fun(x.dot(w))))        return(grad)
    def hessi_mat(self,w,x): #海塞矩阵函数
        hessi_mat = np.zeros([x.shape[1],x.shape[1]])
        x_array=np.array(x)    
        for i in range(x.shape[1]):           
            for j in range(x.shape[1]):
                hessi_mat[i,j] = (np.matrix(x_array[:,i]*x_array[:,j]).dot(np.array(
                        self.log_fun(x.dot(w)))*np.array(1-self.log_fun(x.dot(w)))))[0,0]        
        return(np.matrix(hessi_mat))
    def fit_NT(self,x,y,iter=100): #牛顿法拟合函数
        import numpy as np   
        x = np.matrix(np.hstack((np.ones((x.shape[0],1)),x)))
        y = np.matrix(y).T
        length = x.shape[1]        
        #K_init = np.matrix(np.random.rand(length)-0.5).T #初始化参数
        K_init = np.matrix(np.zeros(length)).T #初始化参数（0向量）
        i = 1
        while i <= iter:
            gi = self.grad_los(K_init,x,y) #产生梯度
            Hi = self.hessi_mat(K_init,x) #产生海塞举证
            #print Hi
            pi = -Hi.I.dot(gi) #计算p
            K_init = K_init + pi #更新模型系数
            i += 1
        self.a = K_init
        return(self.a)    
    def fit_DFP(self,x,y,iter=100,a = 0.01): #DFP拟牛顿算法函数(p的系数取固定值a)
        import numpy as np   
        x = np.matrix(np.hstack((np.ones((x.shape[0],1)),x)))
        y = np.matrix(y).T
        length = x.shape[1]
        G = self.pos_def_ma(x.shape[1])  #初始正定矩阵G（用以替代牛顿法海塞矩阵（逆矩阵））
        #K_init = np.matrix(np.random.rand(length)-0.5).T #初始化参数
        K_init = np.matrix(np.zeros(length)).T
        i = 1
        gi = self.grad_los(K_init,x,y) #梯度计算
        while i <= iter:        
            pi = -G.dot(gi)
            K_init_old = K_init.copy()
            K_init = K_init + a*pi
            gi_old = gi.copy() 
            gi = self.grad_los(K_init,x,y) #梯度更新计算
            gi_dif = gi - gi_old
            K_init_dif = K_init - K_init_old
            Pi = K_init_dif.dot(K_init_dif.T)/(K_init_dif.T.dot(gi_dif)) #正定矩阵G附加项1
            Qi = -G.dot(gi_dif).dot(gi_dif.T).dot(G)/(gi_dif.T.dot(G).dot(gi_dif)) #正定矩阵G附加项2
            G = G + Pi + Qi #更新正定矩阵G
            i += 1
        self.a = K_init
        return(self.a)        
    def predict_prob(self,x): #预测函数（概率）
        return np.array(self.log_fun(x.dot(self.a[1:,:])+self.a[0,0]).T)[0,:]    def predict_type(self,x,thre_var=0.5): #预测函数（类别）
        type_pre = self.predict_prob(x)
        for i in range(type_pre.shape[0]):
            if type_pre[i] > thre_var:
                type_pre[i] = 1
            else:
                type_pre[i] = 0
        return type_pre        
    def predict_accuracy(self,x,y,thre_var=0.5): #预测准确度函数
        y_pre = self.predict_type(x,thre_var)
        accuracy = sum(y_pre==y)/float(y.shape[0])        
        return accuracy

拓展2：正则化的探究

  一般的逻辑回归可能会因为输入的自变量存在多重共线性而影响模型的预测效果，可以在特征工程中对变量进行筛选和降维来避免对模型输入过多的变量，此外，还可以采用正则化的方式来处理。

  常用的正则化方法包括了L1正则化和L2正则化方法，分别是在逻辑回归的损失函数（对数似然函数前加负号）后追加一个L1正则项（各个系数α绝对值和）或L2正则项（各个系数α的平方和）乘以一个系数。从而在优化过程中使得部分α的值趋向于0。

  和L1正则化相比，L2方法不仅可以有效地防止模型的过拟合，提升模型的预测泛化能力，还能使得模型在优化求解的过程中更加稳定和迅速，在统计学的角度，L2正则化方法又可以称为岭回归。本文的代码中对L2正则化方法进行了简单的实现，在梯度下降求解逻辑回归的类中加入下列代码，请读者自行测试。

    def fit_gd_L2(self,x,y,iter=1000,s=0.01,l2=0.1): #逻辑回归参数求解(L2)
        x = np.matrix(np.hstack((np.ones((x.shape[0],1)),x)))
        length = x.shape[1]
        K_init = np.matrix(np.random.rand(length)-0.5).T #初始化参数
        y = np.matrix(y).T
        i = 1
        while i <= iter:#梯度下降迭代(添加了L2正则项)    
            K_init = K_init + x.T.dot(np.array(y - self.log_fun(x.dot(K_init))))*s - l2*K_init*s
            i += 1
        self.a = K_init
        return self.a

本文代码及数据获取方式：下图扫码关注公众号，在公众号中回复逻辑回归即可~

联系作者：当你对文章有任何疑问和建议可以在公众号直接发消息，我们看到都会回复哒，一起交流数据~