機器學習 - 線性回歸與邏輯回歸（實踐部分）

之前對線性回歸和邏輯回歸的理論部分做了較為詳細的論述，下面通過一些例子再來鞏固一下之前所學的內容。
需要說明的是，雖然我們在線性回歸中都是直接通過公式推導求出w和b的精確值，但在實際運用中基本上都會采用梯度下降法作為首選，因為用代碼表示公式會比較繁瑣，而梯度下降法只需要不斷對參數更新公式進行迭代即可，用代碼表示非常簡單，且迭代次數基本上不會對性能有太大的拖累，所以下面我們也將采用梯度下降的方法。

使用線性回歸模擬數據分布

了解數據初始分布情況

現有一個數據集（提取碼：yim9），描述了一系列點的分布情況。我們先將初始的分布情況繪制出來，便於對數據有個整體的了解：

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression

data = np.loadtxt('./data/linear_regression_data.txt', delimiter=',')
print(data[:10])  # 查看部分樣本數據
data.shape        # 查看樣本規模

# 將數據集轉換成 y = w*x 的向量形式
x = np.c_[np.ones(data.shape[0]), data[:, 0]]  # 使用np.c_轉換為列向量
y = np.c_[data[:, 1]]  # 1是給b留的位置

# 展示樣本分布情況
plt.scatter(x[:, 1], y, color='green', marker='x', linewidth=1)
plt.xlabel('X')
plt.ylabel('Y')

可以看到，點在X軸左側分布得比較密集，但是越往右，越呈現出按照直線分布的趨勢，所以可以嘗試使用線性回歸進行模擬。

使用最小二乘法構造目標函數

我們需要求出目標函數，以便求出梯度，並記錄目標函數在每輪迭代中的值

# 目標函數
def target_func(x, y, w=[[0],[0]]):
    # 這里我們設置w和b的初始值為0
    h = x.dot(w)                                   # h表示預測值
    m = y.size                                     # 通常會使用樣本大小m平均整體誤差的值
    target = 1.0/(2*m) * (np.sum(np.square(h-y)))  # 使用最小二乘法得到目標函數
    return target

計算梯度及參數更新公式

# 計算梯度
def grad(x, y, w):
    h = x.dot(w)
    m = y.size
    grad = (1.0/m)*(x.T.dot(h-y))  # 對目標函數求導即是梯度
    return grad

# 使用梯度循環更新參數
def grad_descent(x, y, nums_iter=1000, n=0.01, w=[[0],[0]]):
    # 默認迭代1000次，學習率為0.01
    target_his = np.zeros(nums_iter)    # 存放目標函數每輪迭代后的值
    for i in np.arange(nums_iter):
        # 循環更新w，並記錄目標函數的值
        w = w - n*grad(x, y, w)
        target_his[i] = target_func(x, y, w)
    return w, target_his

w, target_his = grad_descent(x, y)

# 展示迭代效果
plt.plot(target_his)
plt.xlabel('Iter Times')
plt.ylabel('Loss Function')

可以看出，迭代1000次，性能逐漸趨於平穩。
最后與sklearn的結果相對比：

# 樣本
plt.scatter(x[:, 1], y, color='red', marker='x', linewidth=1)
# 自主
xx = np.arange(5, 23)
yy = w[0] + xx*w[1]
plt.plot(xx, yy, label='LD(GR)')

# sklearn
regr = LinearRegression()
# reshape轉換成列向量
regr.fit(x[:, 1].reshape(-1, 1), y.ravel())
yy2 = regr.intercept_ + regr.coef_ * xx
plt.plot(xx, yy2, label='LD(sklearn)')

plt.xlim(4, 24)
plt.xlabel('X')
plt.ylabel('Y')
plt.legend(loc=4)

可以看出效果基本一致。

使用邏輯回歸分類

現有一個考試成績數據集(提取碼：7n04)，該數據集為某班級英語和數學兩個科目的考試成績，前兩列分別表示英語和數學成績，最后一列表示最終是否通過。

加載數據，查看數據樣本及可視化

import pandas as pd
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from scipy.optimize import minimize

# Step 1. 加載數據
def loaddata(data):
     return np.loadtxt(data, delimiter=',')

data = loaddata('./data/sample_1.txt')
# 查看數據
print(data.shape)
print(data[:5, :])

# 運行結果：
(100, 3)
[[34.62365962 78.02469282  0.        ]
 [30.28671077 43.89499752  0.        ]
 [35.84740877 72.90219803  0.        ]
 [60.18259939 86.3085521   1.        ]
 [79.03273605 75.34437644  1.        ]]

# Step 2. 整理訓練樣本
# np.c_是按行連接兩個矩陣，就是把兩矩陣左右相加，要求行數相等
X = np.c_[np.ones((data.shape[0],1)), data[:,0:2]]
y = np.c_[data[:,2]]
print(X[:5])
print(X.shape)
print(y[:5])

# 運行結果：
[[ 1.         34.62365962 78.02469282]
 [ 1.         30.28671077 43.89499752]
 [ 1.         35.84740877 72.90219803]
 [ 1.         60.18259939 86.3085521 ]
 [ 1.         79.03273605 75.34437644]]
(100, 3)
[[0.]
 [0.]
 [0.]
 [1.]
 [1.]]

X的第一列是給偏置b預留的，其余兩列為特征值。

