手写数字识别(机器学习)

一：一对多分类方法(one-vs-all)

这里先上的代码，想看原理可以到代码下面。

在数据集中，y的取值为1~10，y=10表示当前数字为0

首先读取数据,并对数据进行切分。

import pandas as pd
import numpy as np
import matplotlib

matplotlib.use('tKAgg')
import matplotlib.pyplot as plt

file_path = "D:\\JD\\Documents\\大学等等等\\自学部分\\机器学习自学画图\\手写数字识别\\ex3data1.xlsx"
data = pd.read_excel(file_path)
row_X = data.iloc[:, :-1]
row_y = data.iloc[:, -1]
print(row_X.shape, row_y.shape)

(4999, 400) (4999,)

我们来看一下是否显示图像：

def plot_an_image(X):
    pick_one = np.random.randint(len(X))  # 使用len(X)确保随机选择的索引在数据范围内
    image = X.iloc[pick_one, :]
    plt.figure(figsize=(4, 4))  # 设置图像大小
    plt.imshow(image.values.reshape(20, 20).T, cmap='gray_r')  # image.values用于提取数据
    plt.xticks([])
    plt.yticks([])
    plt.show()  # 显示图像


def plot_100_image(X):
    # 随机选择100个样本
    pick_100 = np.random.choice(len(X), 100, replace=False)

    # 创建一个8x8的图形和10x10的子图网格
    fig, axes = plt.subplots(nrows=10, ncols=10, figsize=(8, 8))

    # 确保axes是一个2D数组（10x10）
    axes = axes.flatten()

    for i, idx in enumerate(pick_100):
        image = X.iloc[idx, :]
        # 将图像数据显示到相应的子图中
        axes[i].imshow(image.values.reshape(20, 20).T, cmap='gray_r')
        axes[i].axis('off')  # 关闭坐标轴

    plt.tight_layout()  # 自动调整子图参数
    plt.show()

调用显示100个图片的函数：

plot_100_image(row_X)

构建模型：

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


def costFunction(theta, X, y, lamda):  # 损失函数，要求它最小
    A = sigmoid(X @ theta)  # 预测值
    first = y * np.log(A)  # 损失函数部分
    second = (1 - y) * np.log(1 - A)  # 损失函数部分
    reg = theta[1:] @ theta[1:] * lamda / (2 * len(X))  # 正则项
    return -np.sum(first + second) / len(X) + reg


def gradient_reg(theta, X, y, lamda):
    '''
    :param theta: 要训练的参数
    :param X:    数据集
    :param y:     结果集
    :param lamda:   正则化参数
    :return:
    '''
    reg = theta[1:] * lamda / len(X)
    reg = np.insert(reg, 0, values=0, axis=0)  # 正则项梯度(导数)
    first = (X.T @ (sigmoid(X @ theta) - y)) / len(X)  # 损失函数梯度(导数)
    return reg + first


X = np.insert(row_X, 0, values=1, axis=1)
y = row_y[:]


def one_vs_all(X, y, lamda, K):
    '''
    :param X: 数据集
    :param y:  真实值
    :param lamda: 正则化项参数
    :param K:    标签个数
    :return:
    '''
    n = X.shape[1]  # 维度
    theta_all = np.zeros((K, n))
    for i in range(1, K + 1):
        theta_i = np.zeros(n,)
        res = minimize(costFunction,
                       x0=theta_i,
                       args=(X, y == i, lamda),
                       method='TNC',
                       jac=gradient_reg)
        theta_all[i - 1, :] = res.x
    return theta_all


lamda =1
K=10
theta_final = one_vs_all(X,y,lamda,K)
print((theta_final))


def predict(X,theta_final):
    h = sigmoid(X@theta_final.T)
    h_argmax = np.argmax(h,axis=1)
    return h_argmax+1
y_pred = predict(X,theta_final)
acc = np.mean(y_pred==y)
print(acc)

