Python機器學習的練習四：多元邏輯迴歸
在本系列的第3部分中，我們實現了簡單的和正則化的邏輯迴歸。但我們的解決方法有一個限制—它只適用於二進位制分類。在本文中，我們將在之前的練習中擴充套件我們的解決方案，以處理多級分類問題。
在語法上快速標註，為了顯示語句的輸出，我在程式碼塊中附加了一個“>”，以表明它是執行先前語句的結果。如果結果很長(超過1-2行)，那麼我就把它貼上在程式碼塊的另一個單獨的塊中。希望可以清楚的說明哪些語句是輸入，哪些是輸出。
此練習中的任務是使用邏輯迴歸來識別手寫數字(0-9)。首先載入資料集。與前面的示例不同,我們的資料檔案是MATLAB的本體格式，不能被pandas自動識別，所以把它載入在Python中需要使用SciPy utility。
import numpy as np import pandas as pd import matplotlib.pyplot as plt from scipy.ioimport loadmat %matplotlib inline data= loadmat('data/ex3data1.mat') data {'X': array([[0., 0., 0., ..., 0., 0., 0.], [0., 0., 0., ..., 0., 0., 0.], [0., 0., 0., ..., 0., 0., 0.], ..., [0., 0., 0., ..., 0., 0., 0.], [0., 0., 0., ..., 0., 0., 0.], [0., 0., 0., ..., 0., 0., 0.]]), '__globals__': [], '__header__':'MATLAB 5.0 MAT-file, Platform: GLNXA64, Created on: Sun Oct 16 13:09:09 2011', '__version__':'1.0', 'y': array([[10], [10], [10], ..., [9], [9], [9]], dtype=uint8)}
快速檢查載入到儲存器中的矩陣的形狀
data['X'].shape, data['y'].shape
> ((5000L,400L), (5000L,1L))我們已經載入了我們的資料。影象在martix X 被表現為400維的向量。這400個“特徵”是原始20×20影象中每個畫素的灰度強度。類標籤在向量y中表示影象中數字的數字類。下面的圖片給出了一些數字的例子。每個帶有白色手寫數字的灰色框代表我們資料集中400維的行。
我們的第一個任務是修改邏輯迴歸的實現以完全向量化（即沒有“for”迴圈），這是因為向量化程式碼除了簡潔扼要，還能夠利用線性代數優化，並且比迭代程式碼快得多。我們在練習二中的成本函式實現已經向量化。所以我們在這裡重複使用相同的實現。請注意，我們正在跳到最終的正則化版本。
def sigmoid(z):
return 1 / (1 + np.exp(-z))
def cost(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這個成本函式與我們在先前練習中建立的函式是相同的，如果你不確定我們如何做到這一點，在執行之前檢視以前的文章。
接下來，我們需要計算梯度的函式。我們已經在前面的練習中定義了它，我們在更新步驟中需要去掉“for”迴圈。這是可供參考的原始程式碼：
def gradient_with_loop(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 iin 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梯度函式詳細的闡述瞭如何改變一個引數，以獲得一個比之前好的答案。如果你不熟悉線性代數，這一系列運作背後的數學是很難理解的。
現在我們需要建立一個不使用任何迴圈的梯度函式的版本。在我們的新版本中，我們將去掉“for”迴圈，並使用線性代數（除了截距引數，它需要單獨計算）計算每個引數的梯度。
還要注意，我們將資料結構轉換為NumPy矩陣。這樣做是為了使程式碼看起來更像Octave，而不是陣列，這是因為矩陣自動遵循矩陣運算規則與element-wise運算。我們在下面的例子中使用矩陣。
def gradient(theta, X, y, learningRate):
theta= np.matrix(theta)
X= np.matrix(X)
y= np.matrix(y)
parameters= int(theta.ravel().shape[1])
error= sigmoid(X* theta.T)- y
grad= ((X.T* error)/ len(X)).T+ ((learningRate/ len(X))* theta)
# intercept gradient is not regularized
grad[0,0]= np.sum(np.multiply(error, X[:,0]))/ len(X)
return np.array(grad).ravel()現在我們已經定義了成本和梯度函式，接下來建立一個分類器。對於本章練習的任務，我們有10個可能的分類，由於邏輯迴歸一次只能區分兩個類別，我們需要一個方法去處理多類別場景。在這個練習中我們的任務是實現一對多的分類，其中k個不同類的標籤導致了k個分類器，每個分類器在“class i”和“not class i”之間做決定。我們將在一個函式中完成分類器的訓練，計算10個分類器的最終權重，並將權重返回作為k X（n + 1）陣列，其中n是引數數。
from scipy.optimizeimport minimize
def one_vs_all(X, y, num_labels, learning_rate):
rows= X.shape[0]
params= X.shape[1]
# k X (n + 1) array for the parameters of each of the k classifiers
all_theta= np.zeros((num_labels, params+ 1))
# insert a column of ones at the beginning for the intercept term
X= np.insert(X,0, values=np.ones(rows), axis=1)
# labels are 1-indexed instead of 0-indexed
for iin range(1, num_labels+ 1):
theta= np.zeros(params+ 1)
y_i= np.array([1 if label== ielse 0 for labelin y])
y_i= np.reshape(y_i, (rows,1))
# minimize the objective function
fmin= minimize(fun=cost, x0=theta, args=(X, y_i, learning_rate), method='TNC', jac=gradient)
all_theta[i-1,:]= fmin.x
return all_theta這裡需要注意的幾點：首先，我們為theta添加了一個額外的引數（帶有一列訓練資料）以計算截距項。其次，我們將y從類標籤轉換為每個分類器的二進位制值（要麼是I類，要麼不是I類）。最後，我們使用SciPy的較新優化API來最小化每個分類器的成本函式。API利用目標函式、初始引數集、優化方法和jacobian（梯度）函式，將優化程式找到的引數分配給引數陣列。
實現向量化程式碼的一個更具挑戰性的部分是正確地寫入所有的矩陣互動，所以通過檢視正在使用的陣列/矩陣的形狀來做一些健全性檢查是有用的，我們來看看上面的函式中使用的一些資料結構。
rows= data['X'].shape[0]
params= data['X'].shape[1]
all_theta= np.zeros((10, params+ 1))
X= np.insert(data['X'],0, values=np.ones(rows), axis=1)
theta= np.zeros(params+ 1)
y_0= np.array([1 if label== 0 else 0 for labelin data['y']])
y_0= np.reshape(y_0, (rows,1))
X.shape, y_0.shape, theta.shape, all_theta.shape
> ((5000L,401L), (5000L,1L), (401L,), (10L,401L))注意，theta是一維陣列，所以當它被轉換為計算梯度的程式碼中的矩陣時，它變成一個（1×401）矩陣。 我們還要檢查y中的類標籤，以確保它們看起來像我們期望的。
np.unique(data['y'])
> array([1, 2, 3, 4, 5, 6, 7, 8, 9,10], dtype=uint8)確保函式正常執行，並獲得一些合理的輸出。
all_theta= one_vs_all(data['X'], data['y'],10,1)
all_theta
array([[-5.79312170e+00, 0.00000000e+00, 0.00000000e+00, ...,
1.22140973e-02, 2.88611969e-07, 0.00000000e+00],
[-4.91685285e+00, 0.00000000e+00, 0.00000000e+00, ...,
2.40449128e-01, -1.08488270e-02, 0.00000000e+00],
[-8.56840371e+00, 0.00000000e+00, 0.00000000e+00, ...,
-2.59241796e-04, -1.12756844e-06, 0.00000000e+00],
...,
[-1.32641613e+01, 0.00000000e+00, 0.00000000e+00, ...,
-5.63659404e+00, 6.50939114e-01, 0.00000000e+00],
[-8.55392716e+00, 0.00000000e+00, 0.00000000e+00, ...,
-2.01206880e-01, 9.61930149e-03, 0.00000000e+00],
[-1.29807876e+01, 0.00000000e+00, 0.00000000e+00, ...,
2.60651472e-04, 4.22693052e-05, 0.00000000e+00]])最後一步是使用訓練過的分類器預測每個影象的標籤。對於這一步驟，對於每個訓練例項（使用向量化程式碼），我們將計算每個類的類概率，並將輸出類標籤分配給具有最高概率的類。
def predict_all(X, all_theta):
rows= X.shape[0]
params= X.shape[1]
num_labels= all_theta.shape[0]
# same as before, insert ones to match the shape
X= np.insert(X,0, values=np.ones(rows), axis=1)
# convert to matrices
X= np.matrix(X)
all_theta= np.matrix(all_theta)
# compute the class probability for each class on each training instance
h= sigmoid(X* all_theta.T)
# create array of the index with the maximum probability
h_argmax= np.argmax(h, axis=1)
# because our array was zero-indexed we need to add one for the true label prediction
h_argmax= h_argmax+ 1
return h_argmax現在我們可以使用predict_all函式為每個例項生成類預測，並瞭解分類器的工作情況。
y_pred= predict_all(data['X'], all_theta)
correct= [1 if a== belse 0 for (a, b)in zip(y_pred, data['y'])]
accuracy= (sum(map(int, correct))/ float(len(correct)))
print 'accuracy = {0}%'.format(accuracy* 100)
> accuracy= 97.58%接近98%，相當不錯，邏輯迴歸是一個相對簡單的方法。
本文為編譯文章，作者John Wittenauer，原網址為
http://www.johnwittenauer.net/machine-learning-exercises-in-python-part-4/