# Step 3. 可視化
def plotData(data, label_x, label_y, label_pos, label_neg, axes=None):
    # 獲得正負樣本的下標(正樣本表示考試通過的，負樣本表示考試未通過)
    neg = data[:,2] == 0
    pos = data[:,2] == 1
    if axes == None:
        axes = plt.gca() # 獲取圖的坐標信息

    # 依次繪制正負樣本，使用不同顏色區分
    axes.scatter(data[pos][:,0], data[pos][:,1], marker='*', c='r', s=30, linewidth=2, label=label_pos)
    axes.scatter(data[neg][:,0], data[neg][:,1], c='c', s=30, label=label_neg)
    axes.set_xlabel(label_x)
    axes.set_ylabel(label_y)
    axes.legend(frameon= True, fancybox = True);
    
plotData(data, 'English', 'Math', 'Pass', 'Fail')

從樣本分布可以看出，正負樣本分布的還算整齊，可以使用決策邊界進行划分。

定義Sigmoid函數、損失函數（即目標函數）、梯度函數

# Step 4. 訓練前准備
# Sigmoid函數
def sigmoid(z):
    return(1 / (1 + np.exp(-z)))

# 損失函數
def lossFunction(w, X, y):
    m = y.size
    h = sigmoid(X.dot(w.reshape(-1, 1)))
    # 套公式得到損失函數
    J = -1.0*(1.0/m)*(np.log(h).T.dot(y)+np.log(1-h).T.dot(1-y))
    if np.isnan(J[0]):
        return(np.inf)
    return J[0]

# 梯度計算函數
def gradient(X, y, w):
    grad = np.zeros(w.shape)
    error = (sigmoid(X.dot(w.reshape(-1, 1))).reshape(-1, 1) - y).ravel()
    for j in range(len(w.ravel())):
        # 循環更新w，得到梯度
        term = np.multiply(error, X[:,j])
        grad[j] = np.sum(term) / len(X)
    return grad

# 初始化w
initial_w = np.zeros(X.shape[1])
loss = lossFunction(initial_w, X, y)
grad = gradient(X, y, initial_w)
print('Loss: \n', loss)
print('Grad: \n', grad)

# 輸出結果：
Loss: 
 [0.69314718]
Grad: 
 [ -0.1        -12.00921659 -11.26284221]

關於損失函數，代碼里使用的是下面這個公式，而不是之前推導到最后的公式，因為用代碼更好實現：

使用梯度下降法訓練

# step 5. 使用梯度下降法訓練
# 梯度下降函數
def gradientDescent(X, y, w, alpha=0.01, num_iters=1000):
    m = y.size
    J_history = np.zeros(num_iters)
    for iter in np.arange(num_iters):
        h = X.dot(w.reshape(-1, 1))
        grad = gradient(X, y, w)
        w = w - alpha*(1.0/m)*grad
        J_history[iter] = lossFunction(w, X, y)
    return(w, J_history)

# 畫出每一次迭代和損失函數變化
w , Cost_J = gradientDescent(X, y, initial_w, alpha=0.15, num_iters=100000)

plt.plot(Cost_J)
plt.ylabel('Loss J')
plt.xlabel('Iter')

通過將學習率設置為0.15，並經過10W次的迭代后，損失函數計算出來的誤差控制在了0.32左右。
下面對迭代出來的w進行可視化。

可視化決策邊界

# Step 6. 可視化
plotData(data, 'English', 'Math', 'Pass', 'Fail')
# np.linspace(x1_min, x1_max) 生成的是x軸英語成績的列向量
# np.linspace(x2_min, x2_max) 生成的是y軸數學成績的列向量
# xx1, xx2為meshgrid生成的坐標矩陣
x1_min, x1_max = X[:,1].min(), X[:,1].max(),
x2_min, x2_max = X[:,2].min(), X[:,2].max(),
xx1, xx2 = np.meshgrid(np.linspace(x1_min, x1_max), np.linspace(x2_min, x2_max))
# 將預測值作為等高線繪制出來，即決策邊界
# 將xx1和xx2作為x的特征值，得到預測值h
h = sigmoid(np.c_[np.ones((xx1.ravel().shape[0],1)), xx1.ravel(), xx2.ravel()].dot(w.T))
h = h.reshape(xx1.shape)
plt.contour(xx1, xx2, h, 1, linewidths=1, colors='b');

使用scipy.optimize.minimize實現

scipy.optimize.minimize可以根據現有的損失函數和梯度，優化整個梯度下降的過程，可以自動選擇學習率，從而免去了不斷手動調試的麻煩。

from scipy.optimize import minimize

res = minimize(lossFunction, initial_w, args=(X,y), jac=gradient, options={'maxiter':400})
print(res)

傳入相應的參數，並設置迭代輪數（這里僅設置了400次），就可得到優化后的w。從輸出可以看到，損失函數的誤差已經減小到了0.2，比之前手動實現的要小很多，而且更加省時。

# Step 5. 可視化
plotData(data, 'English', 'Math', 'Pass', 'Fail')
x1_min, x1_max = X[:,1].min(), X[:,1].max(),
x2_min, x2_max = X[:,2].min(), X[:,2].max(),
xx1, xx2 = np.meshgrid(np.linspace(x1_min, x1_max), np.linspace(x2_min, x2_max))

h = sigmoid(np.c_[np.ones((xx1.ravel().shape[0],1)), xx1.ravel(), xx2.ravel()].dot(res.x.T))
h = h.reshape(xx1.shape)
plt.contour(xx1, xx2, h, 1, linewidths=1, colors='b');

從訓練結果來看，手動和調包實現的效果好像差不多，但是明顯后者更加省時省力。