模式识别作业6实验报告

最新推荐文章于 2022-11-26 18:05:48 发布

love_lqz

最新推荐文章于 2022-11-26 18:05:48 发布

阅读量1.5k

点赞数 1

分类专栏：模式识别文章标签：模式识别

本文链接：https://blog.csdn.net/weixin_41553159/article/details/106695297

版权

模式识别专栏收录该内容

2 篇文章 1 订阅

订阅专栏

模式识别作业6实验报告

利用贝叶斯分类器实现手写数字识别

基本原理

对每个手写的数字样品，按bxb(如b=5)方式划分，共有b²份。
对每一份内的象素个数进行累加统计，除以每一份内的象素总数，设定阈值T=0.05(可改)，若每一份内的象素占有率大于T则对应的特征值为1，否则为0。
先计算数字 $i$ 的先验概率 $P(\omega_i)=\frac{N_i}{N}$ ，其中 $N_i$ 为数字 $i$ 的样品数， $N$ 为样品总数。
再计算 $P_j(\omega_i)=\frac{\sum_{k=0,X\in\omega_i}^{N_i}x_{kj}+1}{N_i+2}$ ，再计算类条件概率 $P(X|\omega_i)=\prod_{j=0}^{24}P(x_j=\alpha|X\in\omega_i)$ ，其中 $P_j(\omega_i)$ 表示样品 $X$ 属于 $\omega_i$ 类条件下， $X$ 的第 $j$ 个分量为1的概率估计值， $\alpha=0\,or\,1$ .
利用贝叶斯公式求后验概率 $P(\omega_i|X)=\frac{P(\omega_i)P(X|\omega_i)}{\sum_{k=0}^9P(\omega_k)P(X|\omega_k)}，i=0,1,\cdots,9.$
后验概率的最大值的类别0-9就是手写数字的所属类别。

实验数据采集

从网上下载，来源：
https://blog.csdn.net/qq_25005311/article/details/97910815
有4000张手写数字图片，大小为28×28，数字0-9分别有400张。部分图片如下：

算法流程

算法的思维导图为：
在这里插入图片描述

实验结果分析

取不同的样本个数，设定不同的划分数b和阈值T，得到数字识别正确率的百分比如下两个表格：

表格1

样本总数\参数	b=7,T=0.06	b=7,T=1/16	b=7,T=2/16	b=7,T=3/16	b=7,T=4/16	b=7,T=5/16
100	86%	87%	86%	82%	85%	83%
500	72.4%	73.4%	77.6%	76.4%	74.8%	73.8%
1000	70%	71.7%	74%	75.1%	74.2	75.2%
2000	69.1%	71.65%	73.25%	74.1%	72.25%	72.6%
3000	69.53%	71.37%	73.33%	74.33%	72.97%	73.77%
4000	70.03%	71.92%	73.72%	74.42%	72.72%	73.12%

表格2

样本总数\参数	b=7,T=6/16	b=7,T=7/16	b=7,T=8/16	b=14,T=1/4	b=14,T=2/4	b=28,T=0
100	85%	84%	81%	95%	92%	97%
500	76.6%	77%	70%	85.2%	83.6%	85.6%
1000	76.4	75.8%	72.5%	83.8%	83.1%	84.7%
2000	73.75%	73.6%	69.95%	81.35%	80.65%	82.45%
3000	74.33%	73.7%	69.83%	81.97%	81.03%	82.73%
4000	73.52%	72.9%	69.67%	81.7%	80.65%	82.27%

上表可以看到，当b=28，T=0，即对图片不做分块处理，阈值为0时识别率最高，基本在82%以上。下图为样本总数等于100时的数字识别矩阵，可以看到，错误识别三个数字：

识别错误的图片分别为：
在这里插入图片描述
可以看到，即便肉眼很容易看出的数字，算法还是识别错了，说明Bayes分类器识别正确率不高。下图为样本总数等于200时的数字识别矩阵：

总结：Bayes分类器的优点和缺点为：

优点：贝叶斯算法高效，易于实现。
缺点：分类性能不一定很高，对于4000个样本精度才82%左右。改进之处：提高准确率。

程序代码