[[-2.37986344e+00 0.00000000e+00 0.00000000e+00 ... 1.30436171e-03
-7.36236836e-10 0.00000000e+00]
[-3.18544317e+00 0.00000000e+00 0.00000000e+00 ... 4.45792974e-03
-5.08241654e-04 0.00000000e+00]
[-4.79780226e+00 0.00000000e+00 0.00000000e+00 ... -2.86754653e-05
-2.47292717e-07 0.00000000e+00]
...
[-7.98708564e+00 0.00000000e+00 0.00000000e+00 ... -8.95567756e-05
7.21960016e-06 0.00000000e+00]
[-4.57327215e+00 0.00000000e+00 0.00000000e+00 ... -1.33515788e-03
9.98385662e-05 0.00000000e+00]
[-5.40538138e+00 0.00000000e+00 0.00000000e+00 ... -1.16772407e-04
7.89384343e-06 0.00000000e+00]]
0.9445889177835567

可以看到准确率为0.9445889177835567

这里大家可能不知道minimize函数。这里给大家讲一下：

minimize 函数是 SciPy 库中的一个优化函数，主要用于求解最优化问题。它能够找到给定目标函数的最小值，并返回对应的变量值。理解 minimize 函数的工作原理对于优化模型和算法至关重要。下面，我会详细讲解 minimize 函数的各个参数及其使用方法。

`minimize` 函数的参数详解

fun:
- 这是一个必须的参数，指定了需要最小化的目标函数。该函数必须接受参数并返回一个标量值（即目标函数的值）。
- 在代码中，fun 是 costFunction。它计算了逻辑回归的代价（损失）函数，计算了预测值与真实值之间的误差以及正则化项。
x0:
- 初始猜测值。x0 是一个初始的变量值，用于开始优化过程。优化算法会从这些初始值出发，尝试找到使目标函数最小的最佳值。
- 在代码中，x0 是 theta_i（初始化为全零向量）。它代表了逻辑回归模型的初始参数。
args:
- 这是一个元组，用于将额外的参数传递给目标函数。args 中的参数会被传递给 fun 函数。
- 在你的代码中，args 是 (X, y == i, lamda)，它将数据 X、标签的二值化结果 y == i 和正则化参数 lamda 传递给 costFunction。
method:
- 指定使用的优化算法。SciPy 提供了多种优化算法，例如 'BFGS'、'Nelder-Mead'、'TNC'、'L-BFGS-B' 等等。
- 在代码中，使用的是 'TNC'（Truncated Newton Conjugate Gradient）。这种方法适用于大规模问题，并在处理非线性优化时表现良好。
jac:
- 这是目标函数的梯度（即目标函数对各变量的偏导数）的函数。如果提供，minimize 会利用这些梯度信息来加速优化过程。
- 在代码中，jac 是 gradient_reg，它计算了代价函数的梯度，包括正则化项的梯度。
constraints:
- 用于设置约束条件，通常用于优化问题中需要满足的条件（例如，变量的范围）。这是一个可选参数。
- 代码中没有使用这个参数，因此可以忽略。
bounds:
- 用于设置每个变量的范围（例如，变量的上下界）。这是一个可选参数。
- 代码中没有使用这个参数，因此在你的情况中可以忽略。
tol:
- 优化过程的容差。优化算法在达到设定的精度后停止。容差值越小，算法越精确，但计算可能更耗时。
- 代码中没有设置，使用默认值。
options:
- 这个参数用于设置优化过程中的其他选项，例如最大迭代次数、输出选项等。
- 代码中没有设置这个参数，使用默认值。

`minimize` 函数的工作流程

初始化:
- minimize 函数首先使用 x0 作为初始值来开始优化过程。
计算目标函数:
- 使用 fun 参数提供的函数来计算当前 x 值下的目标函数值。你的目标函数是 costFunction，它计算了逻辑回归的损失。
计算梯度:
- 如果 jac 参数提供了梯度函数，minimize 会计算目标函数的梯度，这有助于加速优化过程。
选择优化算法:
- minimize 根据 method 参数选择的优化算法来进行优化。优化算法会尝试逐步调整 x 的值，以找到使目标函数值最小的点。
迭代优化:
- 优化算法会迭代地更新 x 的值，并计算目标函数和梯度，直到满足停止准则（如最大迭代次数或目标函数变化小于容差值）。
返回结果:
- 优化完成后，minimize 返回一个包含最优解和其他信息的结果对象，包括最优参数 x、目标函数值、退出状态等。

