1、多分类

在本练习中，您将使用逻辑回归和神经网络来识别手写数字（从 0 到 9）。 自动手写数字识别如今已得到广泛应用，从识别邮件信封上的邮政编码（邮政编码）到识别银行支票上的金额。 本练习将向您展示如何将您学到的方法用于此分类任务。
在练习的第一部分中，您将扩展之前的逻辑回归实现并将其应用于一对多分类。

1.1 数据集

import numpy as np  # NumPy 用于数值操作
import pandas as pd  # Pandas 用于数据处理
import matplotlib.pyplot as plt  # Matplotlib 用于绘图
from scipy.io import loadmat  # SciPy 用于加载 MATLAB 文件

ex3data1.mat 中提供了一个数据集，其中包含 5000 个手写数字的训练示例。 .mat 格式意味着数据已以本机 Octave/MATLAB 矩阵格式保存，而不是文本 (ASCII) 格式类似于 csv 文件。 可以使用 load 命令将这些矩阵直接读入您的程序。 加载后，正确尺寸和值的矩阵将出现在程序的内存中。 该矩阵已被命名，因此您无需为其指定名称。

def load_data(path):
    data = loadmat(path)
    X = data['X']
    y = data['y']
    return X,y
    # 这段代码定义了一个名为 load_data 的函数，该函数接受一个参数 path，表示文件路径。
	# 函数的功能是从指定的路径加载数据。首先，它使用 loadmat 函数从指定的路径加载数据文件。
	# 然后，它从加载的数据中提取出两个变量 X 和 y，分别代表输入特征和对应的标签。
	# 最后，函数返回 X 和 y，即输入特征和标签，供后续的数据处理和分析使用。

X, y = load_data('ex3data1.mat')
print(np.unique(y))  # 看下有几类标签
# [ 1  2  3  4  5  6  7  8  9 10]
X.shape, y.shape
# ((5000, 400), (5000, 1))

ex3data1.mat中有5000个训练样本，其中每个训练
示例是数字的 20 像素 x 20 像素灰度图像。 每个像素都由一个浮点数表示，指示该位置的灰度强度。 20 x 20 像素网格被“展开”为 400 维向量。 这些训练示例中的每一个都成为数据矩阵 X 中的一行。这为我们提供了一个 5000 x 400 矩阵 X，其中每一行都是手写数字图像的训练示例。

训练集的第二部分是一个 5000 维向量 y，其中包含训练集的标签。 为了使事情与没有零索引的 Octave/MATLAB 索引更加兼容，我们将数字零映射到值十。 因此，数字“0”被标记为“10”，而数字“1”到“9”则按照其自然顺序被标记为“1”到“9”。
第一个任务是将我们的逻辑回归实现修改为完全向量化（即没有“for”循环）。这是因为向量化代码除了简洁外，还能够利用线性代数优化，并且通常比迭代代码快得多。

1.2 数据可视化

您将首先可视化训练集的子集。 在 ex3.m 的第 1 部分中，代码从 X 中随机选择 100 行，并将这些行传递给 displayData 函数。 此函数将每行映射到 20 像素 x 20 像素的灰度图像，并将图像一起显示。 我们提供了 displayData 函数，我们鼓励您检查代码以了解它是如何工作的。 运行此步骤后，您应该会看到如图 1 所示的图像。

def plot_100_image(X):
    """
    随机打印100个数字
    参数 X: 包含数字图像的数据
    """
    sample_idx = np.random.choice(np.arange(X.shape[0]), 100)
    sample_images = X[sample_idx, :]  # 获取随机选择的数字图像数据
    # 创建一个包含10行10列的图形窗口，每个子图共享 x 轴和 y 轴，图形大小为8x8
    fig, ax_array = plt.subplots(nrows=10, ncols=10, sharey=True, sharex=True, figsize=(8, 8)) 
    # 遍历每个子图
    for row in range(10):
        for column in range(10):
             # 在当前子图中绘制对应索引的数字图像，并将其 reshape 为20x20的矩阵，使用灰度色彩显示
            ax_array[row, column].matshow(sample_images[10 * row + column].reshape((20, 20)),
                                   cmap='gray_r')   
    #具体解释：sample_images[10 * row + column] 表示在第 10 * row + column 行中的样本。然后，使用 .reshape((20, 20)) 将这个一维向量重新塑造为一个 20x20 的二维矩阵，这样可以将其作为图像进行显示。
    plt.xticks([])   # 去除 x 轴的刻度，使图像更美观
    plt.yticks([])   # 去除 y 轴的刻度，使图像更美观
    plt.show() # 显示图形

