機器學習練習二：用Python實現邏輯迴歸

時間 2019-12-02

原文原文鏈接

1.邏輯迴歸

數據介紹：在本部分練習中，您將構建一個邏輯迴歸模型來預測學生是否被大學錄取。假設您是一個大學部門的管理員，您但願根據兩次考試的結果來肯定每一個申請人的錄取機會。您可使用之前申請者的歷史數據做爲邏輯迴歸的訓練集。程序員

1.1 畫出原始數據圖

import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
# 展現原始數據
data = pd.read_csv('ex2data1.txt',names=['score1','score2','admitted'])
# 布爾值索引data['Admitted'].isin([1])：False,True,False...
positive = data[data['admitted'].isin([1])]
negative = data[data['admitted'].isin([0])]
# 子畫布
fig, ax = plt.subplots(figsize=(12,8))
ax.scatter(positive['score1'], positive['score2'], s=50, c='g', marker='o', label='Admitted')
ax.scatter(negative['score1'], negative['score2'], s=50, c='r', marker='x', label='Not Admitted')
ax.legend()  # 標籤
ax.set_xlabel('Score1')
ax.set_ylabel('Score2')
plt.show()
複製代碼

其中data['admitted'].isin([1])可以將admitted列中知足值爲1的位置變爲True，不知足的位置變爲False，即生成這樣的Series: 算法

而後經過布爾值取data中的數據，生成兩個Dataframe：positive、negative。

1.2 數據處理

插入全1列的目的是匹配偏置 $\theta_0$ 數組

# 插入全1列x0
data.insert(0, 'Ones', 1)
# set X (training data) and y (target variable)
cols = data.shape[1]
X = data.iloc[:,0:cols-1]  # X是全部行，去掉最後一列
y = data.iloc[:,cols-1:cols]  # X是全部行，最後一列
theta = np.array([0,0,0])
複製代碼

1.3 sigmoid函數

邏輯迴歸須要使用sigmoid函數來當作目標函數，其中 $h_\theta(x)$ 的值表示樣本爲1類別的機率。

def sigmoid(x):
    return 1 / (1 + np.exp(-x))
複製代碼

# 畫出sigmoid
nums = np.arange(-10,10,step=1)
plt.figure(figsize=(20,8),dpi=100)
plt.plot(nums,sigmoid(nums))
plt.show()
複製代碼

1.4 代價函數和梯度降低

代價函數： $J\left( \theta \right)=\frac{1}{m}\sum\limits_{i=1}^{m}{[-{{y}^{(i)}}\log \left( {{h}_{\theta }}\left( {{x}^{(i)}} \right) \right)-\left( 1-{{y}^{(i)}} \right)\log \left( 1-{{h}_{\theta }}\left( {{x}^{(i)}} \right) \right)]}$ 代碼實現：bash

def cost(theta, X, y):
    # 將參數轉換爲矩陣類型
    theta = np.matrix(theta)
    X = np.matrix(X)
    y = np.matrix(y)
    # X(100,3),y(100,1),theta(1,3)
    first = np.multiply(-y, np.log(sigmoid(X * theta.T)))
    second = np.multiply((1 - y), np.log(1 - sigmoid(X * theta.T)))
    return np.sum(first - second) / (len(X))
複製代碼

梯度降低：函數

代碼實現：

這裏梯度降低函數只實現了 $\frac{\partial J\left( \theta \right)}{\partial {{\theta }_{j}}}=\frac{1}{m}\sum\limits_{i=1}^{m}{({{h}_{\theta }}\left( {{x}^{(i)}} \right)-{{y}^{(i)}})x_{_{j}}^{(i)}}$ 部分且爲向量化實現。測試

# 構造梯度降低函數(只作了J(θ)偏導數部分)
def gradientDescent(theta, X, y):
    X = np.matrix(X)
    y = np.matrix(y)
    theta = np.matrix(theta)
    return (1 / len(X)) * X.T @ (sigmoid(X @ theta.T) - y)
複製代碼

1.5 進行梯度降低並計算準確率

1.5.1 梯度降低計算

這裏使用了scipy庫中的方法，關於該方法參數的意義可使用help查看。優化

# 使用 scipy.optimize.minimize 去尋找參數
import scipy.optimize as opt
res = opt.minimize(fun=cost, x0=theta, args=(X, y), method='TNC', jac=gradientDescent)
res
複製代碼

