周志華《機器學習》課後習題解答系列（四）：Ch3.3 - 程式設計實現對率迴歸

這裡採用Python-sklearn的方式，環境搭建可參考 資料探勘入門：Python開發環境搭建（eclipse-pydev模式）.

相關答案和原始碼託管在我的Github上：PY131/Machine-Learning_ZhouZhihua.

思路概要

程式設計實現對率迴歸：
* 採用sklearn邏輯斯蒂迴歸庫函式實現，通過檢視混淆矩陣，繪製決策區域來檢視模型分類效果；
* 自己程式設計實現，從極大化似然函式出發，採用梯度下降法得到最優引數，然後嘗試了隨機梯度下降法來優化過程。

3.3 程式設計實現對率迴歸

所使用的資料集如下：

本題是本書的第一個程式設計練習，採用了自己程式設計實現和呼叫sklearn庫函式兩種不同的方式，詳細解答和編碼過程：（檢視完整程式碼）：

1.獲取資料、檢視資料、預處理：

觀察資料可知，X包含（密度、含糖量）兩個變數，y為西瓜是否好瓜分類（二分），由此生成.csv資料檔案，在Python中用Numpy讀取資料並採用matplotlib庫視覺化資料：

樣例程式碼：

'''
data importion
'''
import numpy as np
import matplotlib.pyplot as plt  

# load the CSV file as a numpy matrix
dataset = np.loadtxt('../data/watermelon_3a.csv', delimiter=",")

# separate the data from the target attributes
X = dataset[:,1:3]
y = dataset[:,3]

# draw scatter diagram to show the raw data
f1 = plt.figure(1)   
plt.title('watermelon_3a')  
plt.xlabel('density')  
plt.ylabel('ratio_sugar')  
plt.scatter(X[y == 0,0], X[y == 0,1], marker = 'o', color = 'k', s=100, label = 'bad')
plt.scatter(X[y == 1,0], X[y == 1,1], marker = 'o', color = 'g', s=100, label = 'good')
plt.legend(loc = 'upper right')  
plt.show()

資料散點圖：

2.採用sklearn邏輯迴歸庫函式直接擬合：

雖然樣本量很少，這裡還是先劃分訓練集和測試集，採用sklearn.model_selection.train_test_split()實現，然後採用sklearn.linear_model.LogisticRegression，基於訓練集直接擬合出邏輯迴歸模型，然後在測試集上評估模型（檢視混淆矩陣和F1值）。

樣例程式碼：

''' 
using sklearn lib for logistic regression
'''
from sklearn import model_selection
from sklearn.linear_model import LogisticRegression
from sklearn import metrics


# generalization of test and train set
X_train, X_test, y_train, y_test = model_selection.train_test_split(X, y, test_size=0.5, random_state=0)

# model training
log_model = LogisticRegression() 
log_model.fit(X_train, y_train) 

# model testing
y_pred = log_model.predict(X_test)

# summarize the accuracy of fitting
print(metrics.confusion_matrix(y_test, y_pred))
print(metrics.classification_report(y_test, y_pred))

得出混淆矩陣和相關度量（查準率（準確率）、查全率（召回率），F1值）結果如下：

[[4 1]
 [1 3]]
             precision    recall  f1-score   support

        0.0       0.80      0.80      0.80         5
        1.0       0.75      0.75      0.75         4

avg / total       0.78      0.78      0.78         9

由混淆矩陣可以看到，由於樣本本身數量較少，模型擬合效果一般，總體預測精度約為0.78。為提升精度，可以採用自助法進行重抽樣擴充資料集，或是採用交叉驗證選擇最優模型。

下圖是採用matplotlib.contourf繪製的決策區域和邊界，可以看出對率迴歸分類器還是成功的分出了絕大多數類：

3.自己程式設計實現邏輯斯蒂迴歸

程式設計實現邏輯迴歸的主要工作是求取引數w和b（見書p59），最常用的引數估計方法是極大似然法，由於題3.1已經證得對數似然函式（見書3.27）是凸函式，存在最優解，這裡考慮採用梯度下降法來迭代尋優。

回顧一下Sigmoid函式，即邏輯斯蒂迴歸分類器的基礎模型：

目的是基於資料集求出最優引數w和b，最常採用的是極大似然法，引數的似然函式為：

根據書p59，最大化上式等價於最小化下式：

題3.2已證上式為凸函式，一定存在最小值，但按照導數為零的解析求解方式較為困難，於是考慮採用梯度下降法來求解上式最小值時對應的引數。

注：梯度下降法基本知識可參考書中附錄p409頁，也可直接採用書中p60式3.30偏導數公式。書中關於引數迭代改變式子如下：

對於迭代，可每次先根據(B.16)計算出梯度▽f(β)，然後由(B.17)更新得出下一步的Δβ。

接下來程式設計實現基本的梯度下降法：