1.3 矢量化 Logistic 回归

您将使用多个一对多逻辑回归模型来构建多类分类器。 由于有 10 个类，因此您需要训练 10 个单独的逻辑回归分类器。 为了提高训练效率，确保您的代码得到良好的矢量化非常重要。 在本节中，您将实现逻辑回归的矢量化版本，该版本不使用任何 for 循环。 您可以使用上一个练习中的代码作为本练习的起点。

1.3.1 向量化成本函数

我们将首先编写代价函数的矢量化版本。 回想一下，逻辑回归中，代价函数是

为了计算求和中的每个元素，我们必须计算

事实证明，我们可以通过使用矩阵乘法来快速计算所有示例。 让我们将 X 和 θ 定义为

然后，通过计算矩阵乘积 Xθ，我们有

在最后一个等式中，我们用到了一个定理，如果 a 和 b 都是向量，a^Tb=b^Ta，这样我们就可以用一行代码计算出所有的样本。
定义代价函数：

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

def regularized_cost(theta, X, y, l):
    thetaReg = theta[1:]
    first = (-y*np.log(sigmoid(X@theta))) + (y-1)*np.log(1-sigmoid(X@theta))
    reg = (thetaReg@thetaReg)*l / (2*len(X))
    return np.mean(first) + reg

1.3.2 矢量化梯度

回想一下，（非正则化）逻辑回归成本的梯度是一个向量，其中第 j 个元素定义为

为了在数据集上向量化此操作，我们首先明确地写出所有 θ 的所有偏导数

观察到


梯度表示如下：

def regularized_gradient(theta, X, y, l):
    thetaReg = theta[1:]
    first = (1 / len(X)) * X.T @ (sigmoid(X @ theta) - y)
    reg = np.concatenate([np.array([0]),(l/len(X)) * thetaReg])
    return first + reg
    # np.concatenate() 是 NumPy 库中的一个函数，用于沿指定轴连接多个数组或矩阵。
    # 具体来说，它将多个数组按照给定的轴连接起来，生成一个新的数组。

1.4 一对多分类

在练习的这一部分中，您将通过训练多个正则化逻辑回归分类器来实现一对多分类，每个分类器对应数据集中的 K 个类别（图 1）。 在手写数字数据集中，K = 10，但您的代码应该适用于任何 K 值。
您现在应该完成一对多中的代码，为每个类别训练一个分类器。 特别是，您的代码应返回矩阵 θ ∈ R K×(N+1) 中的所有分类器参数，其中 θ 的每一行对应于一个类的学习逻辑回归参数。 您可以使用从 1 到 K 的“for”循环来完成此操作，独立训练每个分类器。

from scipy.optimize import minimize #导入 minimize 函数，它用于最小化（优化）一个目标函数。
def one_vs_all(X, y, l, K):
    '''
    定义了一个名为 one_vs_all 的函数，实现了通用的逻辑回归算法。它有四个参数：
    X：特征矩阵，形状为 (m, n+1)，其中 m 是样本数量，n 是特征数量。这里的特征矩阵包含了截距项，即第一列全为 1。
    y：目标向量，形状为 (m, )，包含了每个样本的分类标签。
    l：正则化参数λ，用于控制正则化的强度。
    K：类别数量，即要进行分类的标签数目。
    '''
    all_theta = np.zeros((K, X.shape[1])) #初始化参数矩阵 all_theta，形状为 (K, n+1)，其中 K 是类别数量，n+1 是特征数量（包括截距项）。初始值为全零。
    '''
    对每个类别循环进行一次逻辑回归训练。
    初始化参数向量 theta，形状与特征矩阵的列数相同。
    将目标向量 y 中的标签转换为二元分类问题。对于当前类别 i，将所有等于 i 的标签设为 1，其余标签设为 0。
    '''
    for i in range(1, K+1):
        theta = np.zeros(X.shape[1])
        y_i = np.array([1 if label == i else 0 for label in y])
        '''
        使用 minimize 函数最小化正则化的成本函数，从而求解当前类别的最优参数 theta。
        fun=regularized_cost 指定成本函数为 regularized_cost，x0=theta 指定初始参数向量为 theta。
        args=(X, y_i, l) 将特征矩阵 X、目标向量 y_i 和正则化参数 l 作为额外参数传递给成本函数。
        method='TNC' 指定优化算法为 TNC 算法（Truncated Newton Conjugate-Gradient）。
        jac=regularized_gradient 指定梯度函数为 regularized_gradient，用于计算成本函数的梯度。
        options={'disp': True} 设置优化选项，其中 'disp': True 表示在优化过程中打印迭代信息。
        将优化得到的最优参数 ret.x 存储在参数矩阵 all_theta 的对应行中。
        '''
        ret = minimize(fun=regularized_cost, x0=theta, args=(X, y_i, l), method='TNC', jac=regularized_gradient, options={'disp': True})
        all_theta[i-1,:] = ret.x
    return all_theta

