Python實現邏輯斯蒂回歸

Python實現邏輯斯蒂回歸

邏輯回歸與加正則化的邏輯斯蒂回歸介紹

本實驗室根據兩次考試成績與是否通過的數據，通過logistic回歸，最後獲得一個分類器。

邏輯斯蒂回歸

導入數據

import numpy as npdef loaddata(file, delimeter): #以delimeter為分隔符導入file數據 data = np.loadtxt(file, delimiter=delimeter) print(Dimensions: ,data.shape) # 列印數據前6行 print(data[1:6,:]) return(data)data = loaddata(data1.txt, ,)

結果為：

Dimensions: (100, 3)

[[30.28671077 43.89499752 0. ]

[35.84740877 72.90219803 0. ]

[60.18259939 86.3085521 1. ]

[79.03273605 75.34437644 1. ]

[45.08327748 56.31637178 0. ]]

可以看見data的數據結構為一個100*3的矩陣，第一列為exam1的成績，第二列為exam2的成績，第三列為是否最終通過（0為否，1為是）。

# 作圖顯示數據分布import matplotlib.pyplot as pltdef 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=k, s=60, linewidth=2, label=label_pos) axes.scatter(data[neg][:,0], data[neg][:,1], c=y, s=60, label=label_neg) axes.set_xlabel(label_x) axes.set_ylabel(label_y) axes.legend(frameon= True, fancybox = True)plotData(data, Exam 1 score, Exam 2 score, Pass, Fail)

讀取數據作為X與y向量：

X = np.c_[np.ones((data.shape[0],1)), data[:,0:2]]y = np.c_[data[:,2]]

則X為：

1 34.6237 78.0247

1 30.2867 43.895

1 35.8474 72.9022

1 60.1826 86.3086

1 79.0327 75.3444

1 45.0833 56.3164

1 61.1067 96.5114

1 75.0247 46.554

1 76.0988 87.4206

1 84.4328 43.5334......

y為：

0

0

0

1

1

0

1

1

1

1......

邏輯斯蒂回歸

邏輯斯蒂回歸函數設為為：

# 定義sigmoid函數def sigmoid(z): return(1 / (1 + np.exp(-z)))

損失函數為：

向量化的損失函數（矩陣形式）：

# 定義損失函數def costFunction(theta, X, y): m = y.size # 此處h為一個列向量 h = sigmoid(X.dot(theta)) J = -1.0*(1.0/m)*(np.log(h).T.dot(y)+np.log(1-h).T.dot(1-y)) return J[0]# 初始化theta為0向量initial_theta = np.zeros(X.shape[1])# 此時邏輯斯蒂回歸的損失函數值為cost = costFunction(initial_theta, X, y)print(Cost:
, cost)

Cost:

0.6931471805599452

求偏導（梯度）

向量化的偏導（梯度）

# 求梯度def gradient(theta, X, y): m = y.size # 此處h為一100*1的列向量 h = sigmoid(X.dot(theta.reshape(-1,1))) # grad為一個3*100的矩陣 grad =(1.0/m)*X.T.dot(h-y) return(grad.flatten())#此時邏輯斯蒂函數在初始化theta處在Θ0、Θ1、Θ2處的梯度分別為：grad = gradient(initial_theta, X, y)print(Grad:
, grad)

Grad:

[ -0.1 -12.00921659 -11.26284221]

最小化損失函數

minimize()函數的介紹

Jacobian矩陣和Hessian矩陣介紹

# 計算損失函數得最小值。from scipy.optimize import minimize# 規定最大迭代次數為400次res = minimize(costFunction, initial_theta, args=(X,y), jac=gradient, options={maxiter:400})res# minimize()返回的格式是固定的，fun為costFunction函數迭代求得的最小值，hess_inv和jac分別為求得最小值時海森矩陣和雅克比矩陣的值，小寫字母x為當costFunction函數最小時函數的解，在costFunction中即為theta的解。即計算損失函數最小時theta的值為[-25.16133284, 0.2062317 , 0.2014716 ]。Out[17]: fun: 0.20349770158944375 hess_inv: array([[ 3.31474479e+03, -2.63892205e+01, -2.70237122e+01], [ -2.63892205e+01, 2.23869433e-01, 2.02682332e-01], [ -2.70237122e+01, 2.02682332e-01, 2.35335117e-01]]) jac: array([ -9.52476821e-09, -9.31921318e-07, -2.82608930e-07]) message: Optimization terminated successfully. nfev: 31 nit: 23 njev: 31 status: 0 success: True x: array([-25.16133284, 0.2062317 , 0.2014716 ])

即損失函數最小時theta的值：

theta = res.x.TthetaOut[20]: array([-25.16133284, 0.2062317 , 0.2014716 ])

咱們來看看考試1得分45，考試2得分85的同學通過概率有多高

sigmoid(np.array([1,45,85]).dot(theta))Out[22]: 0.77629072405889421

即考試1得分45，考試2得分85的同學通過概率約為0.7763。

畫出決策邊界