輸出以下：ui

其中fun參數是最終的損失值，x參數是最終的 $\theta$ 值。

1.5.2 準確度計算

獲取訓練集的預測結果

# 獲取訓練集預測結果
def predict(theta, X):
    theta = np.mat(theta)
    y_predict = sigmoid(X @ theta.T)
    return [1 if x >= 0.5 else 0 for x in y_predict]
複製代碼

判斷準確度

# 判斷準確度
y_pre = np.mat(predict(res.x, X))  # 預測值矩陣
y_true = np.mat(y).ravel()  # 真實值矩陣
# 矩陣進行比較返回各元素比對布爾值矩陣，列表進行比較返回整個列表的比對布爾值
accuracy = np.mean(y_pre == y_true)
print('accuracy = {}%'.format(accuracy * 100))
複製代碼

先將y_predict和y_true轉換爲一維矩陣，一維矩陣的布爾運算返回的是矩陣中各位置元素的布爾值，形如：[True,False,False...]。spa

而np.mean()求得矩陣中True元素佔總元素數量的比例，即預測準確度。3d

1.6 繪製決策邊界

當 $\theta^Tx=0$ 時 $h_\theta(x)=0.5$ ，即全部機率爲類別1的樣本的集合，亦即決策邊界。

顯然本次邏輯迴歸的決策邊界是線性的，即直線 $x_2=\frac{-\theta_0-\theta_1x_1}{\theta_2}$ 。

代碼實現：

# 繪製決策邊界
x1 = np.linspace(data.score1.min(), data.score1.max(), 100)  # 返回在指定範圍內的100個等間隔數
x2 = (- res.x[0] - res.x[1] * x1) / res.x[2]
# 布爾值索引data['Admitted'].isin([1])：False,True,False...
positive = data[data['admitted'].isin([1])]
negative = data[data['admitted'].isin([0])]
# 子畫布
fig, ax = plt.subplots(figsize=(12,8))
ax.scatter(positive['score1'], positive['score2'], s=50, c='g', marker='o', label='Admitted')
ax.scatter(negative['score1'], negative['score2'], s=50, c='r', marker='x', label='Not Admitted')
ax.plot(x1, x2, 'r', label='Prediction')
ax.legend()  # 標籤
ax.set_xlabel('Score1')
ax.set_ylabel('Score2')
plt.show()
複製代碼

其中x1爲在x座標軸間隨機取得的100個點，x2即爲知足要求的y值。

2.正則化邏輯迴歸

數據介紹：在練習的這一部分，您將實現規範的邏輯迴歸來預測來自制造工廠的微芯片是否經過質量保證(QA)。在QA過程當中，每一個微晶片都要通過各類測試以確保其正常工做。假設您是工廠的產品經理，您在兩個不一樣的測試中得到了一些微芯片的測試結果。從這兩個測試中，您但願肯定是否應該接受或拒絕微芯片。爲了幫助您作出決定，您有一個關於過去的微芯片的測試結果的數據集，您能夠從中得到數據。

2.1 數據展現

data2 = pd.read_csv('ex2data2.txt',names=['test1','test2','pass'])
# 布爾值索引data['Admitted'].isin([1])：False,True,False...
positive = data2[data2['pass'].isin([1])]
negative = data2[data2['pass'].isin([0])]
# 子畫布
fig, ax = plt.subplots(figsize=(12,8))
ax.scatter(positive['test1'], positive['test2'], s=50, c='g', marker='o', label='Pass')
ax.scatter(negative['test1'], negative['test2'], s=50, c='r', marker='x', label='Not Pass')
ax.legend()  # 標籤
ax.set_xlabel('test1')
ax.set_ylabel('test2')
plt.show()
複製代碼

如圖所示，咱們的數據集不能經過線性方式劃分爲正例和反例。所以，簡單的邏輯迴歸應用在這個數據集上的效果並很差，由於邏輯迴歸只能找到一個線性決策邊界。更好地適應數據的一種方法是從每一個數據點建立更多的特性，咱們將把特徵映射到全部的多項式項x1和x2的六次方。

2.2 特徵映射與數據分割

多項式展開(polynomial expansion)：

僞代碼：

for i in 0..power+1：
    for p in 0..i+1:
        output x^(i-p) * y^p