clear;clc;
tune_b = 28; % 需要调的参数1
tune_T = 0;  % 需要调的参数2
A = [];N1=9; % 每个数字样本个数为N1+1，需要调的参数3
N = 10*(N1+1); % 样本总数为N
for i=0:9
    B = [];
    for j=0:N1
        x = im2double(imread(strcat('D:\kp_matlab\hw6-number_recog\01\images4000\',...
            num2str(i),'_', num2str(j),'.bmp')));
        b = tune_b; % b*b方式划分，影响正确率
        a = size(x,1)/b; % 每份的行数
        a1 = a*ones(1,b);
        x = mat2cell(x,a1,a1); % 将28*28的数值矩阵分块，分成b*b
        x = reshape(x,b^2,1); % 将矩阵转化为列向量
        C = [];
        for k = 1:b^2
            c = (cell2mat(x(k))~=0);
            count = sum(c(:));
            t = count/a^2;
            T = tune_T; % 修改阈值提高正确率
            if t>T
                t=1;
            else t=0;
            end
            C = [C;t]; % 构造特征指标列向量
        end
        B = [B,C]; % 拼接每个样本的指标列向量
    end
    A = cat(3,A,B); % 每个数字样本个数都相同，构造数字0-9的三维特征矩阵
end

P1 = (N1+1)/N; % 先验概率
for i = 1:10
    for j=1:b^2
        s = sum(A(j,:,i))+1;
        P(i,j) = s/((N1+1)+2); % 计算样品X属于wi类条件下，X的第j个分量为1的概率估计值
    end
end

p = ones(10,1);
for h=1:10
    for k=1:(N1+1)
        for i=1:10
            for j=1:b^2
                if A(j,k,h)==1
                    p(i)=p(i)*P(i,j);
                else p(i)=p(i)*(1-P(i,j));
                end
            end
        end
        p = P1*p; % 利用贝叶斯公式求后验概率，分母都相同故忽略不计
        [maxval,index] = max(p); % 后验概率的最大值的类别0-9就是手写数字的所属类别
        I(h,k) = index; % 行h表示数字0-9，列k表示样本
        p = ones(10,1);
    end
end

I = I-1;
disp(I); % 显示数字识别矩阵
for i=1:10
    z(i) = sum((I(i,:)~=i-1));
end
disp(sum(z(:))); % 计算识别错误的总个数
disp((N-sum(z(:)))/N); % 计算正确率

利用逻辑回归实现手写数字识别

基本原理

数据：从网上下载了一个mnistData.mat文件，将其分为训练集Train_X，大小为 $50000\times784$ ，标签为Train_label，大小为 $50000\times1$ ，测试集为 $10000\times784$ ，标签为 $10000\times1$ ，其中50000为训练集样本个数， $784=28\times28$ 为手写数字图片的大小。标签值为1-10，其中标签10代表数字0，即数字0被标记为10. 将Train_X表示为：
$X=\left[\begin{array}{c} -\left(x^{(1)}\right)^{T}- \\ -\left(x^{(2)}\right)^{T}- \\ \vdots \\ -\left(x^{(m)}\right)^{T}- \end{array}\right]$

因为有10个类，需要训练10个独立的logistic回归分类器。为了提高效率，必须确保代码的矢量化良好。logistic回归的代价函数为
$J(\theta)=\frac{1}{m} \sum_{i=1}^{m}\left[-y^{(i)} \log \left(h_{\theta}\left(x^{(i)}\right)\right)-\left(1-y^{(i)}\right) \log \left(1-h_{\theta}\left(x^{(i)}\right)\right)\right]$

其中 $h_{\theta}\left(x^{(i)}\right)=g(\theta^Tx^{(i)})$ ，g为sigmoid函数。易知
$\theta=\left[\begin{array}{c} -\left(x^{(1)}\right)^{T} \theta- \\ -\left(x^{(2)}\right)^{T} \theta- \\ \vdots \\ -\left(x^{(m)}\right)^{T} \theta- \end{array}\right]=\left[\begin{array}{c} -\theta^{T}\left(x^{(1)}\right)- \\ -\theta^{T}\left(x^{(2)}\right)- \\ \vdots \\ -\theta^{T}\left(x^{(m)}\right)- \end{array}\right]$

其中 $\theta=\left[ \begin{matrix} \theta_0& \theta_1& \cdots& \theta_n\\ \end{matrix} \right]^T$ 为需要求的参数， $\theta_0$ 为偏置。

利用梯度下降法求得使代价函数获得最小值所对应的参数 $\theta$ ，其中梯度为
$\frac{\partial J}{\partial \theta_{j}}=\frac{1}{m} \sum_{i=1}^{m}\left(\left(h_{\theta}\left(x^{(i)}\right)-y^{(i)}\right) x_{j}^{(i)}\right)$

