（1）資料集描述

使用邏輯迴歸來識別手寫數字（0到9）。 將我們之前的邏輯迴歸的實現，擴充套件到多分類的實現。 資料集是MATLAB的本機格式，要載入它到Python，我們需要使用一個SciPy工具。影象在martix X中表示為400維向量（其中有5,000個）, 400維“特徵”是原始20 x 20影象中每個畫素的灰度強度, 類標籤在向量y中作為表示影象中數字的數字類。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.io import loadmat

data = loadmat('ex3data1.mat')
data

data['X'].shape, data['y'].shape

#((5000, 400), (5000, 1))

(2)  邏輯函式

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

（3）損失函式

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

（5）梯度下降

#迴圈函式計算梯度

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 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

#向量化梯度計算

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
    
    # X轉置後 400 * 5000 ，矩陣相加，對應元素相加，輸出[400,]
    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()

（6）計算10個分類模型，分類器訓練

需要計10個可能的類，並且由於邏輯迴歸只能一次在2個類之間進行分類，我們需要多類分類的策略。 其中具有k個不同類的標籤就有k個分類器，每個分類器在“是類別 i”和“不是 i”之間決定，取值為 1  與 0 。 我們將把分類器訓練包含在一個函式中，該函式計算10個分類器中的每個分類器的最終權重，並將權重返回為k ×（n + 1）陣列，其中n是引數數量。

from scipy.optimize import minimize

def all_class_train(X, y, num_labels, learning_rate):
    rows = X.shape[0]
    params = X.shape[1]    
    # k* (n + 1) array for the parameters of each of the k classifiers   10 * 401
    all_theta = np.zeros((num_labels, params + 1))    
    # insert a column of ones at the beginning for the intercept term 在特徵集的第一列插入1
    X = np.insert(X, 0, values=np.ones(rows), axis=1)    
    # labels are 1-indexed instead of 0-indexed  10個標籤，分別計算10個邏輯迴歸模型
    for i in range(1, num_labels + 1):
        #401 個引數
        theta = np.zeros(params + 1)
        #轉換成[1,0,0,10,0,0,0,0 .....]
        y_i = np.array([1 if label == i else 0 for label in 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

(7) 資料預處理

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 label in data['y']])
y_0 = np.reshape(y_0, (rows, 1))

X.shape, y_0.shape, theta.shape, all_theta.shape

#((5000, 401), (5000, 1), (401,), (10, 401))

#theta是一維陣列，因此當它被轉換為計算梯度的程式碼中的矩陣時，它變為（1×401）矩陣。 我們還檢查y中的類標籤，以確保它們看起來像我們想象的一致。

np.unique(data['y'])#看下有幾類標籤

#array([ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10], dtype=uint8)

（8）開始訓練

all_theta = all_class_train(data['X'], data['y'], 10, 1)
all_theta

（9）預測函式

我們現在準備好最後一步 - 使用訓練完畢的分類器預測每個影象的標籤。 對於這一步，我們將計算每個類的類概率，對於每個訓練樣本（使用當然的向量化程式碼），並將輸出類標籤為具有最高概率的類。

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
    # 為了和 data['y']  中的資料比較
    h_argmax = h_argmax + 1
    
    return h_argmax

（10）模型評判

y_pred = predict_all(data['X'], all_theta)
correct = [1 if a == b else 0 for (a, b) in zip(y_pred, data['y'])]
accuracy = (sum(map(int, correct)) / float(len(correct)))
print ('accuracy = {0}%'.format(accuracy * 100))