複製代碼

代碼實現：

def poly_expansion(x1,x2,power,dataframe):
    """ 多項式展開 """
    for i in range(power + 1):
        for j in range(i + 1):
            dataframe['F' + str(i-j) + str(j)] = np.power(x1, i-j) * np.power(x2, j)
    dataframe.drop('test1', axis=1, inplace=True)
    dataframe.drop('test2', axis=1, inplace=True)
複製代碼

映射後data2.shape爲(118,29)，即從兩個特徵拓展爲28個特徵。

數據分割

# 分割提取數據
# set X (training data) and y (target variable)
cols = data2.shape[1]
X = data2.iloc[:,1:cols]  # X是全部行，去掉第一列
y = data2.iloc[:,0:1]  # y是全部行，第一列，注意取值方法
theta = np.zeros(X.shape[1])
複製代碼

2.3 正則化代價函數

J\left( \theta \right)=\frac{1}{m}\sum\limits_{i=1}^{m}{[-{{y}^{(i)}}\log \left( {{h}_{\theta }}\left( {{x}^{(i)}} \right) \right)-\left( 1-{{y}^{(i)}} \right)\log \left( 1-{{h}_{\theta }}\left( {{x}^{(i)}} \right) \right)]}+\frac{\lambda }{2m}\sum\limits_{j=1}^{n}{\theta _{j}^{2}}

前半部分與普通邏輯迴歸一致，後半部分是 $\lambda$ 超參數和 $\theta_j$ 的求和的乘積，其目的是讓每個 $\theta_j$ 的值儘可能小，避免過擬合。

代碼實現：

def cost(theta, X, y, lambd):
    theta = np.matrix(theta)
    X = np.matrix(X)
    y = np.matrix(y)
    # 函數前半部分
    first = np.multiply(-y, np.log(sigmoid(X * theta.T)))
    second = np.multiply((1 - y), np.log(1 - sigmoid(X * theta.T)))
    front = np.sum(first - second) / (len(X))
    # 函數後半部分
    last = lambd / (2 * len(X)) * np.sum(np.power(theta, 2))
    return front + last
複製代碼

2.4 梯度降低

原理

由於咱們未對 ${{\theta }_{0}}$ 進行正則化，因此梯度降低算法將分兩種情形：

$\begin{align} & Repeat\text{ }until\text{ }convergence\text{ }\!\!\{\!\!\text{ } \\ & \text{ }{{\theta }_{0}}:={{\theta }_{0}}-a\frac{1}{m}\sum\limits_{i=1}^{m}{[{{h}_{\theta }}\left( {{x}^{(i)}} \right)-{{y}^{(i)}}]x_{_{0}}^{(i)}} \\ & \text{ }{{\theta }_{j}}:={{\theta }_{j}}-a\frac{1}{m}\sum\limits_{i=1}^{m}{[{{h}_{\theta }}\left( {{x}^{(i)}} \right)-{{y}^{(i)}}]x_{j}^{(i)}}-\frac{\lambda }{m}{{\theta }_{j}} \\ & \text{ }\!\!\}\!\!\text{ } \\ & Repeat \\ \end{align}$

對上面的算法中 j=1,2,...,n 時的更新式子進行調整可得：

${{\theta }_{j}}:={{\theta }_{j}}(1-a\frac{\lambda }{m})-a\frac{1}{m}\sum\limits_{i=1}^{m}{({{h}_{\theta }}\left( {{x}^{(i)}} \right)-{{y}^{(i)}})x_{j}^{(i)}}$

在正則化的梯度降低中，咱們對除了 $\theta_0$ 之外的參數進行優化，由於 $\theta_0$ 對特徵值的權重沒有影響（）。

仍然只對 $\frac{\partial J\left( \theta \right)}{\partial {{\theta }_{j}}}=\frac{1}{m}\sum\limits_{i=1}^{m}{({{h}_{\theta }}\left( {{x}^{(i)}} \right)-{{y}^{(i)}})x_{_{j}}^{(i)}} + \frac{\lambda}{m}\theta_j$ 這一部分進行代碼實現。

代碼實現：

