邏輯迴歸模型介紹和程式實現

阿新 • • 發佈：2019-02-16

　　雖然叫做“迴歸”，但是這個演算法是用來解決分類問題的。迴歸與分類的區別在於：迴歸所預測的目標量的取值是連續的（例如房屋的價格）；而分類所預測的目標變數的取值是離散的（例如判斷郵件是否為垃圾郵件）。當然，為了便於理解，我們從二值分類（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}為分類結果。 clip_image006

對於訓練資料集，特徵資料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    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: [email protected] \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(）

以下為一元迴歸線性方式，其中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）來表示購買與未購買的比值。