训练完一对多分类器后，您现在可以使用它来预测给定图像中包含的数字。 对于每个输入，您应该使用经过训练的逻辑回归分类器计算它属于每个类的“概率”。 您的一对多预测函数将选择相应逻辑回归分类器输出最高概率的类，并返回类标签（1、2、…或 K）作为输入示例的预测。

def predict_all(X, all_theta):
    # 计算每个训练实例的每个类别的类别概率
    h = sigmoid(X @ all_theta.T)  # 注意的这里的all_theta需要转置
    # 创建具有最大概率索引的数组
    # 返回沿着指定轴的最大值的索引。
    h_argmax = np.argmax(h, axis=1)
    # 因为我们的数组是从零开始的，所以我们需要为真实标签预测添加1
    h_argmax = h_argmax + 1
    
    return h_argmax

raw_X, raw_y = load_data('E:\python work\machine learning\ML-homework\ex3-neural network\ex3data1.mat') #调用 load_data 函数加载数据集。
#ex3data1.mat 文件中包含了原始特征矩阵 raw_X 和标签向量 raw_y。
'''
在原始特征矩阵 raw_X 中插入一列全为1的偏置项，以便进行逻辑回归。这是为了对应模型中的截距项（intercept term）。
np.insert() 函数用于在数组中插入值，第一个参数是原始特征矩阵 raw_X，第二个参数是要插入的位置，这里是0，表示在第一列插入。
第三个参数是要插入的值，这里是1，表示偏置项。axis=1 表示在列方向插入。
最终得到的特征矩阵 X 的形状为 (5000, 401)，其中5000是样本数量，401是特征数量（包括偏置项）。
'''
X = np.insert(raw_X, 0, 1, axis=1) # (5000, 401)
y = raw_y.flatten()  # 这里消除了一个维度，方便后面的计算 or .reshape(-1) （5000，）
'''
调用 one_vs_all 函数训练一个一对所有（One-vs-All）逻辑回归模型，用于多类别分类任务。
参数 X 是特征矩阵，y 是标签向量，1 是正则化参数 λ，10 是类别数量。
函数返回每个类别的训练参数，存储在参数矩阵 all_theta 中。
'''
all_theta = one_vs_all(X, y, 1, 10)
all_theta  # 每一行是一个分类器的一组参数
'''
使用训练好的参数 all_theta 对训练集特征矩阵 X 进行预测，得到预测的类别标签向量 y_pred。
predict_all 函数根据输入特征矩阵 X 和参数矩阵 all_theta 进行预测，并返回预测的类别标签。
'''
y_pred = predict_all(X, all_theta)
'''
计算预测的准确率。准确率定义为正确预测的样本数量占总样本数量的比例。
使用 NumPy 的 mean() 函数计算预测结果 y_pred 与真实标签 y 相等的比例，
因为当预测值等于真实值时，对应的值为 True，否则为 False。然后求平均值，即准确率。
'''
accuracy = np.mean(y_pred == y)
print ('accuracy = {0}%'.format(accuracy * 100))

结果如下：

2.神经网络

