逻辑回归模型介绍和程序实现

逻辑回归原理及实现

 　　虽然叫做“回归”，但是这个算法是用来解决分类问题的。回归与分类的区别在于：回归所预测的目标量的取值是连续的（例如房屋的价格）；而分类所预测的目标变量的取值是离散的（例如判断邮件是否为垃圾邮件）。当然，为了便于理解，我们从二值分类（binary classification）开始，在这类分类问题中，y只能取0或1。更好的理解问题，先举个小例子：假如我们要制作一个垃圾邮件过滤系统，如果一封邮件是垃圾系统，y=1，否则y=0 。给定训练样本集，当然它们的特征和label都已知，我们就是要训练一个分类器，将它们分开。

         回归分析用来描述自变量x和因变量Y之间的关系，或者说自变量X对因变量Y的影响程度，并对因变量Y进行预测。其中因变量是我们希望获得的结果，自变量是影响结果的潜在因素，自变量可以有一个，也可以有多个。一个自变量的叫做一元回归分析，超过一个自变量的叫做多元回归分析。下面是一组广告费用和曝光次数的数据，费用和曝光次数一一对应。其中曝光次数是我们希望知道的结果，费用是影响曝光次数的因素，我们将费用设置为自变量X，将曝光次数设置为因变量Y，通过一元线性回归方程和判定系数可以发现费用(X)对曝光次数(Y)的影响。

1、 逻辑回归模型

    回归是一种极易理解的模型，就相当于y=f(x)，表明自变量x与因变量y的关系。最常见问题有如医生治病时的望、闻、问、切，之后判定病人是否生病或生了什么病，其中的望闻问切就是获取自变量x，即特征数据，判断是否生病就相当于获取因变量y，即预测分类。

    最简单的回归是线性回归，在此借用Andrew NG的讲义，有如图1.a所示，X为数据点——肿瘤的大小，Y为观测值——是否是恶性肿瘤。通过构建线性回归模型，如h_θ(x)所示，构建线性回归模型后，即可以根据肿瘤大小，预测是否为恶性肿瘤h_θ(x)≥.05为恶性，h_θ(x)<0.5为良性。

图1 线性回归示例

    然而线性回归的鲁棒性很差，例如在图1.b的数据集上建立回归，因最右边噪点的存在，使回归模型在训练集上表现都很差。这主要是由于线性回归在整个实数域内敏感度一致，而分类范围，需要在[0,1]。逻辑回归就是一种减小预测范围，将预测值限定为[0,1]间的一种回归模型，其回归方程与回归曲线如图2所示。逻辑曲线在z=0时，十分敏感，在z>>0或z<<0处，都不敏感，将预测值限定为(0,1)。

图2 逻辑方程与逻辑曲线

    逻辑回归其实仅为在线性回归的基础上，套用了一个逻辑函数，但也就由于这个逻辑函数，逻辑回归成为了机器学习领域一颗耀眼的明星，更是计算广告学的核心。对于多元逻辑回归，可用如下公式似合分类，其中公式(4)的变换，将在逻辑回归模型参数估计时，化简公式带来很多益处，y={0,1}为分类结果。 

    对于训练数据集，特征数据x={x₁, x₂, … , x_m}和对应的分类数据y={y₁, y₂, … , y_m}。构建逻辑回归模型f(θ)，最典型的构建方法便是应用极大似然估计。首先，对于单个样本，其后验概率为：

   

那么，极大似然函数为：

log似然是：

2、 梯度下降

    由第1节可知，求逻辑回归模型f(θ)，等价于：

   

采用梯度下降法：

     从而迭代θ至收敛即可：

 

 

2.1 梯度下降算法

梯度下降算法的伪代码如下：

################################################

初始化回归系数为1

重复下面步骤直到收敛{

        计算整个数据集的梯度

        使用alpha x gradient来更新回归系数

}

返回回归系数值