(1)首先程式設計實現物件式3.27：

def likelihood_sub(x, y, beta):
    '''
    @param x: one sample variables
    @param y: one sample label
    @param beta: the parameter vector in 3.27
    @return: the sub_log-likelihood of 3.27  
    ''' 
    return -y * np.dot(beta, x.T) + np.math.log(1 + np.math.exp(np.dot(beta, x.T)))   

def likelihood(X, y, beta):
    '''
    @param X: the sample variables matrix
    @param y: the sample label matrix
    @param beta: the parameter vector in 3.27
    @return: the log-likelihood of 3.27  
    '''
    sum = 0
    m,n = np.shape(X)  

    for i in range(m):
        sum += likelihood_sub(X[i], y[i], beta)

    return sum  

(2)然後基於訓練集（注意x->[x,1]），給出基於3.27似然函式的定步長梯度下降法，注意這裡的偏梯度實現技巧：

'''
@param X: X is the variable matrix 
@param y: y is the label array
@return: the best parameter estimate of 3.27
'''
def gradDscent_1(X, y):  #implementation of basic gradDscent algorithms

    h = 0.1  # step length of iteration
    max_times= 500  # give the iterative times limit
    m, n = np.shape(X)

    beta = np.zeros(n)  # parameter and initial to 0
    delta_beta = np.ones(n)*h  # delta parameter and initial to h
    llh = 0
    llh_temp = 0

    for i in range(max_times):
        beta_temp = beta.copy()

        # for partial derivative 
        for j in range(n): 
            beta[j] += delta_beta[j]
            llh_tmp = likelihood(X, y, beta)
            delta_beta[j] = -h * (llh_tmp - llh) / delta_beta[j]                
            beta[j] = beta_temp[j]

        beta += delta_beta            
        llh = likelihood(X, y, beta)

    return beta

通過追蹤引數，檢視其收斂曲線，然後來調節相關引數（步長h，迭代次數max_times）。下圖是在當前引數取值下的beta曲線，可以看到其收斂良好：

(3)最後建立Sigmoid預測函式，對測試集資料進預測，得到混淆矩陣如下：

[[ 4.  1.]
 [ 1.  3.]]

可以看出其總體預測精度（7/9 ≈ 0.78）與呼叫sklearn庫得出的結果相當。

(4)採用隨機梯度下降法來優化：上面採用的是全域性定步長梯度下降法（稱之為批量梯度下降），這種方法在可能會面臨收斂過慢和收斂曲線波動情況的同時，每次迭代需要全域性計算，計算量隨資料量增大而急劇增大。所以嘗試採用隨機梯度下降來改善引數迭代尋優過程。

隨機梯度下降法的核心思想是增量學習：一次只用一個新樣本來更新迴歸係數，從而形成線上流式處理。

同時為了加快收斂，採用變步長的策略，h隨著迭代次數逐漸減小。

給出變步長隨機梯度下降法的程式碼如下：

def gradDscent_2(X, y):  #implementation of stochastic gradDscent algorithms
  '''
   @param X: X is the variable matrix 
   @param y: y is the label array
   @return: the best parameter estimate of 3.27
   '''
   import matplotlib.pyplot as plt  

   m, n = np.shape(X)
   h = 0.5  #  step length of iterator and initial
   beta = np.zeros(n)  # parameter and initial
   delta_beta = np.ones(n) * h
   llh = 0
   llh_temp = 0

   for i in range(m):
       beta_temp = beta.copy()  # for partial derivative 

       for j in range(n): 
           h = 0.5 * 1 / (1 + i + j)  # change step length of iterator 
           beta[j] += delta_beta[j]
           llh_tmp = likelihood_sub(X[i], y[i], beta)
           delta_beta[j] = -h * (llh_tmp - llh) / delta_beta[j]   
           beta[j] = beta_temp[j]  

       beta += delta_beta    
       llh = likelihood_sub(X[i], y[i], beta)

   return beta

得出混淆矩陣：

[[ 3.  2.]
 [ 0.  4.]]

從結果看到的是：由於這裡的西瓜資料集並不大，所以隨機梯度下降法採用一次遍歷所得的結果不太好，引數也沒有完成收斂。這裡只是給出隨機梯度下降法的實現樣例，這種方法在大資料集下相比批量梯度法應會有明顯的優勢。

---

參考連結：
由於這是本書第一個程式設計，索引資料較多，擇其重要的一些列出如下：

• Introduction to Machine Learning with Python and Scikit-Learn
• scikit-learn官方主頁
• matplotlib官方主頁
• 隨機梯度下降（Stochastic gradient descent）和批量梯度下降（Batch gradient descent）的公式對比、實現對比
• 機器學習演算法與Python實踐之（七）邏輯迴歸（Logistic Regression）
• plot decision boundary matplotlib
  Ask（matplotlib決策區域和邊界繪製
• Python資料視覺化——散點圖