# 標註考試1得分45，考試2得分85的同學plt.scatter(45, 85, s=60, c=r, marker=v, label=(45, 85))# plotData為之前定義的畫分類點的函數plotData(data, Exam 1 score, Exam 2 score, Admitted, Not admitted)# 生成網格數據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的值h = sigmoid(np.c_[np.ones((xx1.ravel().shape[0],1)), xx1.ravel(), xx2.ravel()].dot(res.x))h = h.reshape(xx1.shape)# 作等高線，可理解為在這條線上h值為0.5plt.contour(xx1, xx2, h, [0.5], linewidths=1, colors=b)

加正則化項的邏輯斯蒂回歸

導入數據

data2 = loaddata(data2.txt, ,)

data2的格式為：

Dimensions: (118, 3)

[[-0.092742 0.68494 1. ]

[-0.21371 0.69225 1. ]

[-0.375 0.50219 1. ]

[-0.51325 0.46564 1. ]

[-0.52477 0.2098 1. ]]

作分布圖

y = np.c_[data2[:,2]]X = data2[:,0:2]plotData(data2, Microchip Test 1, Microchip Test 2, y = 1, y = 0)

在上一個邏輯回歸試驗中，我們把sigmoid函數（即這裡的g函數）設置的為簡單的一次多項式。

在這個邏輯回歸實驗里，因為樣本的分布比較複雜，可以採用多次多項式來代替ΘTX。這裡取最高六次項。

取高階多項式放入sigmoid函數進行模擬

from sklearn.preprocessing import PolynomialFeatures# 生成一個六次多項式poly = PolynomialFeatures(6)# XX為生成的六次項的數據XX = poly.fit_transform(data2[:,0:2])# 六次項後有28個特徵值了。即，之前我們只有兩個特徵值x1、x2，取六次項多項式後我們會有x1、x2、x1^2、x2^2、x1*x2、x1^2*x2、……，總共28項。XX.shapeOut[12]: (118, 28)

正則化

因為取得的多項式最高項為6次，容易發生過擬合情況。將損失函數採取「正則化」處理，引入懲罰項。

正則化後損失函數：

向量化的損失函數（矩陣形式）：

# 定義損失函數def costFunctionReg(theta,reg ,XX, y): m = y.size h = sigmoid(XX.dot(theta)) J = -1.0*(1.0/m)*(np.log(h).T.dot(y)+np.log(1-h).T.dot(1-y)) +(reg/(2.0*m))*np.sum(np.square(theta[1:])) if np.isnan(J[0]): return(np.inf) return(J[0])

與之對應的偏導（梯度）：

向量化的偏導（梯度）：

注意，我們另外自己加的參數 θ0 不需要被正則化

# 定義正則化損失函數的偏導def gradientReg(theta, reg, XX, y): m = y.size h = sigmoid(XX.dot(theta.reshape(-1,1))) grad = (1.0/m)*XX.T.dot(h-y) + (reg/m)*np.r_[[[0]],theta[1:].reshape(-1,1)] return(grad.flatten())

設初始化theta為0向量，計算此時初始損失值

initial_theta = np.zeros(XX.shape[1])costFunctionReg(initial_theta, 1, XX, y)Out[9]: 0.69314718055994529

畫出決策邊界

定義預測函數，用來統計準確率。分類的閾值定為0.5，即計算的h(x)>0.5則分到1類（即通過），h(x)<0.5則分到0類（即不通過）：

def predict(theta, X, threshold=0.5): h = sigmoid(X.dot(theta.T)) >= threshold # 返回的h值只會有兩種，1或0 return(h.astype(int))

決策邊界，咱們分別來看看正則化係數lambda太大太小分別會出現什麼情況

• Lambda = 0 : 就是沒有正則化，這樣的話，就過擬合咯
• Lambda = 1 : 這才是正確的打開方式
• Lambda = 100 : 卧槽，正則化項太激進，導致基本就沒擬合出決策邊界

fig, axes = plt.subplots(1,3, sharey = True, figsize=(17,5))# 分別取lambda為0、1、100for i, C in enumerate([0.0, 1.0, 100.0]): # 最優化 costFunctionReg res2 = minimize(costFunctionReg, initial_theta, args=(C, XX, y), jac=gradientReg, options={maxiter:3000}) # 準確率 accuracy = 100.0*sum(predict(res2.x, XX) == y.ravel())/y.size # 對X,y的散列繪圖 plotData(data2, Microchip Test 1, Microchip Test 2, y = 1, y = 0, axes.flatten()[i]) # 畫出決策邊界 x1_min, x1_max = X[:,0].min(), X[:,0].max(), x2_min, x2_max = X[:,1].min(), X[:,1].max(), xx1, xx2 = np.meshgrid(np.linspace(x1_min, x1_max), np.linspace(x2_min, x2_max)) h = sigmoid(poly.fit_transform(np.c_[xx1.ravel(), xx2.ravel()]).dot(res2.x)) h = h.reshape(xx1.shape) axes.flatten()[i].contour(xx1, xx2, h, [0.5], linewidths=1, colors=g); axes.flatten()[i].set_title(Train accuracy {}% with Lambda = {}.format(np.round(accuracy, decimals=2), C))

推薦閱讀：