使用逻辑回归进行多类分类识别手写数字-优快云博客

本文链接：https://blog.youkuaiyun.com/qq_41476257/article/details/124695197

编程作业 3 - 多类分类

对于此练习，我们将使用逻辑回归来识别手写数字（0到9），我们将扩展我们在练习2中写的逻辑回归的实现，并将其应用于一对一的分类，让我们开始加载数据集，它是在MATLAB的本机格式，所以要加载它在Python，我们需要使用SciPy工具、

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

data = loadmat('ex3data1.mat') # 加载数据
data

{'__header__': b'MATLAB 5.0 MAT-file, Platform: GLNXA64, Created on: Sun Oct 16 13:09:09 2011',
 '__version__': '1.0',
 '__globals__': [],
 '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.]]),
 'y': array([[10],
        [10],
        [10],
        ...,
        [ 9],
        [ 9],
        [ 9]], dtype=uint8)}

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

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

好的，我们已经加载了我们的数据，图像在martirx X中表示为400维向量（其中有5000个）。400维‘特征’是原始20$\times$20图像中每个像素的灰度强度。类标签在向量y中作为表示图像中数字的数字类。

第一个任务是将我们的逻辑回归实现修改为完全向量化（即没有’for’循环）。这是因为向量化代码除了简洁外，还能够利用线性代数优化，并且通常比迭代代码快得多。但是，如果从练习2中看到我们的代价函数已经完全向量化实现了，所以我们可以在这里重复使用相同的实现。

sigmoid 函数

g 代表一个常用的逻辑函数（logistic function）为S形函数（Sigmoid function），公式为：

$g(z)=\frac {1}{1+e^{-z}}$

合起来，我们得到逻辑回归模型的假设函数：

$h_\theta(x)=\frac {1}{1+e^{-\theta^TX}}$

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

代价函数：

$J(\theta)=\frac {1}{m}\sum_{i=1}^m[-y^{(i)}\log(h_\theta(x^{(i)}))-(1-y^{(i)})\log(1-h_\theta(x^{(i)}))]+\frac {\lambda}{2m}\sum_{j=1}^{n}\theta_j^2$
Tip：可以通过np.matrix()函数将一个变量转换为numpy型矩形

def cost(theta, X, y, learningRate):
    # INPUT：参数theta，参数X，标签y，正则化参数 学习率作业中设为1
    # OUTPUT：当前参数值下的交叉损失
    # TODO：根据参数和输入的数据计算交叉损失函数
    
    # STEP1：将theta，X，y，转换为numpy类型的矩阵
    theta = np.matrix(theta)
    X = np.matrix(X)
    y = np.matrix(y)
    
    # STEP2：根据公式计算损失函数（不含正则化）
    cross_cost = np.multiply(-y, np.log(sigmoid(X * theta.T))) - np.multiply((1 - y), np.log(1 - sigmoid(X * theta.T)))
    
    # STEP3：根据公式计算损失函数中的正则化部分
    reg = (learningRate / (2 * len(X))) * np.sum(np.power(theta[1:], 2))
    
    # STEP4：把上两步当中的结果加起来得到整体损失函数
    whole_cost = np.sum(cross_cost) / len(X) + reg # sum 不要忘掉，否则最后准确率差别很大
    
    return whole_cost

如果我们需要使用梯度下降法令这个代价函数最小化，因为我们未对 $\theta_0$ 进行正则化，所以梯度下降算法将分两种情况：

$\theta_0 := \theta_0 - a\frac {1}{m} \sum_{i=1}^{m}[h_\theta(x^{(i)})-y^{(i)}]x_0^{(i)}$
$\theta_j := \theta_j - a\frac {1}{m} \sum_{i=1}^{m}[h_\theta(x^{(i)})-y^{(i)}]x_j^{(i)}+\frac {\lambda}{m}\theta_j$
向量化的梯度函数

def gradient(theta, X, y, learningRate):
    # INPUT：参数theta，参数X，标签y，正则化参数
    # OUTPUT：当前参数值下的梯度
    # TODO：根据参数和输入的数据计算梯度
    
    # STEP1：将theta，X，y，转换为numpy类型的矩阵
    theta = np.matrix(theta)
    X = np.matrix(X)
    y = np.matrix(y)
    
    # STEP2：将theta矩阵拉直（转换为一个向量）
    parameters = int(theta.ravel().shape[1]) # theta数量
    
    # STEP3：计算预测的误差
    error = sigmoid(X * theta.T) - y
    
    # STEP4：根据上面的公式计算梯度
    grad = ((X.T * error) / len(X)).T + ((learningRate / len(X)) * theta)
    
    # STEP5：由j=0时不需要正则化，所以这里重置一下
    grad[0, 0] = np.sum(np.multiply(error, X[:,0])) / len(X)
    
    return np.array(grad).ravel()

