吳恩達機器學習作業2- 邏輯回歸與正則化作業（python實現）

機器學習練習2 python復現- 邏輯回歸

在此練習中，需要實現邏輯回歸應用於分類任務。還通過將正則化加入訓練算法中來提高算法的魯棒性，並用更復雜的情形進行測試。

邏輯回歸

在訓練的初始階段，我們將要構建一個邏輯回歸模型來預測，某個學生是否被大學錄取。設想你是大學相關部分的管理者，想通過申請學生兩次測試的評分，來決定他們是否被錄取。現在你擁有之前申請學生的可以用於訓練邏輯回歸的訓練樣本集。對於每一個訓練樣本，你有他們兩次測試的評分和最后是被錄取的結果。為了完成這個預測任務，我們准備構建一個可以基於兩次測試評分來評估錄取可能性的分類模型。

數據檢查

#導包
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

#數據導入與查看
path = 'ex2data1.txt'
data = pd.read_csv(path,header = None,names = ['Exam 1','Exam 2','Admitted'])
data.head()

	Exam 1	Exam 2	Admitted
0	34.623660	78.024693	0
1	30.286711	43.894998	0
2	35.847409	72.902198	0
3	60.182599	86.308552	1
4	79.032736	75.344376	1

#創建兩次測試分數對應散點圖，並使用顏色編碼進行可視化，樣本為正表示被錄取，負表示未被錄取
#取出被錄取的兩次測試數據
positive = data[data['Admitted'].isin([1])]
#取出未被錄取的兩次測試數據
negative = data[data['Admitted'].isin([0])]
#繪制散點圖
fig,ax = plt.subplots(figsize = (12,8))
ax.scatter(positive['Exam 1'],positive['Exam 2'],s = 50,c = 'b',marker = 'o',label = 'Admitted')
ax.scatter(negative['Exam 1'],negative['Exam 2'],s = 50,c = 'r',marker = 'x',label = 'Not Admitted')
ax.legend()
ax.set_xlabel('Exam 1 Score')
ax.set_ylabel('Exam 2 Score')
plt.show()

​

​

觀察可以發現兩類間有清晰的決策邊界，現在我們需要實現邏輯回歸，就可以訓練一個模型來預測結果。

sigmoid函數

g代表一個常用的邏輯函數為S形函數，公式為：

合起來，我們得到邏輯回歸模型的假設函數：

def sigmoid(z):
    return 1 / (1+np.exp(-z))

#sigmoid函數的簡單運用
nums = np.arange(-10,10,step = 1)
fig,ax = plt.subplots(figsize = (12,8))
ax.plot(nums,sigmoid(nums),'r')
plt.show()

​

​

編寫代價函數評估結果。

def cost(theta,X,y):
    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)))
    return np.sum(first-second)/ (len(X))

#調整原始數據矩陣
data.insert(0,'Ones',1)

cols = data.shape[1]
X = data.iloc[:,0:cols - 1]
y = data.iloc[:,cols-1:cols]

X = np.array(X.values)
y = np.array(y.values)
theta = np.zeros(3)

theta

array([0., 0., 0.])

X.shape,theta.shape,y.shape

((100, 3), (3,), (100, 1))

計算初始化參數的代價函數(theta為0)

cost(theta,X,y)

0.6931471805599453

梯度下降(gradient descent)

def gradient(theta,X,y):
    theta = np.matrix(theta)
    X = np.matrix(X)
    y = np.matrix(y)
    
    parameters = int(theta.ravel().shape[1])
    grad = np.zeros(parameters)
    
    error = sigmoid(X*theta.T) - y
    
    for i in range(parameters):
        term = np.multiply(error,X[:,i])
        grad[i] = np.sum(term) / len(X)
        
    return grad

注意，我們實際上沒有在這個函數中執行梯度下降，我們僅僅在計算一個梯度步長。在練習中，一個稱為“fminunc”的Octave函數是用來優化函數來計算成本和梯度參數。由於我們使用Python，我們可以用SciPy的“optimize”命名空間來做同樣的事情。

使用數據和初始參數為0的梯度下降的結果。

gradient(theta,X,y)

array([ -0.1       , -12.00921659, -11.26284221])

使用SciPy's truncated newton（TNC）實現尋找最優參數。

import scipy.optimize as opt

result = opt.fmin_tnc(func = cost,x0=theta,fprime=gradient,args=(X,y))
result

(array([-25.16131872,   0.20623159,   0.20147149]), 36, 0)

代價函數計算的結果

cost(result[0],X,y)

0.20349770158947425

接下來，我們需要編寫一個函數，用我們所學的參數theta來為數據集X輸出預測。然后，我們可以使用這個函數來給我們的分類器的訓練精度打分。 邏輯回歸模型的假設函數：

當 ℎ𝜃 大於等於0.5時，預測 y=1

當 ℎ𝜃 小於0.5時，預測 y=0 。

def predict(theta,X):
    probability = sigmoid(X*theta.T)
    return [1 if x>= 0.5 else 0 for x in probability]

theta_min = np.matrix(result[0])
predictions = predict(theta_min,X)
correct = [1 if ((a==1 and b == 1) or (a==0 and b == 0)) else 0 for (a, b) in zip(predictions, y)]
accuracy = (sum(correct) % len(correct))
print('accuracy = {0}%'.format(accuracy))

accuracy = 89%

我們的邏輯回歸分類器預測正確，如果一個學生被錄取或沒有錄取，達到89%的精確度。不壞！記住，這是訓練集的准確性。我們沒有保持住了設置或使用交叉驗證得到的真實逼近，所以這個數字有可能高於其真實值（這個話題將在以后說明）。