向量化后得到：
$\begin{aligned}\left[\begin{array}{c} \frac{\partial J}{\partial \theta_{0}} \\ \frac{\partial J}{\partial \theta_{1}} \\ \frac{\partial J}{\partial \theta_{2}} \\ \vdots \\ \frac{\partial J}{\partial \theta_{n}} \end{array}\right] &=\frac{1}{m}\left[\begin{array}{c} \sum_{i=1}^{m}\left(\left(h_{\theta}\left(x^{(i)}\right)-y^{(i)}\right) x_{0}^{(i)}\right) \\ \sum_{i=1}^{m}\left(\left(h_{\theta}\left(x^{(i)}\right)-y^{(i)}\right) x_{1}^{(i)}\right) \\ \sum_{i=1}^{m}\left(\left(h_{\theta}\left(x^{(i)}\right)-y^{(i)}\right) x_{2}^{(i)}\right) \\ \vdots \\ \sum_{i=1}^{m}\left(\left(h_{\theta}\left(x^{(i)}\right)-y^{(i)}\right) x_{n}^{(i)}\right) \end{array}\right]\\ &=\frac{1}{m} \sum_{i=1}^{m}\left(\left(h_{\theta}\left(x^{(i)}\right)-y^{(i)}\right) x^{(i)}\right) \\ &=\frac{1}{m} X^{T}\left(h_{\theta}(x)-y\right) \end{aligned}$

其中 $h_{\theta}(x)-y=\left[ \begin{array}{c} h_{\theta}(x^{(1)})-y^{(1)}\\ h_{\theta}(x^{(2)})-y^{(2)}\\ \vdots\\ h_{\theta}(x^{(m)})-y^{(m)}\\ \end{array} \right]$

正则化后的代价函数为：
$J(\theta)=\frac{1}{m} \sum_{i=1}^{m}\left[-y^{(i)} \log \left(h_{\theta}\left(x^{(i)}\right)\right)-\left(1-y^{(i)}\right) \log \left(1-h_{\theta}\left(x^{(i)}\right)\right)\right]+\frac{\lambda}{2 m} \sum_{j=1}^{n} \theta_{j}^{2}$

正则化后的梯度为：
$\begin{array}{ll} \frac{\partial J(\theta)}{\partial \theta_{0}}=\frac{1}{m} \sum_{i=1}^{m}\left(h_{\theta}\left(x^{(i)}\right)-y^{(i)}\right) x_{j}^{(i)} & \text { for } j=0 \\ \frac{\partial J(\theta)}{\partial \theta_{j}}=\left(\frac{1}{m} \sum_{i=1}^{m}\left(h_{\theta}\left(x^{(i)}\right)-y^{(i)}\right) x_{j}^{(i)}\right)+\frac{\lambda}{m} \theta_{j} & \text { for } j \geq 1 \end{array}$

于是梯度下降算法为：
Repeat until convergence{
$\theta_{0}:=\theta_{0}-a \frac{1}{m} \sum_{i=1}^{m}\left(\left(h_{\theta}\left(x^{(i)}\right)-y^{(i)}\right) x_{0}^{(i)}\right)$
$\theta_{j}:=\theta_{j}-a\left[\frac{1}{m} \sum_{i=1}^{m}\left(h_{\theta}\left(x^{(i)}\right)-y^{(i)}\right) x_{j}^{(i)}+\frac{\lambda}{m} \theta_{j}\right]$
for $\ldots n$
}

通过梯度下降迭代训练多个正则化logistic回归分类器来实现one vs all分类，每个分类器对应于数据集中的10个类，代码应该返回矩阵中的所有分类器参数 $\theta$ ，这里使用fmincg函数来迭代计算使得Cost函数取得最小值所对应的参数 $\theta$ .

训练完one-vs-all分类器之后，现在可以使用它来预测给定（测试集）图像中包含的数字。对于每个输入，使用经过训练的logistic回归分类器计算它属于每个类的“概率”。one-vs-all预测函数将选择对应的logistic回归分类器输出最大概率的类，并返回类标签（1、2、…、10）作为输入示例的预测。

实验结果分析

数字的展示：

准确率为：90.92%
在这里插入图片描述
可以发现准确率比Bayes方法提升大约8.65%，但运行时间比Bayes慢约20s. 改进之处：继续提高准确率，缩短运行时间。

Logistic	Bayes
92.10s	72.20s

程序代码

主文件

%% 初始化
clear ; close all; clc
tic;
%% 设置参数，划分训练集和测试集
num_labels = 10;          % 10个labels, 从1到10，label 10代表数字0  
load('mnistData.mat')
X = Data(:,1:784);
y = Data(:,785);
X = double(X);
y = double(y);
y(y==0)=10;
% 分为50000组训练集和10000组测试集
for i = 1:10
   Train_X((i-1)*5000+(1:5000),:) = X((i-1)*6000+(1:5000),:); 
   Train_label((i-1)*5000+(1:5000),:) = y((i-1)*6000+(1:5000),:); 
   Test_X((i-1)*1000+(1:1000),:) = X((i-1)*6000+(5001:6000),:);
   Test_label((i-1)*1000+(1:1000),:) = y((i-1)*6000+(5001:6000),:); 