################################################

      注：因为本文中是求解的Logit回归的代价函数是似然函数，需要最大化似然函数。所以我们要用的是梯度上升算法。但因为其和梯度下降的原理是一样的，只是一个是找最大值，一个是找最小值。找最大值的方向就是梯度的方向，最小值的方向就是梯度的负方向。不影响我们的说明，所以当时自己就忘了改过来了，谢谢评论下面@wxltt的指出。另外，最大似然可以通过取负对数，转化为求最小值。代码里面的注释也是有误的，写的代码是梯度上升，注销成了梯度下降，对大家造成的不便，希望大家海涵。

 

2.2 随机梯度下降SGD (stochastic gradient descent)

      梯度下降算法在每次更新回归系数的时候都需要遍历整个数据集（计算整个数据集的回归误差），该方法对小数据集尚可。但当遇到有数十亿样本和成千上万的特征时，就有点力不从心了，它的计算复杂度太高。改进的方法是一次仅用一个样本点（的回归误差）来更新回归系数。这个方法叫随机梯度下降算法。由于可以在新的样本到来的时候对分类器进行增量的更新（假设我们已经在数据库A上训练好一个分类器h了，那新来一个样本x。对非增量学习算法来说，我们需要把x和数据库A混在一起，组成新的数据库B，再重新训练新的分类器。但对增量学习算法，我们只需要用新样本x来更新已有分类器h的参数即可），所以它属于在线学习算法。与在线学习相对应，一次处理整个数据集的叫“批处理”。

        随机梯度下降算法的伪代码如下：

###############################################

初始化回归系数为1

重复下面步骤直到收敛{

        对数据集中每个样本

               计算该样本的梯度

                使用alpha xgradient来更新回归系数

 }

返回回归系数值

###############################################

 

2.3 改进的随机梯度下降

　　1）在每次迭代时，调整更新步长alpha的值。随着迭代的进行，alpha越来越小，这会缓解系数的高频波动（也就是每次迭代系数改变得太大，跳的跨度太大）。当然了，为了避免alpha随着迭代不断减小到接近于0（这时候，系数几乎没有调整，那么迭代也没有意义了），我们约束alpha一定大于一个稍微大点的常数项，具体见代码。

　　2）每次迭代，改变样本的优化顺序。也就是随机选择样本来更新回归系数。这样做可以减少周期性的波动，因为样本顺序的改变，使得每次迭代不再形成周期性。

 改进的随机梯度下降算法的伪代码如下：

################################################

初始化回归系数为1

重复下面步骤直到收敛{

       对随机遍历的数据集中的每个样本

              随着迭代的逐渐进行，减小alpha的值

              计算该样本的梯度

              使用alpha x gradient来更新回归系数

    }

返回回归系数值

#################################################

 

3、Python实现

训练数据：

 1 -0.017612    14.053064    0 2 -1.395634    4.662541    1 3 -0.752157    6.538620    0 4 -1.322371    7.152853    0 5 0.423363    11.054677    0 6 0.406704    7.067335    1 7 0.667394    12.741452    0 8 -2.460150    6.866805    1 9 0.569411    9.548755    010 -0.026632    10.427743    011 0.850433    6.920334    112 1.347183    13.175500    013 1.176813    3.167020    114 -1.781871    9.097953    015 -0.566606    5.749003    116 0.931635    1.589505    117 -0.024205    6.151823    118 -0.036453    2.690988    119 -0.196949    0.444165    120 1.014459    5.754399    121 1.985298    3.230619    122 -1.693453    -0.557540    123 -0.576525    11.778922    024 -0.346811    -1.678730    125 -2.124484    2.672471    126 1.217916    9.597015    027 -0.733928    9.098687    028 -3.642001    -1.618087    129 0.315985    3.523953    130 1.416614    9.619232    031 -0.386323    3.989286    132 0.556921    8.294984    133 1.224863    11.587360    034 -1.347803    -2.406051    135 1.196604    4.951851    136 0.275221    9.543647    037 0.470575    9.332488    038 -1.889567    9.542662    039 -1.527893    12.150579    040 -1.185247    11.309318    041 -0.445678    3.297303    142 1.042222    6.105155    143 -0.618787    10.320986    044 1.152083    0.548467    145 0.828534    2.676045    146 -1.237728    10.549033    047 -0.683565    -2.166125    148 0.229456    5.921938    149 -0.959885    11.555336    050 0.492911    10.993324    051 0.184992    8.721488    052 -0.355715    10.325976    053 -0.397822    8.058397    054 0.824839    13.730343    055 1.507278    5.027866    156 0.099671    6.835839    157 -0.344008    10.717485    058 1.785928    7.718645    159 -0.918801    11.560217    060 -0.364009    4.747300    161 -0.841722    4.119083    162 0.490426    1.960539    163 -0.007194    9.075792    064 0.356107    12.447863    065 0.342578    12.281162    066 -0.810823    -1.466018    167 2.530777    6.476801    168 1.296683    11.607559    069 0.475487    12.040035    070 -0.783277    11.009725    071 0.074798    11.023650    072 -1.337472    0.468339    173 -0.102781    13.763651    074 -0.147324    2.874846    175 0.518389    9.887035    076 1.015399    7.571882    077 -1.658086    -0.027255    178 1.319944    2.171228    179 2.056216    5.019981    180 -0.851633    4.375691    1

 

