注：最近開始學習《人工智能》選修課，老師提綱挈領的介紹了一番，聽完課只了解了個大概，剩下的細節只能自己繼續摸索。
從本質上講：機器學習就是一個模型對外界的刺激（訓練樣本）做出反應，趨利避害（評價標准）。
1. 什么是邏輯回歸？
許多人對線性回歸都比較熟悉，但知道邏輯回歸的人可能就要少的多。從大的類別上來說，邏輯回歸是一種有監督的統計學習方法，主要用於對樣本進行分類。
在線性回歸模型中，輸出一般是連續的，例如$$y = f(x) = ax + b$$，對於每一個輸入的x，都有一個對應的y輸出。模型的定義域和值域都可以是[-∞, +∞]。但是對於邏輯回歸，輸入可以是連續的[-∞, +∞]，但輸出一般是離散的，即只有有限多個輸出值。例如，其值域可以只有兩個值{0, 1}，這兩個值可以表示對樣本的某種分類，高/低、患病/健康、陰性/陽性等，這就是最常見的二分類邏輯回歸。因此，從整體上來說，通過邏輯回歸模型，我們將在整個實數范圍上的x映射到了有限個點上，這樣就實現了對x的分類。因為每次拿過來一個x，經過邏輯回歸分析，就可以將它歸入某一類y中。
邏輯回歸與線性回歸的關系
邏輯回歸也被稱為廣義線性回歸模型，它與線性回歸模型的形式基本上相同，都具有 ax+b，其中a和b是待求參數，其區別在於他們的因變量不同，多重線性回歸直接將ax+b作為因變量，即y = ax+b，而logistic回歸則通過函數S將ax+b對應到一個隱狀態p，p = S(ax+b)，然后根據p與1-p的大小決定因變量的值。這里的函數S就是Sigmoid函數
$$S(t) = \frac{1}{1 + e^{-t}}$$
將t換成ax+b，可以得到邏輯回歸模型的參數形式：$$p(x; a,b) = \frac{1}{1 + e^{-(ax+b)}} ……（1）$$
圖1：sigmoid函數的圖像
通過函數S的作用，我們可以將輸出的值限制在區間[0， 1]上，p(x)則可以用來表示概率p(y=1|x)，即當一個x發生時，y被分到1那一組的概率。可是，等等，我們上面說y只有兩種取值，但是這里卻出現了一個區間[0, 1]，這是什么鬼？？其實在真實情況下，我們最終得到的y的值是在[0, 1]這個區間上的一個數，然后我們可以選擇一個閾值，通常是0.5，當y>0.5時，就將這個x歸到1這一類，如果y<0.5就將x歸到0這一類。但是閾值是可以調整的，比如說一個比較保守的人，可能將閾值設為0.9，也就是說有超過90%的把握，才相信這個x屬於1這一類。了解一個算法，最好的辦法就是自己從頭實現一次。下面是邏輯回歸的具體實現。
邏輯回歸模型的代價函數
邏輯回歸一般使用交叉熵作為代價函數。關於代價函數的具體細節，請參考代價函數，這里只給出交叉熵公式：
$$J(\theta) = -\frac{ 1 }{ m }[\sum_{ i=1 }^{ m } ({y^{(i)} \log h_\theta(x^{(i)}) + (1-y^{(i)}) \log (1-h_\theta(x^{(i)})})]$$
m：訓練樣本的個數；
hθ(x)：用參數θ和x預測出來的y值；
y：原訓練樣本中的y值，也就是標准答案
上角標(i)：第i個樣本
2. 數據准備
下面的數據來自《機器學習實戰》中的示例：
-0.017612 14.053064 0 -1.395634 4.662541 1 -0.752157 6.538620 0 -1.322371 7.152853 0 0.423363 11.054677 0 0.406704 7.067335 1 0.667394 12.741452 0 -2.460150 6.866805 1 0.569411 9.548755 0 -0.026632 10.427743 0
上面的數據一共是3列10行，其中前兩列為x1和x2的值，第3列表示y的值；10行表示取了10個樣本點。我們可以將這些數據當做訓練模型參數的訓練樣本。
見到訓練樣本就可以比較直觀的理解算法的輸入，以及我們如何利用這些數據來訓練邏輯回歸分類器，進而用訓練好的模型來預測新的樣本（檢測樣本）。
從邏輯回歸的參數形式，式子（1）我們可以看到邏輯回歸模型中有兩個待定參數a（x的系數）和b（常數項），我們現在給出來的數據有兩個特征x1, x2，因此整個模型就增加了一項：ax1 + cx2 + b。為了形式上的統一，我們使用帶下標的a表示不同的參數（a0表示常數項b並作x0的參數<x0=1>，a1、a2分別表示x1和x2的參數），就可以得到：
$$ a_0x_0 + a_1x_1 + a_2x_2 $$
這樣統一起來后，就可以使用矩陣表示了（比起前面展開的線性表示方式，用矩陣表示模型和參數更加簡便，而且矩陣運算的速度也更快）：
$$ \begin{bmatrix} a_0 & a_1 & a_2 \end{bmatrix} \begin{bmatrix} x_0 \\ x_1 \\ x_2 \end{bmatrix} = a^{ \mathrm{ T } }X$$
將上面的式子帶入到（1）式，我們就可以得到邏輯回歸的另一種表示形式了：
$$p(x; a) = \frac{1}{1 + e^{-a^{ \mathrm{ T } }X}} ……（2）$$
此時，可以很清楚的看到，我們后面的行動都是為了確定一個合適的a（一個參數向量），使得對於一個新來的X（也是一個向量），我們可以盡可能准確的給出一個y值，0或者1.
注：數據是二維的，也就是說這組觀察樣本中有兩個自變量，即兩個特征（feature）。
3. 訓練分類器
就像上面說的，訓練分類器的過程，就是根據已經知道的數據（訓練樣本）確定一個使得代價函數的值最小的a（參數向量/回歸系數）的過程。邏輯回歸模型屬於有監督的學習方法，上面示例數據中的第3列其實是訓練樣本提供的"標准答案"。也就是說，這些數據是已經分好類的（兩類，0或者1）。在訓練階段，我們要做的就是利用訓練樣本和（2）式中的模型，估計一個比較合適的參數a，使得僅通過前面兩列數據（觀察值/測量值）就可以估計一個值h(a)，這個值越接近標准答案y，說明我們的模型預測的越准確。
下面是估計回歸系數a的值的過程，還是借鑒了《機器學習實戰》中的代碼，做了少量修改：
其中計算參數梯度，即代價函數對每個參數的偏導數（下面代碼中的第36-38行），的詳細推導過程可以參考這里
1 ''' 2 Created on Oct 27, 2010 3 Logistic Regression Working Module 4 @author: Peter 5 ''' 6 from numpy import * 7 import os 8 9 path = 'D:\MechineLearning\MLiA_SourceCode\machinelearninginaction\Ch05' 10 training_sample = 'trainingSample.txt' 11 testing_sample = 'testingSample.txt' 12 13 # 從文件中讀入訓練樣本的數據，同上面給出的示例數據 14 # 下面第20行代碼中的1.0表示x0 = 1 15 def loadDataSet(p, file_n): 16 dataMat = []; labelMat = [] 17 fr = open(os.path.join(p, file_n)) 18 for line in fr.readlines(): 19 lineArr = line.strip().split() 20 dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])]) # 三個特征x0, x1, x2 21 labelMat.append(int(lineArr[2])) # 標准答案y 22 return dataMat,labelMat 23 24 def sigmoid(inX): 25 return 1.0/(1+exp(-inX)) 26 27 # 梯度下降法求回歸系數a，由於樣本量少，我將迭代次數改成了1000次 28 def gradAscent(dataMatIn, classLabels): 29 dataMatrix = mat(dataMatIn) #convert to NumPy matrix 30 labelMat = mat(classLabels).transpose() #convert to NumPy matrix 31 m,n = shape(dataMatrix) 32 alpha = 0.001 # 學習率 33 maxCycles = 1000 34 weights = ones((n,1)) 35 for k in range(maxCycles): # heavy on matrix operations 36 h = sigmoid(dataMatrix*weights) # 模型預測值, 90 x 1 37 error = h - labelMat # 真實值與預測值之間的誤差, 90 x 1 38 temp = dataMatrix.transpose()* error # 交叉熵代價函數對所有參數的偏導數, 3 x 1 39 weights = weights - alpha * temp # 更新權重 40 return weights 41 42 # 下面是我自己寫的測試函數 43 def test_logistic_regression(): 44 dataArr, labelMat = loadDataSet(path, training_sample) # 讀入訓練樣本中的原始數據 45 A = gradAscent(dataArr, labelMat) # 回歸系數a的值 46 h = sigmoid(mat(dataArr)*A) #預測結果h(a)的值 47 print(dataArr, labelMat) 48 print(A) 49 print(h) 50 # plotBestFit(A) 51 52 test_logistic_regression()
上面代碼的輸出如下：
- 一個元組，包含兩個數組：第一個數組是所有的訓練樣本中的觀察值，也就是X，包括x0, x1, x2；第二個數組是每組觀察值對應的標准答案y。
([[1.0, -0.017612, 14.053064], [1.0, -1.395634, 4.662541], [1.0, -0.752157, 6.53862], [1.0, -1.322371, 7.152853], [1.0, 0.423363, 11.054677], [1.0, 0.406704, 7.067335], [1.0, 0.667394, 12.741452], [1.0, -2.46015, 6.866805], [1.0, 0.569411, 9.548755], [1.0, -0.026632, 10.427743]], [0, 1, 0, 0, 0, 1, 0, 1, 0, 0])
- 本次預測出來的回歸系數a，包括a0, a1, a2
[[ 1.39174871] [-0.5227482 ] [-0.33100373]]
- 根據回歸系數a和（2）式中的模型預測出來的h(a)。這里預測得到的結果都是區間(0, 1)上的實數。
[[ 0.03730313] [ 0.64060602] [ 0.40627881] [ 0.4293251 ] [ 0.07665396] [ 0.23863652] [ 0.0401329 ] [ 0.59985228] [ 0.11238742] [ 0.11446212]]
標准答案是{0, 1}，如何將預測到的結果與標准答案y進行比較呢？取0.5作為閾值，大於該值的樣本就划分到1這一組，小於等於該值的樣本就划分到0這一組，這樣就可以將數據分為兩類。檢查一下結果可以看到，我們現在分出來的1這一類中包括原來y=1的兩個樣本，另一類包括原來y=0的所有樣本和一個y=1的樣本（分錯了）。鑒於我們選擇取的樣本比較少（只有10個），這樣的效果其實還算非常不錯的！
4. 結果展示
上面已經求出了一組回歸系數，它確定了不同類別數據之間的分割線。可以利用X內部（x1與x2之間的關系）的關系畫出該分割線，從而更直觀的感受到分類的效果。
添加下面一段代碼：
1 # 分類效果展示，參數weights就是回歸系數 2 def plotBestFit(weights): 3 import matplotlib.pyplot as plt 4 dataMat,labelMat=loadDataSet(path, training_sample) 5 dataArr = array(dataMat) 6 n = shape(dataArr)[0] 7 xcord1 = []; ycord1 = [] 8 xcord2 = []; ycord2 = [] 9 for i in range(n): 10 if int(labelMat[i])== 1: 11 xcord1.append(dataArr[i,1]); ycord1.append(dataArr[i,2]) 12 else: 13 xcord2.append(dataArr[i,1]); ycord2.append(dataArr[i,2]) 14 fig = plt.figure() 15 ax = fig.add_subplot(111) 16 ax.scatter(xcord1, ycord1, s=30, c='red', marker='s') 17 ax.scatter(xcord2, ycord2, s=30, c='green') 18 x = arange(-3.0, 3.0, 0.1) 19 y = (-weights[0]-weights[1]*x)/weights[2] # x2 = f(x1) 20 ax.plot(x.reshape(1, -1), y.reshape(1, -1)) 21 plt.xlabel('X1'); plt.ylabel('X2'); 22 plt.show()
將上面的test_logistic_regression()函數中的最后一句注釋去掉，調用plotBestFit函數就可以看到分類的效果了。
這里說明一下上面代碼中的第19行，這里設置了sigmoid函數的取值為1/2，也就是說取閾值為0.5來划分最后預測的結果。這樣可以得到$$e^{-a^{ \mathrm{ T } }X} = 1$$，即$-a^TX=0$，可以推出$x_2 = (-a_0x_0 - a_1x_1)/a_2$，同第19行，也就是說這里的$y$實際上是$x_1$，而$x$是$x_1$。因此下圖表示的是$x_1$與$x_2$之間的關系。
分類效果圖如下：
三個紅色的點是原來$y=1$的樣本，有一個分錯了。這里相當於將所有的數據用二維坐標(x1, x2)表示了出來，而且根據回歸參數畫出的線將這些點一分為二。如果有新的樣本，不知道在哪一類，只用將該點畫在圖上，看它在這條直線的哪一邊就可以分類了。
下面是使用90個訓練樣本得到的結果：
可以看出一個非常明顯的規律是，$y=1$的這一類樣本（紅色的點）具有更小的$x_2$值，當$x_2$相近時則具有更大的$x_1$值。
此時計算出來的回歸系數a為：
[[ 5.262118 ] [ 0.60847797] [-0.75168429]]
5. 預測新樣本
添加一個預測函數，如下：
直接將上面計算出來的回歸系數a拿來使用，測試數據其實也是《機器學習實戰》這本書中的訓練數據，我拆成了兩份，前面90行用來做訓練數據，后面10行用來當測試數據。
1 def predict_test_sample(): 2 A = [5.262118, 0.60847797, -0.75168429] # 上面計算出來的回歸系數a 3 dataArr, labelMat = loadDataSet(path, testing_sample) 4 h_test = sigmoid(mat(dataArr) * mat(A).transpose()) # 將讀入的數據和A轉化成numpy中的矩陣 5 print(h_test) # 預測的結果
調用上面的函數，可以得到以下結果，即h(a)：
[[ 0.99714035] [ 0.04035907] [ 0.12535895] [ 0.99048731] [ 0.98075409] [ 0.97708653] [ 0.09004989] [ 0.97884487] [ 0.28594188] [ 0.00359693]]
下面是我們的測試數據（原來的訓練樣本后十行的數據，包括標准答案y）：
0.089392 -0.715300 1 1.825662 12.693808 0 0.197445 9.744638 0 0.126117 0.922311 1 -0.679797 1.220530 1 0.677983 2.556666 1 0.761349 10.693862 0 -2.168791 0.143632 1 1.388610 9.341997 0 0.317029 14.739025 0
比較我們預測得到的h(a)和標准答案y，如果按照0.5為分界線的話，我們利用前90個樣本訓練出來的分類器對后面10個樣本的類型預測全部正確。
附件：
github上的代碼更新到python3.6, 2019-1-6
完整代碼：https://github.com/OnlyBelter/MachineLearning_examples/tree/master/de_novo/regression
訓練數據：https://github.com/OnlyBelter/MachineLearning_examples/blob/master/de_novo/data/Logistic_Regression-trainingSample.txt
測試數據：https://github.com/OnlyBelter/MachineLearning_examples/blob/master/de_novo/data/Logistic_Regression-testingSample.txt
參考：
http://baike.baidu.com/item/logistic%E5%9B%9E%E5%BD%92
https://en.wikipedia.org/wiki/Sigmoid_function
《機器學習實戰》，哈林頓著，李銳等譯，人民郵電出版社，2013年6月第一版