假设你有一个简单的二次函数需要最小化，例如 $f(x)=(x-3)^{2}$

from scipy.optimize import minimize

# 定义目标函数
def objective_function(x):
    return (x - 3)**2

# 初始猜测值
x0 = [0]

# 调用 minimize 函数
result = minimize(objective_function, x0)

# 输出结果
print("最优解:", result.x)
print("目标函数值:", result.fun)

注意minimize要求参数的顺序，x0是objective_function的第一个参数。

假设我们有一个目标函数，目的是找到使得以下函数值最小化的 x：

其中 a 和 b 是我们在计算函数值时需要用到的额外参数。假设 a 和 b 是已知的常数，我们希望优化 x。

1. 定义目标函数

首先，我们定义目标函数 f。这个函数需要额外的参数 a 和 b：

def objective_function(x, a, b):
    return (x - a)**2 + b

2. 传递额外参数

为了优化 x，我们需要将 a 和 b 传递给 objective_function。这里我们使用 scipy.optimize.minimize 函数，其中 args 参数用于传递这些额外的参数。

例如，我们希望将 a 设置为 5，b 设置为 2，初始值 x0 为 0。加上损失函数的梯度可能更快一些。

from scipy.optimize import minimize

# 额外参数
a = 5
b = 2

#定义目标函数导数(关于目标参数的梯度)
# 定义梯度函数
def gradient_function(x):
    return 2 * (x - a)


# 初始猜测值
x0 = 0

# 调用 minimize 函数
result = minimize(objective_function, x0, args=(a, b),method='TNC', jac=gradient_function)

print("最优解:", result.x)
print("最小值:", result.fun)

3. 结果解释

objective_function: 我们定义的目标函数，计算 (x - a)^2 + b。
x0: 优化算法的初始猜测值。
args=(a, b): 传递给目标函数的额外参数 a 和 b。

通过上述代码，minimize 会自动调用 objective_function，并将 x、a 和 b 传递给它。优化过程会在 x 上进行，以最小化目标函数的值。

这里一定要注意：当我们使用 minimize 函数时，args 中的参数顺序需要与目标函数定义中的顺序一致。例如，在上述目标函数中，args 应该按照 (a, b) 的顺序传递，因为目标函数期望 a 和 b 分别作为第二和第三个参数：

一对多原理

可能大家也不知道为什么要用“y == i”表达式。

大家应该理解二分类，如果上面是预测（“是”，“否”），那么X不变，y的值域为（“是”，“否”），只需要写成args=(X,y,lamda)就行了。

那如果有多个呢（假如有3个A，B，C）？我们是不是可以先预测A，讲A=1，BorC=0。这样以此类推。下面给出一些推理过程：

设 $y$ 为标签数组，其中 $y$ 的每个元素表示样本的真实类别

我们要训练的类别是 i。例如，若 i = 1，则我们要训练一个模型来区分类别 1 和其他类别。

布尔数组 $y == i$ 生成一个表示每个样本是否属于类别 i的数组。数学上可以写作：

$\text{label}_{i} = \begin{cases} 1 & \text{if } y_j = i \\ 0 & \text{if } y_j \neq i \end{cases}$

代价函数 $J(\theta)$ 用于计算当前类别 i的分类误差，包括正则化项。可以表示为：

$J(\theta) = -\frac{1}{m} \sum_{j=1}^m \left[ y_j \log(h(\mathbf{x}_j)) + (1 - y_j) \log(1 - h(\mathbf{x}_j)) \right] + \frac{\lambda}{2m} \sum_{k=1}^n \theta_k^2$