测试数据：

 1 -1.510047    6.061992    0 2 -1.076637    -3.181888    1 3 1.821096    10.283990    0 4 3.010150    8.401766    1 5 -1.099458    1.688274    1 6 -0.834872    -1.733869    1 7 -0.846637    3.849075    1 8 1.400102    12.628781    0 9 1.752842    5.468166    110 0.078557    0.059736    111 0.089392    -0.715300    112 1.825662    12.693808    013 0.197445    9.744638    014 0.126117    0.922311    115 -0.679797    1.220530    116 0.677983    2.556666    117 0.761349    10.693862    018 -2.168791    0.143632    119 1.388610    9.341997    020 0.317029    14.739025    0

 

python 代码：

  1 #!/usr/bin/python  2 import sys  3 import copy  4 import math  5 import time  6 import random  7 import getopt  8   9 def usage(): 10     print '''Help Information: 11     -h, --help: show help information; 12     -r, --train: train file; 13     -t, --test: test file; 14     -k, --ratio: study ratio; 15     -i, --iter: iter num; 16     -p, --type: optimize type:"gradDescent","stocGradDescent","smoothStocGradDescent"; 17     ''' 18  19 def getparamenter(): 20     try: 21         opts, args = getopt.getopt(sys.argv[1:], "hr:t:k:i:p:", ["help","train=","test=","kst=","iter=","type="]) 22     except getopt.GetoptError, err: 23         print str(err) 24         usage() 25         sys.exit(1) 26  27     sys.stderr.write("\ntrain.py : a python script for perception training.\n") 28     sys.stderr.write("Copyright 2016 sxron, search, Sogou. \n") 29     sys.stderr.write("Email: [e-mail] \n\n") 30  31     train = '' 32     test = '' 33     kst = 0.01 34     iter = 100 35     type = 'gradDescent' 36     for i, f in opts: 37         if i in ("-h", "--help"): 38             usage() 39             sys.exit(1) 40         elif i in ("-r", "--train"): 41             train = f 42         elif i in ("-t", "--test"): 43             test = f 44         elif i in ("-k", "--ratio"): 45             kst = float(f) 46         elif i in ("-i", "--iter"): 47             iter = int(f) 48         elif i in ("-p", "--type"): 49             type = f 50         else: 51             assert False, "unknown option" 52    53     print "start trian parameter\ttrain:%s\ttest:%s\tkst:%f\titer:%d\ttype:%s" % (train,test,kst,iter,type) 54  55     return train,test,kst,iter,type 56  57 def loadData(file): 58     data = [] 59     label = [] 60     fin = open(file,'r') 61     while 1: 62         line = fin.readline() 63         if not line: 64             break 65         tokens = line.strip().split('\t') 66         fea = [] 67         try: 68             lab = float(tokens[-1]) 69             fea.append(1.0) 70             for i in range(0,len(tokens)-1,1): 71                 value = float(tokens[i]) 72                 fea.append(value) 73         except: 74             continue 75         label.append(lab) 76         data.append(fea) 77     return data,label 78  79 def sigmoid(inX):   80     return 1.0/(1+math.exp(-inX)) 81  82 def getMatResult(data,weights): 83     result = 0.0 84     for i in range(0,len(data),1): 85         result += data[i]*weights[i] 86     return result 87  88 def trainLogRegress(data,label,iter,kst,type): 89     weights = [] 90     for i in range(0,len(data[0]),1): 91         weights.append(1.0) 92  93     for i in range(0,iter,1): 94         errors = [] 95         if type=="gradDescent": 96             for k in range(0,len(data),1): 97                 result = getMatResult(data[k],weights) 98                 error = label[k] - sigmoid(result) 99                 errors.append(error)100             for k in range(0,len(weights),1):101                  updata = 0.0102                  for idx in range(0,len(errors),1):103                      updata += errors[idx]*data[idx][k]104                  weights[k] += kst*updata105 106         elif type=="stocGradDescent":107             for k in range(0,len(data),1):108                 result = getMatResult(data[k],weights)109                 error = label[k] - sigmoid(result)110                 for idx in range(0,len(weights),1):111                     weights[idx] += kst*error*data[k][idx]112 113         elif type=="smoothStocGradDescent":114             dataIndex = range(len(data))115             for k in range(0,len(data),1):116                 randIndex = int(random.uniform(0,len(dataIndex)))117                 result = getMatResult(data[randIndex],weights)118                 error = label[randIndex] - sigmoid(result)119                 for idx in range(0,len(weights),1):120                     weights[idx] += kst*error*data[randIndex][idx]121         else:122             print "Not support optimize method type!"123     return weights124 125 def testLogRegress(weights,data,label):126     testNum = 0127     matchNum = 0128     for i in range(0,len(data),1):129         result = getMatResult(data[i],weights)130         predict = 0131         if sigmoid(result)>0.5:132             predict = 1133         testNum += 1134         if predict==int(label[i]):135             matchNum += 1136     print "testNum:%d\tmatchNum:%d\tratio:%f" % (testNum,matchNum,float(matchNum)/testNum)137 138 def main():139     #set parameter140     train,test,kst,iter,type = getparamenter()141 142     #load train data143     trnData,trnLabel = loadData(train)144     testData,testLabel = loadData(test)145 146     #train logregress147     weights = trainLogRegress(trnData,trnLabel,iter,kst,type)148 149     #test logregress150     testLogRegress(weights,testData,testLabel)151 152 if __name__=="__main__":153     main(）

