吳恩達|機器學習作業3.0.邏輯回歸解決多元分類

3.0.邏輯回歸解決多元分類

1）題目：

在本次練習中，你将使用邏輯回歸和神經網絡來識别手寫數字(從0到9)。今天，自動手寫數字識别被廣泛使用，從識别信封上的郵政編碼到識别銀行支票上的金額。這個練習将向您展示如何将您所學的方法用于此分類任務。在第一部分中，将擴充以前的邏輯回歸，并将其應用于one-vs-all分類。

關于資料：本次的資料是以.mat格式儲存的，mat格式是matlab的資料存儲格式，按照矩陣儲存，與numpy資料格式相容，适合于各種數學運算，是以這次主要使用numpy進行運算。ex3data1中有5000個訓練樣例，其中每個訓練樣例是一個20像素×20像素灰階圖像的數字，每個像素由一個浮點數表示，該浮點數表示該位置的灰階強度。每個20×20像素的網格被展開成一個400維的向量。這些每個訓練樣例都變成資料矩陣X中的一行。這就得到了一個5000×400矩陣X，其中每一行都是手寫數字圖像的訓練樣例。訓練集的第二部分是一個包含訓練集标簽的5000維向量y，“0”的數字标記為“10”，而“1”到“9”的數字按自然順序标記為“1”到“9”。

資料集連結: https://pan.baidu.com/s/1cEgQIvehUcLxZ0WVhxcPuQ 提取碼: xejn

2）知識點概括：

• 西瓜書上說多元分類的基本思路為拆解法，即将多分類任務拆為若幹個二分類任務求解。最經典的拆分政策有三種：一對一（OvO）、一對其餘（OvR）和多對多（MvM）。我們這裡用的是One-vs-all，也就是一對其餘。

• 大概意思是把每個分類分别當做1，不是這個分類的當做0，進行上次作業中的邏輯回歸的二分類任務，除了這一點其餘完全和上次作業類似。如果具體算法忘了可以參考

  吳恩達|機器學習作業2.1正則化的Logistic回歸

3）大緻步驟：

首先用scio.loadmat讀取mat資料，分别将訓練樣例和訓練集标簽設定為x和y。再從x中随機選擇100行資料進行可視化。之後再分别計算不帶正則化和帶正則化的代價函數和梯度，這裡可以直接用ex2.1中寫好的函數。然後進行一對多分類，利用for循環對每種數字習得一個帶正則的邏輯回歸分類器，然後将10個分類器的參數組成一個參數矩陣all_theta傳回。最後得到一個5000乘10的預測機率矩陣，找到每一行的機率最大的值，看這個值的位置，得到預測的類别，再和期望值y進行比較，得到精度。

4）關于Python：

• scipy.io庫可以進行mat格式的資料讀寫，讀取mat後，data為dict格式，主要另外添加了一些附加頭資訊。
• np.random.choice隻能從一維數組裡面随機選，可能是因為我自己沒有找到從高維數組裡面選的函數（要是某位大神看到可以留言教一下我~），是以先用np.random.permutation産生一個随機打亂的數組而不重新洗牌x（shuffle函數是重新洗牌），再進行選取。
• plt.figure( ) 相當于建立了一個畫布，括号裡面代表了畫布的編号，使用subplot( nrows, ncols ) 可以在畫布上劃分為nrows * ncols個子圖，imshow可以在每個子圖上畫出數字，注意這裡reshape了每個訓練樣例之後需要進行轉置，圖像才能顯示正常，可能和資料有關。
• 我之前一直以為array和matrix可以通用，因為都可以傳回shape，但是在寫一對多分類函數時發現不行，(401, ) 和 (401,1) 不是一樣的，前面一種情況一般是對一維數組求shape，這時對它求轉置是不起作用的，如果想要轉置的話必須換成矩陣的格式，可以用np.mat轉換。
• np.unravel_index( )的作用是傳回查找值的坐标，np.argmax( ) 的作用是傳回矩陣中的最大值，是以套用起來np.unravel_index(np.argmax(prob[i,:]), prob[i,:].shape)就可以查找第i行的最大值并傳回它所在的位置坐标。

5）代碼與結果：