现在我们已经定义了代价函数和梯度函数，现在是构建分类器的时候了，对于这个任务。我们有10个可能的类，并且由于逻辑回归只能一次在2个类之间进行分类，我们需要多类分类的策略。在本练习中，我们的任务是实现一对一全分类方法，其中具有k个不同类的标签就有k个分类器，每个分类器在“类别 i”和“不是 i”之间决定。我们将把分类器训练包含在一个函数中，该函数计算10个分类器中的每个分类器的最终权重，并将权值返回为k $\times$ (n + 1) 数组，其中n是参数数量。

from scipy.optimize import minimize

def one_vs_all(X, y, num_labels, learing_rate):
    rows = X.shape[0] # 行数，样本个数
    params = X.shape[1] # 列数，样本的维度
    
    all_theta = np.zeros((num_labels, params + 1)) 
    
    X = np.insert(X, 0, values = np.ones(rows), axis = 1) # 添加一列1
    
    for i in range(1, num_labels + 1):
        theta = np.zeros(params + 1)
        y_i = np.array([1 if label == i else 0 for label in y])
        y_i = np.reshape(y_i, (rows, 1)) # reshape()是数组对象中的方法，用于改变数组的形状
        
        fmin = minimize(fun = cost, x0 = theta, args = (X, y_i, learing_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))

# insert(arr, obj, values, axis):arr原始数组，可一可多，obj插入元素位置，values是插入内容，axis是按行按列插入
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$\times$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([[-2.38144072e+00,  0.00000000e+00,  0.00000000e+00, ...,
         1.30199873e-03, -1.67892252e-10,  0.00000000e+00],
       [-3.19355840e+00,  0.00000000e+00,  0.00000000e+00, ...,
         4.35068309e-03, -4.96682254e-04,  0.00000000e+00],
       [-4.79165620e+00,  0.00000000e+00,  0.00000000e+00, ...,
        -2.87354946e-05, -2.48232928e-07,  0.00000000e+00],
       ...,
       [-8.02400188e+00,  0.00000000e+00,  0.00000000e+00, ...,
        -8.63372814e-05,  7.12725320e-06,  0.00000000e+00],
       [-4.58136316e+00,  0.00000000e+00,  0.00000000e+00, ...,
        -1.35269796e-03,  1.01718435e-04,  0.00000000e+00],
       [-5.30435880e+00,  0.00000000e+00,  0.00000000e+00, ...,
        -1.13434749e-04,  9.23080649e-06,  0.00000000e+00]])

我们现在准备好最后一步，使用训练完毕的分类器预测每个图像的标签。对于这一步，我们将计算每个类的类概率，对于每个训练样本（使用当然的向量化代码），并将输出类标签为具有最高概率的类。

Tip：可以使用np.argmax(),numpy.argmax(array, axis) ⽤于返回⼀个numpy数组中最⼤值的索引值。当⼀组中同时出现⼏个最⼤值时，返回第⼀个最⼤值的索引。

def predict_all(X, all_theta):
    # INPUT：参数值theta，测试数据X
    # OUTPUT：预测值
    # TODO：对测试数据进行预测
    
    # STEP1：获取矩阵的维度信息
    rows = X.shape[0]
    params = X.shape[1]
    num_labels = all_theta.shape[0]
    
    # STEP2：把矩阵X加入一行零元素
    X = np.insert(X, 0, values = np.ones(rows), axis = 1)
    
    # STEP3：把矩阵X和all_theta转换为numpy型矩阵
    X = np.matrix(X)
    all_theta = np.matrix(all_theta)
    
    # STEP4：计算样本属于每一类的概率
    h = sigmoid(X * all_theta.T)
    
    # STEP5：找到每个样本中预测概率最大的值
    h_argmax = np.argmax(h, axis = 1)
    
    # STEP6：因为我们的数组是零索引的，所以我们需要为真正的标签+1
    h_argmax = h_argmax + 1
    
    return h_argmax

现在我们可以使用predict_all函数为每个实例生成类预测，看看我们的分类器是如何工作的。

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