下面的程序取自：http://www.tuicool.com/articles/vMzaEvj

以下为一元回归线性方式，其中y是因变量，X是自变量，我们只需求出截距b0和斜率b1就可以获得费用和曝光次数之间的关系，并对曝光次数进行预测。这里我们使用最小二乘法来计算截距b0和斜率b1。最小二乘法通过最小化误差的平方和寻找数据的最佳函数匹配。

下表中是使用最小二乘法计算回归方程的一些必要的计算过程。在表中最左侧的两列分别为自变量X和因变量Y，我们首先计算出自变量和因变量的均值，然后计算每一个观测值与均值的差，以及用于计算回归方程斜率b1所需的数据。

根据表中的数据按公式计算出了回归方程的斜率b1，计算过程如下。斜率表示了自变量和因变量间的关系，斜率为正表示自变量和因变量正相关，斜率为负表示自变量和因变量负相关，斜率为0表示自变量和因变量不相关。

求得斜率b1后，按下面的公式可以求出Y轴的截距b0。

将斜率b1和截距b0代入到回归方程中，通过这个方程我们可以获得自变量和因变量的关系，费用每增加1元，曝光次数会增长7437次。以下为回归方程和图示。

在回归方程的图示中，还有一个R平方，这个值叫做判定系数，用来衡量回归方程是否很好的拟合了样本的数据。判定系数在0-1之间，值越大说明拟合的越好，换句话说就是自变量对因变量的解释度越高。判定系数的计算公式为SST=SSR+SSE，其中SST是总平方和，SSR是回归平方和，SSE是误差平方和。下表为计算判定系数所需三个指标的一些必要的计算过程。