import numpy as np
import matplotlib.pyplot as plt
import scipy.io as scio #讀取mat檔案
import scipy.optimize as opt

data = scio.loadmat('ex3data1.mat')
x = data['X']
y = data['y']


'''資料可視化'''
s = np.random.permutation(x) #随機重排，但不打亂x中的順序
a = s[:100,:] #選前100行,(100, 400)

#定義可視化函數
def displayData(x):
    plt.figure()
    n = np.round(np.sqrt(x.shape[0])).astype(int)
    #定義10*10的子畫布
    fig, a = plt.subplots(nrows=n, ncols=n, sharex=True, sharey=True, figsize=(6, 6))
    #在每個子畫布中畫出一個數字
    for row in range(n):
        for column in range(n):
            a[row, column].imshow(x[10 * row + column].reshape(20,20).T, cmap='gray')
    plt.xticks([]) #去掉坐标軸
    plt.yticks([])       
    plt.show()

displayData(a)


'''計算代價函數和梯度'''
#sigmoid函數
def sigmoid(x):
    return 1/(1+np.exp(-x))

#原本代價函數
def cost_func(theta, x, y):
    m = y.size
    return -1/m*([email protected].log(sigmoid([email protected]))+(1-y)@np.log(1-sigmoid([email protected])))

#正則化的代價函數，不懲罰第一項theta[0]
def cost_reg(theta, x, y, l=0.1):
    m = y.size
    theta_ = theta[1:] #選取第二項以後的
    return cost_func(theta, x, y) + l/(2*m)*np.sum(theta_*theta_)

#原本的梯度
def gradient_func(theta, x, y):
    m = y.size
    return 1/m*((sigmoid([email protected]))-y).[email protected]

#正則化的的梯度，不懲罰第一項theta[0]
def gradient_reg(theta, x, y, l=0.1):
    theta_ = l/(y.size)*theta
    theta_[0] = 0 #第一項不懲罰設為0
    return gradient_func(theta, x, y) + theta_


'''一對多分類'''
def one_vs_all(x, y, l, K=10):
    all_theta = np.zeros((x.shape[1], K)) #應該是(10, 401)
    for i in range(K):
        iteration_y = np.array([1 if j==i+1 else 0 for j in y]) #第0列到第9列分别對應類别1到10
        p = opt.fmin_ncg(f=cost_reg, fprime=gradient_reg, x0=all_theta[:, i:i+1], args=(x, iteration_y), maxiter=400)
        all_theta[:, i:i+1] = np.mat(p).T
    return all_theta
    
#為x添加了一列常數項 1 ，以計算截距項（常數項）
x = np.column_stack((np.ones(x.shape[0]), x))
lmd = 0.1
all_theta = one_vs_all(x, y, l=lmd)


'''預測與評價'''
#預測的機率矩陣 (5000, 10)
def probability(x, theta):
    return sigmoid([email protected])

#預測的y值
def predict(prob):
    y_predict = np.zeros((prob.shape[0],1))
    for i in range(prob.shape[0]):
        #查找第i行的最大值并傳回它所在的位置,再加1就是對應的類别
        y_predict[i] = np.unravel_index(np.argmax(prob[i,:]), prob[i,:].shape)[0]+1
    return y_predict
        
#精度
def accuracy(y_predict, y=y):
    m = y.size
    count = 0
    for i in range(y.shape[0]):
        if y_predict[i] == y[i]:
            j = 1 
        else:
            j = 0
        count = j+count #計數預測值和期望值相等的項
    return count/m
    
        
prob = probability(x, all_theta)
y_predict = predict(prob)
accuracy(y_predict)
print ('accuracy = {0}%'.format(accuracy(y_predict) * 100))


           

下圖為随機選取100個訓練樣本可視化之後的圖像，和作業裡的那個圖還是有點差别，顔色我調不到那個顔色。gray_r又太淺了。。

這是采用牛頓共轭梯度的結果，不知道為什麼隻有第一次和最後一次是顯示優化成功，其餘均不收斂，但是這種情況下最後的精度還有94.46%，和Andrew Ng給出來的94.9%差不多，也不知道這個結果可不可信。。。