正則化邏輯回歸

加入正則項提升邏輯回歸算法。

設想你是工廠的生產主管，你有一些芯片在兩次測試中的測試結果。對於這兩次測試，你想決定是否芯片要被接受或拋棄。為了幫助你做出艱難的決定，你擁有過去芯片的測試數據集，從其中你可以構建一個邏輯回歸模型。

path = 'ex2data2.txt'
data2 = pd.read_csv(path,header=None,names=['Test 1', 'Test 2', 'Accepted'])
data2.head()

	Test 1	Test 2	Accepted
0	0.051267	0.69956	1
1	-0.092742	0.68494	1
2	-0.213710	0.69225	1
3	-0.375000	0.50219	1
4	-0.513250	0.46564	1

positive = data2[data2['Accepted'].isin([1])]
negative = data2[data2['Accepted'].isin([0])]

fig, ax = plt.subplots(figsize=(12,8))
ax.scatter(positive['Test 1'], positive['Test 2'], s=50, c='b', marker='o', label='Accepted')
ax.scatter(negative['Test 1'], negative['Test 2'], s=50, c='r', marker='x', label='Rejected')
ax.legend()
ax.set_xlabel('Test 1 Score')
ax.set_ylabel('Test 2 Score')
plt.show()

​

​

注意到其中沒有線性決策界限，來良好的分開兩類數據。一個方法是用像邏輯回歸這樣的線性技術來構造從原始特征的多項式中得到的特征。讓我們通過創建一組多項式特征入手吧。

degree = 5
x1 = data2['Test 1']
x2 = data2['Test 2']

data2.insert(3, 'Ones', 1)

for i in range(1, degree):
    for j in range(0, i):
        data2['F' + str(i) + str(j)] = np.power(x1, i-j) * np.power(x2, j)

data2.drop('Test 1', axis=1, inplace=True)
data2.drop('Test 2', axis=1, inplace=True)

data2.head()

	Accepted	Ones	F10	F20	F21	F30	F31	F32	F40	F41	F42	F43
0	1	1	0.051267	0.002628	0.035864	0.000135	0.001839	0.025089	0.000007	0.000094	0.001286	0.017551
1	1	1	-0.092742	0.008601	-0.063523	-0.000798	0.005891	-0.043509	0.000074	-0.000546	0.004035	-0.029801
2	1	1	-0.213710	0.045672	-0.147941	-0.009761	0.031616	-0.102412	0.002086	-0.006757	0.021886	-0.070895
3	1	1	-0.375000	0.140625	-0.188321	-0.052734	0.070620	-0.094573	0.019775	-0.026483	0.035465	-0.047494
4	1	1	-0.513250	0.263426	-0.238990	-0.135203	0.122661	-0.111283	0.069393	-0.062956	0.057116	-0.051818

此時需要修改第1部分的成本和梯度函數，包括正則化項，首先是代價函數

def costReg(theta,X,y,learningRate):
    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)))
    reg = (learningRate / (2 * len(X))) * np.sum(np.power(theta[:,1:theta.shape[1]], 2))
    return np.sum(first - second) / len(X) + reg

添加正則化梯度函數：

如果我們需要梯度下降法令這個代價函數最小化，因為我們未對𝜃0進行正則化，所以梯度下降分為兩種情形：

整理可得：

def gradientReg(theta, X, y, learningRate):
    theta = np.matrix(theta)
    X = np.matrix(X)
    y = np.matrix(y)
    
    parameters = int(theta.ravel().shape[1])
    grad = np.zeros(parameters)
    
    error = sigmoid(X * theta.T) - y
    
    for i in range(parameters):
        term = np.multiply(error, X[:,i])
        
        if (i == 0):
            grad[i] = np.sum(term) / len(X)
        else:
            grad[i] = (np.sum(term) / len(X)) + ((learningRate / len(X)) * theta[:,i])
    
    return grad

初始化變量

# set X and y (remember from above that we moved the label to column 0)
cols = data2.shape[1]
X2 = data2.iloc[:,1:cols]
y2 = data2.iloc[:,0:1]

# convert to numpy arrays and initalize the parameter array theta
X2 = np.array(X2.values)
y2 = np.array(y2.values)
theta2 = np.zeros(11)

learningRate = 1

costReg(theta2, X2, y2, learningRate)

0.6931471805599454

gradientReg(theta2, X2, y2, learningRate)

array([0.00847458, 0.01878809, 0.05034464, 0.01150133, 0.01835599,
       0.00732393, 0.00819244, 0.03934862, 0.00223924, 0.01286005,
       0.00309594])

使用優化函數計算優化結果

result2 = opt.fmin_tnc(func=costReg, x0=theta2, fprime=gradientReg, args=(X2, y2, learningRate))
result2

(array([ 0.53010248,  0.29075567, -1.60725764, -0.58213819,  0.01781027,
        -0.21329508, -0.40024142, -1.37144139,  0.02264304, -0.9503358 ,
         0.0344085 ]),
 22,
 1)

最后使用第1部分中的預測函數查看我們方案在訓練數據上的准確度

theta_min = np.matrix(result2[0])
predictions = predict(theta_min, X2)
correct = [1 if ((a == 1 and b == 1) or (a == 0 and b == 0)) else 0 for (a, b) in zip(predictions, y2)]
accuracy = (sum(map(int, correct)) % len(correct))
print ('accuracy = {0}%'.format(accuracy))

accuracy = 78%