end

%% =========== Part 1: 可视化部分手写数字 =============
X = Train_X;
y = Train_label;
m = size(X, 1); % 总样本数

% 展示100个图片
rand_indices = randperm(m);
sel = X(rand_indices(1:100), :);

displayData(sel); % 调用函数displayData

%% ============ Part 2: 向量化逻辑回归并迭代求参 ============
lambda = 0.1;
[all_theta] = oneVsAll(X, y, num_labels, lambda); % 调用oneVsAll函数

%% ================ Part 3: 对测试集进行测试 ================
pred = predictOneVsAll(all_theta, Test_X);
fprintf('\nTest Set Accuracy: %f\n', mean(double(pred == Test_label)) * 100); % 计算准确率
toc;

展示图片

function [h, display_array] = displayData(X, example_width)
% Set example_width automatically if not passed in
if ~exist('example_width', 'var') || isempty(example_width) 
	example_width = round(sqrt(size(X, 2))); % 对浮点数进行四舍五入
end

% 灰度图像
colormap(gray); % 输出一个灰色系的曲面图，用gray矩阵映射当前图形的色图

% 计算行和列
[m,n] = size(X);
example_height = (n / example_width);

% Compute number of items to display
display_rows = floor(sqrt(m)); % 向下舍入为最接近的整数
display_cols = ceil(m / display_rows); % 向上舍入为最接近的整数

% Between images padding
pad = 1;

% Setup blank display
display_array = - ones(pad + display_rows * (example_height + pad), ...
                       pad + display_cols * (example_width + pad)); % 1代表边缘或边界线

% Copy each example into a patch on the display array
curr_ex = 1;
for j = 1:display_rows
	for i = 1:display_cols
        if curr_ex > m
            break;
        end
		
		% Get the max value of the patch
		max_val = max(abs(X(curr_ex, :))); % 归一化
        X1 = reshape(X(curr_ex, :), example_height, example_width);
		display_array(pad + (j - 1) * (example_height + pad) + (1:example_height), ...
		              pad + (i - 1) * (example_width + pad) + (1:example_width)) = ...
						X1' / max_val;
		curr_ex = curr_ex + 1;
	end
	if curr_ex > m 
		break; 
	end
end

% 展示图片
h = imagesc(display_array, [-1 1]); % 颜色范围[-1,1]

% 不显示坐标轴
axis image off

drawnow;

end

逻辑回归代价函数

function [J, grad] = lrCostFunction(theta, X, y, lambda)
% 初始化一些有用的参数
m = length(y); % 训练样本个数
J = 0;
grad = zeros(size(theta));

theta1 = [0;theta(2:end)]; % 先把theta(1)拿掉，不参与正则化
J1 = -y.*log(sigmoid(X * theta))-(1-y).*(log(1-sigmoid(X * theta))); % logistic回归的Cost函数
J = 1/m * sum(J1(:)) + lambda/(2*m) * theta1' * theta; % 正则化的logistic回归Cost函数
grad = (X'*(sigmoid(X * theta)-y))/m + lambda/m * theta1; % 正则化梯度

grad = grad(:);
end

调参

function [all_theta] = oneVsAll(X, y, num_labels, lambda)
m = size(X, 1);
n = size(X, 2);

all_theta = zeros(num_labels, n + 1);

% 增加一列全是1
X = [ones(m, 1) X]; % 使得维度一致，避免乘积出错

initial_theta = zeros(n + 1, 1); % 初始化theta
options = optimset('GradObj', 'on', 'MaxIter', 50);
% 在fmincg函数中，使用自定义的代价函数，即lrCostFunction，并定义了最大迭代次数为50
% 使用Matlab的fmincg库函数来求解使得代价函数取最小值的模型参数theta
for c = 1:num_labels
    all_theta(c,:) = fmincg (@(t)(lrCostFunction(t, X, (y == c), lambda)), ...
        initial_theta, options);
end
end

预测（测试集）

function p = predictOneVsAll(all_theta, X)
m = size(X, 1);
num_labels = size(all_theta, 1);

p = zeros(size(X, 1), 1); % 初始化

X = [ones(m, 1) X]; % 增加一列全是1

H = sigmoid(X * all_theta');
[~,p] = max(H,[],2); % 返回每行最大值的索引位置，也就是预测的数字

end