在本练习的前一部分中，您实现了多类逻辑回归来识别手写数字。 然而，逻辑回归无法形成更复杂的假设，因为它只是一个线性分类器。3 在练习的这一部分中，您将使用与之前相同的训练集实现一个神经网络来识别手写数字。 神经网络将能够表示形成非线性假设的复杂模型。 本周，您将使用我们已经训练过的神经网络中的参数。 您的目标是实现前馈传播算法以使用我们的权重进行预测。 在下周的练习中，您将编写用于学习神经网络参数的反向传播算法。提供的脚本 ex3.m 将帮助您逐步完成此练习。

'''
load_weight 是一个函数，用于加载神经网络的权重参数。
loadmat(path) 函数从指定路径加载数据，假设数据是以Matlab格式保存的文件。
加载的数据被赋值给变量 data。
data['Theta1'] 返回加载的数据中的键名为 'Theta1' 的值，即神经网络的第一层到第二层之间的权重参数矩阵。
data['Theta2'] 返回加载的数据中的键名为 'Theta2' 的值，即神经网络的第二层到输出层之间的权重参数矩阵。
函数最终返回两个参数，分别是神经网络的两个权重参数矩阵 Theta1 和 Theta2。
'''
def load_weight(path):
    data = loadmat(path)
    return data['Theta1'], data['Theta2']
theta1, theta2 = load_weight('ex3weights.mat')
theta1.shape, theta2.shape

# 使用 load_data 函数从指定路径加载数据集。假设数据集是以Matlab格式保存的文件。
# 加载的数据集包含特征矩阵 X 和标签向量 y。
X, y = load_data('ex3data1.mat')
# 将标签向量 y 进行扁平化操作，将多维数组降为一维数组。这样做是为了确保标签向量的形状与后续计算的一致性。
y = y.flatten()
# 在特征矩阵 X 中插入一列全为1的偏置项，以便进行神经网络训练。这是为了对应模型中的截距项（intercept term）。
# 使用 NumPy 的 insert() 函数在 X 的第一列（索引为0）插入值为1的偏置项。参数 axis=1 表示在列方向插入。
# 最终得到的特征矩阵 X 的形状为 (m, n+1)，其中 m 是样本数量，n 是特征数量（不包括偏置项）。
X = np.insert(X, 0, values=np.ones(X.shape[0]), axis=1)  # intercept
# 打印输出特征矩阵 X 和标签向量 y 的形状，以确保数据预处理的正确性。
X.shape, y.shape

a1 = X #将输入特征矩阵 X 赋值给变量 a1，表示神经网络的第一层的激活值（也即是输入层的值）。
#计算神经网络的第二层输入值 z2，即第二层（隐藏层）的激活前的值。这一层的输入值是由第一层的激活值 a1 和第一层到第二层之间的权重参数矩阵 theta1 计算得到的。
# @ 表示矩阵乘法，theta1.T 是参数矩阵 theta1 的转置，确保参数矩阵的列数与特征矩阵 a1 的行数相匹配。这样，z2 的每一行对应一个样本，每一列对应一个隐藏单元。
z2 = a1 @ theta1.T 
# 打印输出隐藏层输入值 z2 的形状，以便检查计算结果的正确性。z2 的形状应该是 (m, hidden_units)，其中 m 是样本数量，hidden_units 是隐藏层单元的数量。
z2.shape

# 在隐藏层的输入值矩阵z2中插入一列全为1的偏置项，以便进行神经网络的计算。这是为了对应模型中的截距项。
z2 = np.insert(z2, 0, 1, axis=1) 
# 计算隐藏层的激活值a2，使用 S 型函数（sigmoid 函数）对隐藏层的输入值进行激活。
a2 = sigmoid(z2)
a2.shape
z3 = a2 @ theta2.T
z3.shape
a3 = sigmoid(z3)
a3.shape

'''
算预测类别的向量 y_pred。由于输出层的激活值a3是对每个类别的预测概率，我们选择具有最高预测概率的类别作为最终预测结果。
np.argmax(a3, axis=1) 返回每个样本中预测概率最高的类别的索引，axis=1 表示沿着每行的方向进行操作，即对每个样本找到概率最高的类别。
由于索引是从0开始的，而类别是从1开始的，因此需要加1。
'''
y_pred = np.argmax(a3, axis=1) + 1 
accuracy = np.mean(y_pred == y)
print ('accuracy = {0}%'.format(accuracy * 100))  # accuracy = 97.52%