# 構造梯度降低函數(只作了J(θ)偏導數部分)
def gradient(theta, X, y, lambd):
    X = np.matrix(X)
    y = np.matrix(y)
    theta = np.matrix(theta)
    first = (1 / len(X)) * (X.T @ (sigmoid(X @ theta.T) - y))
    # theta0不須要優化
    second = (lambd / len(X)) * theta[:,1:].T
    second = np.concatenate((np.array([[0]]),second),axis=0)
    return first + second
複製代碼

second = np.concatenate((np.array([[0]]),second),axis=0)這一句代碼是爲了在不優化 $\theta_0$ 的同時使得先後兩部分矩陣對齊（加了一個0行）。

進行梯度降低優化

和上一章方法同樣，再也不贅述。

# 使用 scipy.optimize.minimize 去尋找參數
import scipy.optimize as opt
res = opt.minimize(fun=cost, x0=theta, args=(X, y, 100), method='TNC', jac=gradient)
res
複製代碼

準確度計算

# 判斷準確度
y_pre = np.mat(predict(res.x, X))  # 預測值矩陣
y_true = np.mat(y).ravel()  # 真實值矩陣
# 矩陣進行比較返回各元素比對布爾值矩陣，列表進行比較返回整個列表的比對布爾值
accuracy = np.mean(y_pre == y_true)
print('accuracy = {}%'.format(accuracy * 100))
複製代碼

2.5 決策邊界

與上一章不一樣的是，如今咱們將要進行不規則決策邊界的繪製。顯然上文所使用的的繪製直線方程的方法再也不適用。

從數據圖中能夠看出，紅藍點的區分界限並非一條直線，而是一個不規則的形狀，這就是不規則決策邊界。那麼繪製不規則決策邊界的方法其實也很簡單，就是將特徵平面上的每個點都用咱們訓練出的模型判斷它屬於哪一類，而後將判斷出的分類顏色繪製出來，就獲得了上圖所示的效果，那麼不規則決策邊界天然也就出來了，這個原理相似繪製地形圖的等高線，在同一等高範圍內的點就是同一類。

引自：程序員說

代碼實現:

def plot_decision_boundary(axis,axe):
    # meshgrid函數用兩個座標軸上的點在平面上畫格，返回座標矩陣
    X0, X1 = np.meshgrid(
        # 隨機兩組數，起始值和密度由座標軸的起始值決定
        np.linspace(axis[0], axis[1], int((axis[1] - axis[0]) * 100)).reshape(-1, 1),
        np.linspace(axis[2], axis[3], int((axis[3] - axis[2]) * 100)).reshape(-1, 1),
    )
    # ravel()方法將高維數組降爲一維數組，c_[]將兩個數組以列的形式拼接起來，造成矩陣
    X_grid_matrix = np.c_[X0.ravel(), X1.ravel()]
    # 將特徵矩陣轉換爲DataFrame進行多項式展開
    X_grid_matrix = pd.DataFrame(X_grid_matrix,columns=["test1","test2"])
    x1 = X_grid_matrix.test1
    x2 = X_grid_matrix.test2
    poly_expansion(x1,x2,6,X_grid_matrix)
    # 經過訓練好的邏輯迴歸模型，預測平面上這些點的分類
    y_predict = np.array(predict(res.x, X_grid_matrix))
    # 將預測結果轉換爲與X0和X1形狀相同的矩陣
    y_predict_matrix = y_predict.reshape(X0.shape)
    
    # 設置色彩表
    from matplotlib.colors import ListedColormap
    my_colormap = ListedColormap(['#0000CD', '#40E0D0', '#FFFF00'])
    
    # 繪製等高線，而且填充等高區域的顏色
    ax.contourf(X0, X1, y_predict_matrix, cmap=my_colormap)
複製代碼

繪製邊界：

# 布爾值索引data['Admitted'].isin([1])：False,True,False...
positive = data2[data2['pass'].isin([1])]
negative = data2[data2['pass'].isin([0])]
# 子畫布
fig, ax = plt.subplots(figsize=(12,8))
plot_decision_boundary(axis=[data2.F10.min(), data2.F10.max(), data2.F01.min(), data2.F01.max()], axe=ax)
ax.scatter(positive['F10'], positive['F01'], s=50, c='g', marker='o', label='Pass')
ax.scatter(negative['F10'], negative['F01'], s=50, c='r', marker='x', label='Not Pass')
ax.legend()  # 標籤
ax.set_xlabel('test1')
ax.set_ylabel('test2')
plt.show()
複製代碼