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

逻辑回归原理及实现

 　　虽然叫做“回归”，但是这个算法是用来解决分类问题的。回归与分类的区别在于：回归所预测的目标量的取值是连续的（例如房屋的价格）；而分类所预测的目标变量的取值是离散的（例如判断邮件是否为垃圾邮件）。当然，为了便于理解，我们从二值分类（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={x1, x2, … , xm}和对应的分类数据y={y1, y2, … , ym}。构建逻辑回归模型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    0
10 -0.026632    10.427743    0
11 0.850433    6.920334    1
12 1.347183    13.175500    0
13 1.176813    3.167020    1
14 -1.781871    9.097953    0
15 -0.566606    5.749003    1
16 0.931635    1.589505    1
17 -0.024205    6.151823    1
18 -0.036453    2.690988    1
19 -0.196949    0.444165    1
20 1.014459    5.754399    1
21 1.985298    3.230619    1
22 -1.693453    -0.557540    1
23 -0.576525    11.778922    0
24 -0.346811    -1.678730    1
25 -2.124484    2.672471    1
26 1.217916    9.597015    0
27 -0.733928    9.098687    0
28 -3.642001    -1.618087    1
29 0.315985    3.523953    1
30 1.416614    9.619232    0
31 -0.386323    3.989286    1
32 0.556921    8.294984    1
33 1.224863    11.587360    0
34 -1.347803    -2.406051    1
35 1.196604    4.951851    1
36 0.275221    9.543647    0
37 0.470575    9.332488    0
38 -1.889567    9.542662    0
39 -1.527893    12.150579    0
40 -1.185247    11.309318    0
41 -0.445678    3.297303    1
42 1.042222    6.105155    1
43 -0.618787    10.320986    0
44 1.152083    0.548467    1
45 0.828534    2.676045    1
46 -1.237728    10.549033    0
47 -0.683565    -2.166125    1
48 0.229456    5.921938    1
49 -0.959885    11.555336    0
50 0.492911    10.993324    0
51 0.184992    8.721488    0
52 -0.355715    10.325976    0
53 -0.397822    8.058397    0
54 0.824839    13.730343    0
55 1.507278    5.027866    1
56 0.099671    6.835839    1
57 -0.344008    10.717485    0
58 1.785928    7.718645    1
59 -0.918801    11.560217    0
60 -0.364009    4.747300    1
61 -0.841722    4.119083    1
62 0.490426    1.960539    1
63 -0.007194    9.075792    0
64 0.356107    12.447863    0
65 0.342578    12.281162    0
66 -0.810823    -1.466018    1
67 2.530777    6.476801    1
68 1.296683    11.607559    0
69 0.475487    12.040035    0
70 -0.783277    11.009725    0
71 0.074798    11.023650    0
72 -1.337472    0.468339    1
73 -0.102781    13.763651    0
74 -0.147324    2.874846    1
75 0.518389    9.887035    0
76 1.015399    7.571882    0
77 -1.658086    -0.027255    1
78 1.319944    2.171228    1
79 2.056216    5.019981    1
80 -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    1
10 0.078557    0.059736    1
11 0.089392    -0.715300    1
12 1.825662    12.693808    0
13 0.197445    9.744638    0
14 0.126117    0.922311    1
15 -0.679797    1.220530    1
16 0.677983    2.556666    1
17 0.761349    10.693862    0
18 -2.168791    0.143632    1
19 1.388610    9.341997    0
20 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: shixiang08abc@gmail.com \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.0
102                  for idx in range(0,len(errors),1):
103                      updata += errors[idx]*data[idx][k]
104                  weights[k] += kst*updata
105
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 weights
124
125 def testLogRegress(weights,data,label):
126     testNum = 0
127     matchNum = 0
128     for i in range(0,len(data),1):
129         result = getMatResult(data[i],weights)
130         predict = 0
131         if sigmoid(result)>0.5:
132             predict = 1
133         testNum += 1
134         if predict==int(label[i]):
135             matchNum += 1
136     print "testNum:%d\tmatchNum:%d\tratio:%f" % (testNum,matchNum,float(matchNum)/testNum)
137
138 def main():
139     #set parameter
140     train,test,kst,iter,type = getparamenter()
141
142     #load train data
143     trnData,trnLabel = loadData(train)
144     testData,testLabel = loadData(test)
145
146     #train logregress
147     weights = trainLogRegress(trnData,trnLabel,iter,kst,type)
148
149     #test logregress
150     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的三维图表来绘制年龄和注册天数对购买概率的影响。



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