机器学习与模式识别:手写数字识别

实验内容

1. 数据集选择
2. Bayes判别分类
3. Fisher 线性判别
4. SVM的线性与非线性分类
5. 不同分类器之间的比较

原始数据集上传至网盘: https://pan.baidu.com/s/1uqmJg7EGxpKR62j-Qbr1ow?pwd=jjrc

实验数据特征提取方法

手写数字样本。每个数字有 $50$ 张图片，选择其中 $40$ 个作为训练集，$10$ 个作为测试集。

首先将含有全部特征信息的手写数字图像从坐标轴中提取出来，将提取出来的书写数字图像进行二值化处理; 将处理后的每个数字图像提取 $5\times 5$ 块模板，每个模块中 1 值像素点与总像素点的比值就是这个模块的特征值。将所有特征值放入 $5\times 5$ 的矩阵。设定阈值 $T = 0.05$，每块内所对应的元素白像素占有率大于 $T$ ，则该块特征取1;否则取0。

选择minst手写数字数据集，因为图片尺寸为 $28\times 28$ 不为5的倍数，使用中心裁剪法将图像裁剪为25x25的大小再进行分块操作与特征提取。

数据处理(特征提取) dataset.m 代码如下:

% ----------
%
% 数据集处理
%
% ----------

function [train_X, train_Y, test_X, test_Y] = load_datasets(train_pc)
    % 参数设置
    T = 0.05;
    kernelSize = 5;
    imgSize = 28;
    sub_counts = floor(imgSize / kernelSize);
    newImgSize = kernelSize * sub_counts;
    train_X = [];
    train_Y = [];
    test_X = [];
    test_Y = [];
    
    for digit = 0:9
        digitFolderPath = fullfile('./mnist', num2str(digit));
        imageFiles = dir(fullfile(digitFolderPath, '*.png'));
        % 读取当前数字的所有图像数据
        images = length(imageFiles);
        for i = 1:images
            imgPath = fullfile(digitFolderPath, imageFiles(i).name);
            % 读取图像
            img = double(imread(imgPath)); 
            % 中心裁剪图像
            croppedImg = centerCropImage(img, newImgSize);
            % 提取图像特征
            features = getFeatures(croppedImg, kernelSize, T);
            %features = img(:);
            if i <= images * train_pc
                % 划分为训练集
                train_X = [train_X; features'];
                train_Y = [train_Y; digit];
            else
                % 划分为测试集
                test_X = [test_X; features'];
                test_Y = [test_Y; digit];
            end
        end
    end
    
    end
    
    
    function croppedImg = centerCropImage(img, newImgSize)
        [rows, cols] = size(img);
        startRow = floor((rows - newImgSize) / 2) + 1;
        startCol = floor((cols - newImgSize) / 2) + 1;
        croppedImg = img(startRow:startRow + newImgSize - 1, startCol:startCol + newImgSize - 1);
    end
    
    function features = getFeatures(img, kernelSize, T)
        [rows, cols] = size(img);
        % 计算分块数目
        numBlocksRow = round(rows / kernelSize);
        numBlocksCol = round(cols / kernelSize);
        % 初始化特征向量
        features = zeros(1, numBlocksRow * numBlocksCol);
    
        blockIndex = 1;
        for i = 1:numBlocksRow
            for j = 1:numBlocksCol
                % 计算分块的起始和结束位置
                startRow = (i - 1) * kernelSize + 1;
                endRow = i * kernelSize;
                startCol = (j - 1) * kernelSize + 1;
                endCol = j * kernelSize;
                % 计算分块的总像素数
                totalPixels = kernelSize * kernelSize;
                % 计算分块内白像素的个数
                whitePixels = sum(sum(img(startRow:endRow, startCol:endCol) >= 250));
                
                % 根据阈值T判断特征取值
                if whitePixels / totalPixels > T
                    features(blockIndex) = 1;
                else
                    features(blockIndex) = 0;
                end
                blockIndex = blockIndex + 1;
            end
        end
        features = features';
    end

朴素Bayes判别分类

理论基础

设 $B_i$ 表示事件: 图片为数字 $i$ ，则由 Bayes 公式有

$$\begin{align*} P(B_i|A) = \frac{P(A|B_i)P(B_i)}{\sum_{j=0}^{9} P(A|B_j)P(B_j)} \end{align*}$$

其中，$P(B_i)$ 在这里是先验概率，在这里等于 0.1。$P(B_i|A)$ 是后验概率，在这里是对于一张手写数字图片(事件 $A$)上的数字是 $d$ ($0\sim 9$ 对应事件 $B_0\sim B_9$)的概率，由于这里是设计基于最小错误率的贝叶斯分类器，故而认为该数字为后验概率最大的数字。

令 $\mathbf{X}$ 表示图片集合，$\mathbf{Y}$ 表示标签集合，则训练数据集可以表示为:

$$\mathbf{T} = \{(\bar{x}_1, y_1), (\bar{x}_2,y_2),\cdots, (\bar{x}_n,y_n)\}$$

其中，$\bar{x}_1,\cdots,\bar{x}_n\in \mathbf{X}$ , $y_1,\cdots,y_n\in \mathbf{Y}$ ，对于任意的 $\bar{x}_i$ 有 $\bar{x}_i = \{x_i^1, x_i^2, \cdots, x_i^m\}$ 意为第 $i$ 张图片的 $m$ 个特征。

对于训练数据集 $P(X,Y)$ 独立同分布，所以有

$$P(X|Y) = \frac{P(X,Y)}{P(Y)}$$

又有先验概率 $P(Y = c_k) = 0.1$ , $k=0,1,\cdots,9$ , 而条件概率

$$P(X=x|Y=c_k) = P(X^1=x^1,X^2=x^2,\cdots,X^m=x^m|Y=c_k)$$

又因为这里数据的条件概率分布是特征条件独立，所以进一步地可以表示为

$$P(X=x|Y=c_k) = \prod_{j=1}^m P(X^j=x^j|Y=c_k)$$

在这里问题里的实际含义是: 对于测试集的任意一张 $28\times 28$ 大小的手写数字图片，最后得到 25 个特征，每一个特征对应于每一个模块的取值。在朴素贝叶斯的假设条件下，这张图片是 1 的概率就是每一个特征都是 1 的特征的概率的累乘。

进一步得到后验概率的计算公式:

$$P(Y=c_k|X=x) = \frac{P(Y=c_k)\prod_j^m P(X^j=x^j|Y=c_k)}{\sum_kP(Y=c_k)\prod_j^m P(X^j=x^j|Y=c_k)}$$

因为是要取最大值，所以可以去掉公分母，得到朴素贝叶斯分类器的判别式

$$y = f(x) = \max_{c_k} P(Y=c_k)\prod_j P(X^j=x^j|Y=c_k)$$

对每一个实验样本，选取前 40 个作为训练集，后 10 个作为测试集。

编写代码 bayesClassifier.m

% ----------
%
% 贝叶斯判别代码
%
% ----------

clc, clear;
% 重新加载数据
load('datasets.mat');

% v = mnistData{8}{1, 2};
% A = reshape(v, 5, 5)';
% disp(A);

% 从每个类别中选择40个样本作为训练集，10个样本作为测试集
train_samples = cell(10, 40);
test_samples = cell(10, 10);

error_count = 0;
for digit = 1:10
    % 从当前类别中随机选择40个样本作为训练集
    all_samples = mnistData{digit};
%     random_indices = randperm(length(all_samples), 40);
%     % disp(random_indices);
%     train_samples{digit} = all_samples(random_indices, :);
    
    train_samples{digit} = all_samples(1:40, :);
    % 剩余的10个样本作为测试集
%     test_indices = setdiff(1:length(all_samples), random_indices);
%     test_samples{digit} = all_samples(test_indices, :);
    test_samples{digit} = all_samples(41:50, :);
    %(test_indices);
end

% 重新计算先验概率和类条件概率
num_classes = 10; % 数字类别数量
num_features = 25; % 特征数量
num_train = 40; % 训练样本的数量
num_test = 10; % 测试样本的数量

true_positives = zeros(1, num_classes); % 正类别被正确分类的样本数量
false_positives = zeros(1, num_classes); % 负类别被错误分类成正类别的样本数量
false_negatives = zeros(1, num_classes); % 正类别被错误分类成负类别的样本数量

for x = 1:10
    for y = 1:10
        % 获取第i个测试样本的特征向量
        test_sample = test_samples{x}{y, 2};
        
        prior_prob = zeros(1, num_classes); % 先验概率
        class_cond_prob = zeros(num_features, num_classes); % 类条件概率
        pij = []; % i类的样本第j个特征为1的概率
        for i = 1:num_classes
            % 计算先验概率
            prior_prob(i) = 0.1;
            for j = 1:num_features % 每个数字图片提取出来的特征数
                sum = 0;
                for k = 1:num_train % 每个类别下训练样本的个数
                    i_feature = train_samples{i}{k, 2}; % 获取第k个训练样本的特征向量
                    sum = sum + i_feature(j);
                end
                % disp(sum);
                pij(i,j) = (sum + 1) / (num_train + 2); % 计算概率估计值即Pj(ωi)，注意拉普拉斯平滑处理
            end
        end
        for i = 1:num_classes
            multi = 1;
            for j = 1:num_features % 每个数字图片提取出来的特征数
                if(test_sample(j) == 1)
                    multi = multi * pij(i,j);
                else
                    multi = multi * (1 - pij(i,j));
                end
            end
            class_cond_prob(i) = multi;
        end
        %计算后验概率
        p_class = []; % 后验概率
        sum = 0;
        for i=1:num_classes%数字类别个数
            sum = sum + prior_prob(i) * class_cond_prob(i);
        end
        for i = 1:num_classes % 数字类别个数
            p_class(i) = prior_prob(i) * class_cond_prob(i) / sum;
        end
        [maxval, maxpos] = max(p_class);
        if maxpos == x
            true_positives(x) = true_positives(x) + 1;
        else
            error_count = error_count + 1;
            false_positives(maxpos) = false_positives(maxpos) + 1;
            false_negatives(x) = false_negatives(x) + 1;
        end
    end
end

% 计算准确率（acc）、精确率（precision）、召回率（recall）、F1-score
precision = true_positives ./ (true_positives + false_positives);
recall = true_positives ./ (true_positives + false_negatives);
f1_score = 2 * (precision .* recall) ./ (precision + recall);

% 计算错误率和正确率
error_rate = error_count / 100;
accuracy = 1 - error_rate;

disp(['Accuracy: ', num2str(accuracy)]);
disp(['Precision: ', num2str(mean(precision))]); % Changed here
disp(['Recall: ', num2str(mean(recall))]); % Changed here
disp(['F1 Score: ', num2str(mean(f1_score))]); % Changed here

得到结果:
Accuracy: 0.66
Precision: 0.67991
Recall: 0.66
F1 Score: 0.65476

Fisher判别分类

使用 Fisher 线性判别方法求分类器的步骤:

1. 计算各类的均值向量: $\mu_i = \frac{1}{N_i}\sum_{x\in X_i}x$ ;
2. 计算各类的类内离散矩阵: $S_{wi} = \sum_{x\in X_i}(x-\mu_i)(x-\mu_i)^T$ ;
3. 计算类内总离散矩阵: $S_w = S_{w0}+S_{w1}+\cdots$ ;
4. 计算总离散矩阵的逆矩阵: $S_w^{-1}$ ;
5. 求出向量 $w^* = S_w^{-1}(\mu_1-\mu_0)$ ;
6. 判别函数为: $y=(w^*)^Tx$ ;
7. 求出判别函数的阈值: $w_0 = \frac{(w^*)^T(\mu_0+\mu_1+\cdots)}{2}$ ;
8. 比较 $y$ 值与阈值的大小得出分类。

% ----------
%
% Fisher分类代码
%
% ----------

clc,clear;
load('datasets.mat');
% 初始化

numClasses = 10; % 类别数
numImages = 50; % 每个类别的图像数
numFeatures = 25; % 特征数
trainingSize = 40; % 训练集大小

% 创建训练集和测试集
trainingData = zeros(numClasses * trainingSize, numFeatures);
trainingLabels = zeros(numClasses * trainingSize, 1);
testData = zeros(numClasses * (numImages - trainingSize), numFeatures);
testLabels = zeros(numClasses * (numImages - trainingSize), 1);

for i = 1:numClasses
    for j = 1:numImages
        if j <= trainingSize
            trainingData((i-1)*trainingSize + j, :) = mnistData{i,1}{j,2};
            trainingLabels((i-1)*trainingSize + j) = i;
        else
            testData((i-1)*(numImages - trainingSize) + j - trainingSize, :) = mnistData{i,1}{51-j,2};
            testLabels((i-1)*(numImages - trainingSize) + j - trainingSize) = i;
        end
    end
end

% 使用Fisher线性判别方法进行训练
MdlLinear = fitcdiscr(trainingData, trainingLabels, 'DiscrimType', 'pseudoLinear');
% 对测试集进行预测

predictedLabels = predict(MdlLinear, testData);

% 计算错误率
errorRate = sum(predictedLabels ~= testLabels) / length(testLabels);
fprintf('Error Rate: %.2f%%\n', errorRate * 100);

% 初始化
prior = ones(1, numClasses) / numClasses; % 先验概率

% 使用Fisher线性判别方法进行训练
MdlLinear = fitcdiscr(trainingData, trainingLabels, 'DiscrimType', 'pseudoLinear', 'Prior', prior);

% 对测试集进行预测
predictedLabels = predict(MdlLinear, testData);

% 计算错误率
errorRate = sum(predictedLabels ~= testLabels) / length(testLabels);
fprintf('Error Rate: %.2f%%\n', errorRate * 100);

% 计算混淆矩阵
C = confusionmat(testLabels, predictedLabels);

% 计算准确率（accuracy）
accuracy = sum(diag(C)) / sum(C(:));
fprintf('Accuracy: %.2f%%\n', accuracy * 100);

% 计算精确率（precision）
precision = diag(C) ./ sum(C, 2);
fprintf('Precision: %.2f%%\n', mean(precision) * 100);

% 计算召回率（recall）
recall = diag(C) ./ sum(C, 1)';
fprintf('Recall: %.2f%%\n', mean(recall) * 100);

% 计算F1-score
f1score = 2 * (precision .* recall) ./ (precision + recall);
fprintf('F1-score: %.2f%%\n', mean(f1score) * 100);

多分类支持向量机

二分类支持向量机介绍

对于线性不可分情况引入惩罚因子 $C$ ，于是广义最优分类面问题模型如下:

$$\max_{a} \sum_{j=1}^N a_j - \frac{1}{2}\sum_{i=1}^{N}\sum_{j=1}^{N} y_iy_ja_ia_jK(x_i,x_j) , s.t.\sum_{ i=1}^N y_ia_i=0$$

其中 $0\le a_i\le C$ 。

特征选择

用于训练SVM的特征使用的是图像的完整像素特征，即一张 $28\times 28$ 的图像，它的特征向量的大小为 $1\times 784$ 。将该特征进行标准化处理后即可用于训练SVM。

ECOC编码与多分类SVM

ECOC(Error-Correcting Output Codes)编码是一种纠错输出编码用于将多分类任务高效地转换为多个二分类任务。Mnist数据集有0~9个数字共10分类。对应的ECOC编码如下图:

代码如下:

% ----------
%
% 支持向量机代码
%
% ----------

clc,clear;

[train_X, train_Y, test_X, test_Y] = load_datasets(0.8);

% 核函数选择，可选：'linear','gaussian','rbf','polynomial'
KernelFunction = 'polynomial';
% 惩罚参数C确认
C = 1000;

template = templateSVM(...
    'KernelFunction', KernelFunction, ...
    'PolynomialOrder', 3, ...
    'KernelScale', 'auto', ...
    'BoxConstraint', C, ...
    'Standardize', true);
svm_model = fitcecoc(...
    train_X, ...
    train_Y, ...
    'Learners', template);

% spy(svm_model.BinaryY(1:40:400,:));
% title('ECOC编码');
% yticks(1:10);
% yticklabels(0:9);
% xlabel('分类器数目');

predicted_labels = predict(svm_model, test_X);

% 计算混淆矩阵
C = confusionmat(test_Y, predicted_labels);

% 计算准确率（accuracy）
accuracy = sum(diag(C)) / sum(C(:));
fprintf('Accuracy: %.2f%%\n', accuracy * 100);

% 计算精确率（precision）
precision = diag(C) ./ sum(C, 2);
fprintf('Precision: %.2f%%\n', mean(precision) * 100);

% 计算召回率（recall）
recall = diag(C) ./ sum(C, 1)';
fprintf('Recall: %.2f%%\n', mean(recall) * 100);

% 计算F1-score
f1score = 2 * (precision .* recall) ./ (precision + recall);
fprintf('F1-score: %.2f%%\n', mean(f1score) * 100);

下表为使用全部特征进行训练、测试得到的结果。

核函数	惩罚参数C	准确率(%)	精确率(%)	召回率(%)	F1-score(%)
线性	1	86.00	86.00	87.83	85.77
线性	100	82.00	82.00	84.35	81.44
线性	1000	82.00	82.00	84.35	81.44
高斯	1	71.00	71.00	83.30	73.64
高斯	100	72.00	72.00	82.89	73.80
高斯	1000	72.00	72.00	82.89	73.80
三次多项式	1	87.00	87.00	87.18	86.58
三次多项式	100	87.00	87.00	87.18	86.58
三次多项式	1000	87.00	87.00	87.18	86.58

由上表得，三次多项式作为核函数效果最佳，且乘法参数C取值对评估结果没有影响。但如果选取原始特征提取方法(图像被分为5x5个块，一共提取了25个特征)，准确率将有所下降，仅能达到 60% 左右。

各分类器间的比较

无论是贝叶斯判别还是Fisher分类，两者改为多分类方法比较容易。SVM是一个性能很好的二分类算法，然而在进行多分类任务时需要多个SVM才能进行，这导致SVM在多分类任务中的准确率下降。本次实验如果全部使用提取特征后的数据来训练模型，Fisher判别表现最好，准确率在74%；其次是贝叶斯判别，准确率在66%. SVM分类效果最差，准确率为61%.

总的来说，分类算法的选择更多取决于数据集。如果数据集规模较大，且基本线性可分，使用贝叶斯或是Fisher判别效率更高，反之应使用SVM处理更加复杂的非线性分类任务。