其中 $h(\mathbf{x}_j) = \sigma(\mathbf{x}_j^T \theta)$ 是模型的预测概率， $\sigma$ 为 sigmoid 函数。

梯度 $\nabla_\theta J(\theta)$ 的计算公式为：

$\nabla_\theta J(\theta) = \frac{1}{m} \sum_{j=1}^m (h(\mathbf{x}_j) - y_j) \mathbf{x}_j + \frac{\lambda}{m} \theta$

优化算法（例如梯度下降）用于最小化代价函数 $J(\theta)$ ，得到优化后的参数 $\theta$ ：

$\theta_{i} = \text{argmin}_{\theta} \left[ J(\theta) \right]$

对于新的样本，使用优化后的参数进行预测：

$\hat{y} = \text{argmax}_i \left( \text{softmax}(\mathbf{x} \cdot \theta_i^T) \right)$

注意我们得到的 $\theta$ 是一个矩阵：

$\begin{bmatrix} \theta_{1,1} & ... & \theta_{1,n} \\ \vdots &\vdots & \vdots \\ \theta_{10,1} & \cdots&\theta_{10,n} \end{bmatrix}$

一共10行，每i行是对于第i个标签的判定。一行共有n个元素，对应着一条记录的n个特征。

通俗的说，第i行是用来预测这个记录属于i类标签的概率，因为一共有10类标签，所以每一记录都会训练10次，得出每一记录属于每一种标签的概率。(或者说，每一行都会训练一次所有的样本，属于这类标签的是正向，否则负向。最终得到每一类标签的权重参数)

注意：

$X_{m\times n}\cdot \theta_{n \times 10} = M_{m \times 10}$

这个M是结果。每一个记录都会有10列，找到每一列中最大的就是预测值。

def predict(X, theta_final):
    h = sigmoid(X @ theta_final.T)
    h_argmax = np.argmax(h, axis=1)
    return h_argmax + 1


y_pred = predict(X, theta_final)
acc = np.mean(y_pred == y)
print(acc)

np.argmax(h, axis=1) 返回每行中最大值的索引。对于每个样本，h_argmax 是该样本最可能的类别的索引（从 0 开始）

h_argmax 返回的是类别索引（从 0 开始），所以需要加 1 将其转换为从 1 开始的类别标签。

二：Softmax方法

对于这个我上一篇文章有介绍原理

Softmax函数：

$h_\theta(x)_k = \frac{e^{\theta_k^T x}}{\sum_{j=1}^K e^{\theta_j^T x}}$

交叉熵损失度量的是预测的概率分布与实际标签之间的差异。对于样本 i和类别 k，交叉熵损失为：

$-\left[ y_{i,k} \log(h_\theta(x_i)_k) \right]$

当样本 i属于类别 k（即 $y_{i,k} = 1$ ）时，损失函数为：

$-\log(h_\theta(x_i)_k)$

否则（即 $y_{i,k} = 0$ ），损失函数为 0。

将每个样本的损失求和得到总损失:

$\text{Total Loss} = \sum_{i=1}^{m} \sum_{k=1}^{K} \left[ y_{i,k} \log(h_\theta(x_i)_k) \right]$

为了得到每个样本的平均损失，最终的代价函数为：

$J(\Theta) = -\frac{1}{m} \sum_{i=1}^{m} \sum_{k=1}^{K} \left[ y_{i,k} \log(h_\theta(x_i)_k) \right]$

import pandas as pd
import numpy as np
import matplotlib
from scipy.optimize import minimize

matplotlib.use('tkAgg')
import matplotlib.pyplot as plt

file_path = "D:\\JD\\Documents\\大学等等等\\自学部分\\机器学习自学画图\\手写数字识别\\ex3data1.xlsx"
data = pd.read_excel(file_path)
row_X = data.iloc[:, :-1]
row_y = data.iloc[:, -1]
print(row_X.shape, row_y.shape)