根据前面求得的回归平方和(SSR)和总平方和(SST)求得判定系数为0.94344。

以上为回归方程的计算过程，在根据费用预测曝光数量的场景下，我们可以通过回归方程在已知费用的情况下计算出曝光数量。逻辑回归与回归方程相比在线性回归的基础上增加了一个逻辑函数。例如通过用户的属性和特征来判断用户最终是否会进行购买。其中购买的概率是因变量Y，用户的属性和特征是自变量X。Y值越大说明用户购买的概率越大。这里我们使用事件发生的可能性（odds）来表示购买与未购买的比值。

使用E作为购买事件，P(E)是购买的概率，P(E’)是未购买的概率，Odds(E)是事件E(购买)发生的可能性。

Odds是一个从0到无穷的数字，Odds的值越大，表明事件发生的可能性越大。下面我们要将Odds转化为0-1之间的概率函数。首先对Odds取自然对数，得到logit方程，logit是一个范围在负无穷到正无穷的值。

基于上面的logit方程，获得以下公式：

其中使用π替换了公式中的P(E),π=P(E)。根据指数函数和对数规则获得以下公式：

并最终获得逻辑回归方程：

下面根据逻辑回归方程来计算用户购买的概率，下表是用户注册天数和是否购买的数据，其中注册天数是自变量X，是否购买是自变量Y。我们将购买标记为1，将未购买标记为0。接下来我们将在Excel中通过8个步骤计算出逻辑回归方程的斜率和截距。并通过方程预测新用户是否会购买。

• 第一步，使用Excel的排序功能对原始数据按因变量Y进行排序，将已购买和未购买的数据分开，使得数据特征更加明显。
• 第二步，按照Logit方程预设斜率b1和截距b0的值，这里我们将两个值都预设为0.1。后续再通过Excel求最优解。
• 第三步，按照logit方程，使用之前预设的斜率和截距值计算出L值。

• 第四步，将L值取自然对数，
• 第五步，计算P(X)的值，P(X)为事件发生的可能性(Odds)。具体的计算步骤和过程见下图。

• 第六步，计算每个值的对数似然函数估计值（Log-Likelihood）。方法和过程见下图。
• 第七步，将对数似然函数值进行汇总。

• 第八步，使用Excel的规划求解功能，计算最大对数似然函数值。方法和过程见下图。设置汇总的对数似然函数值LL为最大化的目标，预设的斜率b1和截距b0是可变单元格，取消”使无约束变量为非负数”的选项。进行求解。

Excel将自动求出逻辑回归方程中斜率和截距的最优解，结果如下图所示。

求得逻辑回归方程的斜率和截距以后，我们可以将值代入方程，获得一个注册天数与购买概率的预测模型，通过这个模型我们可以对不同注册天数(X)用户的购买概率(Y)进行预测。以下为计算过程。

• 第一步，输入自变量注册天数(X)的值，这里我们输入50天。
• 第二步，将输入的X值，以及斜率和截距套入Logit方程，求出L值。
• 第三步，对L值取自然对数。
• 第四步，求时间发生可能性P(X)的概率值。

注册天数为50天的用户购买的概率约为17.60%。

我们将所有注册天数的值代入到购买概率预测模型中，获得了一条注册天数对购买概率影响的曲线。从曲线中可以发现，注册天数在较低和较高天数的用户购买概率较为平稳。中间天数用户的购买概率变化较大。

我们继续在上面的计算结果中增加新的自变量“年龄”。以下是原始数据的截图。现在有年龄和注册天数两个自变量和一个因变量。

依照前面的方法计算斜率和截距的最优解，并获得逻辑回归方程，将不同的年龄和注册天数代入到方程中，获得了用户年龄和注册天数对购买的预测模型。我们通过Excel的三维图表来绘制年龄和注册天数对购买概率的影响。

从图中可以看出，购买概率随着注册天数的增加而增长，并且在相同的注册天数下，年龄较小的用户购买概率相对较高。