def plot_an_image(X):
    pick_one = np.random.randint(len(X))  # 使用len(X)确保随机选择的索引在数据范围内
    image = X.iloc[pick_one, :]
    plt.figure(figsize=(4, 4))  # 设置图像大小
    plt.imshow(image.values.reshape(20, 20).T, cmap='gray_r')  # image.values用于提取数据
    plt.xticks([])
    plt.yticks([])
    plt.show()  # 显示图像


def plot_100_image(X):
    # 随机选择100个样本
    pick_100 = np.random.choice(len(X), 100, replace=False)

    # 创建一个8x8的图形和10x10的子图网格
    fig, axes = plt.subplots(nrows=10, ncols=10, figsize=(8, 8))

    # 确保axes是一个2D数组（10x10）
    axes = axes.flatten()

    for i, idx in enumerate(pick_100):
        image = X.iloc[idx, :]
        # 将图像数据显示到相应的子图中
        axes[i].imshow(image.values.reshape(20, 20).T, cmap='gray_r')
        axes[i].axis('off')  # 关闭坐标轴

    plt.tight_layout()  # 自动调整子图参数
    plt.show()


def softmax(z):
    e_z = np.exp(z - np.max(z, axis=1, keepdims=True))  # 稳定计算softmax
    return e_z / np.sum(e_z, axis=1, keepdims=True)


def costFunction(theta, X, y, lamda, K):
    m = X.shape[0]
    Theta = theta.reshape((K, X.shape[1]))  # 转换theta为K x n的矩阵
    A = softmax(X @ Theta.T)  # 计算预测值
    Y = np.eye(K)[y]  # one-hot 编码真实标签
    reg = lamda * np.sum(Theta[:, 1:] ** 2) / (2 * m)  # 正则项
    return -np.sum(Y * np.log(A)) / m + reg


def gradient_reg(theta, X, y, lamda, K):
    m = X.shape[0]
    Theta = theta.reshape((K, X.shape[1]))
    A = softmax(X @ Theta.T)
    Y = np.eye(K)[y]
    grad = (A - Y).T @ X / m
    reg = lamda * np.concatenate([np.zeros((K, 1)), Theta[:, 1:]], axis=1) / m
    return (grad + reg).ravel()


def one_vs_all(X, y, lamda, K):
    n = X.shape[1]  # 维度
    initial_theta = np.zeros((K * n,))
    res = minimize(costFunction,
                   x0=initial_theta,
                   args=(X, y, lamda, K),
                   method='TNC',
                   jac=gradient_reg)
    return res.x.reshape((K, n))


def predict(X, theta_final):
    h = softmax(X @ theta_final.T)
    return np.argmax(h, axis=1)


X = np.insert(row_X, 0, values=1, axis=1)  # 添加偏置项
y = row_y - 1  # 转换标签从1-10到0-9

lamda = 1
K = 10
theta_final = one_vs_all(X, y, lamda, K)
print(theta_final)

y_pred = predict(X, theta_final)
acc = np.mean(y_pred == y)
print(acc)

(4999, 400) (4999,)
[[ 6.02649610e-01 0.00000000e+00 0.00000000e+00 ... 1.10517407e-03
-3.37000911e-10 0.00000000e+00]
[ 1.60407084e-01 0.00000000e+00 0.00000000e+00 ... 2.69379180e-03
-3.10863107e-04 0.00000000e+00]
[-1.01989081e+00 0.00000000e+00 0.00000000e+00 ... -1.40438599e-04
1.22128124e-07 0.00000000e+00]
...
[-3.11440704e+00 0.00000000e+00 0.00000000e+00 ... -1.46182140e-04
8.23433590e-06 0.00000000e+00]
[-3.81197348e-01 0.00000000e+00 0.00000000e+00 ... -1.62596349e-03
1.21163807e-04 0.00000000e+00]
[-1.26535166e+00 0.00000000e+00 0.00000000e+00 ... -7.32969625e-05
7.50789066e-06 0.00000000e+00]]
0.9623924784956991