s4_5_q9.m

Duda《模式分类》第二版第1、3、5章部分课后习题和上机题的解答和程序代码

function nearest_classify
% 画出数据点的位置
clear all;close all;clc;
w1=[10,0;0,-10;5,-2];
w2=[5,10;0,5;5,5];
w3=[2,8;-5,2;10,-4];

region=[-12 12 -12 12 100];

subplot(221);
patterns=[w1;w2]';
D=voronoi_regions(patterns, region);
subplot(222);
patterns=[w1;w3]';
D=voronoi_regions(patterns, region);
subplot(223);
patterns=[w2;w3]';
D=voronoi_regions(patterns, region);
subplot(224);
patterns=[w1;w2;w3]';
D=voronoi_regions(patterns, region);


function D = voronoi_regions(patterns, region)
% Make a Voronoi diagram from sample points
% Inputs:
%	patterns	- Input data patterns
%	targets	- Input data targets
%	region	- Decision region vector: [-x x -y y number_of_points]
% 原理：给定已知类别的数据点 Xi （i=1,...,n）的集合 patterns（2*n 阶矩阵）
% 对决策区域均匀分布的N^2个点 Pj（i=1,...,N^2） ，计算出 Pj 与 Xi 的距离
% 根据最近邻规则， Pj 与最靠近的点 Xi 同类
% 用矩阵 D 来标记离 Pj 最近的点 Xi 的编号 i

% 点 Pj 的坐标由向量 x、y 中的元素组合而成
N		= region(5);
x		= linspace (region(1),region(2),N);
y		= linspace (region(3),region(4),N);
[r,c] = size(patterns);

w = c/3;
if w==2
    m1=mean(patterns(:,1:3),2)';
    m2=mean(patterns(:,4:6),2)';
    m=[m1;m2]';
elseif w==3
    m1=mean(patterns(:,1:3),2)';
    m2=mean(patterns(:,4:6),2)';
    m3=mean(patterns(:,7:9),2)';
    m=[m1;m2;m3]';
end

% 用矩阵 D 来标记离 Pj 最近的点 Xi 的编号 i
D = vr(patterns,x,y,N);
% 用矩阵 Dm 来标记离 Pj 最近的类别均值 mi 的编号 i
Dm = vr(m,x,y,N);

% 分别以‘列首尾相接’和‘行首尾相接’的方式，将矩阵 D 转换为向量 d
% 以避免使用 for 循环，减少程序运行时间
i=1:N^2;
% 列首尾相接
% 把 N*N 阶矩阵各点的坐标（i,j）按‘列首尾相接’的方式组成一个 (N^2)*2 坐标阵 [rx,cy]
rx=mod(i,N);rx(find(rx==0))=N;cy=ceil(i/N);
% 将元素为整数坐标值的 [rx,cy] 转换为实际坐标值的 [xn,yn]
xn_c=x(rx)';yn_c=y(cy)';
% 行首尾相接
% 把 N*N 阶矩阵各点的坐标（i,j）按‘行首尾相接’的方式组成一个 (N^2)*2 坐标阵 [cx,ry]
cx=ceil(i/N);ry=mod(i,N);ry(find(ry==0))=N;
xn_r=x(cx)';yn_r=y(ry)';
% 把行列坐标转换为序列编号
ind=sub2ind([N,N],cx,ry);
% 向量 dr 由矩阵 D 按‘行首尾相接’方式转换而成
dr=D(ind);
% 向量 dc 由矩阵 D 按‘列首尾相接’方式转换而成
dc=D(:);

% 画出 Voronoi 网格
vline=0;
if vline
in=1:N^2-1;
lp_c=find(dc(in)~=dc(in+1));
lp_c=lp_c(find(mod(lp_c,N)~=0));
xp_c=xn_c(lp_c);yp_c=yn_c(lp_c);
plot(xp_c,yp_c,'.b');hold on;

lp_r=find(dr(in)~=dr(in+1));
lp_r=lp_r(find(mod(lp_r,N)~=0));
xp_r=xn_r(lp_r);yp_r=yn_r(lp_r);
plot(xp_r,yp_r,'.b');
end

% 画出决策边界
dr_d=ceil(dr/3);
dc_d=ceil(dc/3);
in=1:N^2-1;
lc_d=find(dc_d(in)~=dc_d(in+1));
lc_d=lc_d(find(mod(lc_d,N)~=0));
xc_d=xn_c(lc_d);yc_d=yn_c(lc_d);
plot(xc_d,yc_d,'.k');hold on;

lr_d=find(dr_d(in)~=dr_d(in+1));
lr_d=lr_d(find(mod(lr_d,N)~=0));
xr_d=xn_r(lr_d);yr_d=yn_r(lr_d);
plot(xr_d,yr_d,'.k');

% 画出均值决策边界
drm=Dm(ind);
dcm=Dm(:);
in=1:N^2-1;
lcm=find(dcm(in)~=dcm(in+1));
lcm=lcm(find(mod(lcm,N)~=0));
xcm=xn_c(lcm);ycm=yn_c(lcm);
plot(xcm,ycm,'.g');hold on;

lrm=find(drm(in)~=drm(in+1));
lrm=lrm(find(mod(lrm,N)~=0));
xrm=xn_r(lrm);yrm=yn_r(lrm);
plot(xrm,yrm,'.g');
color=['r';'g';'b';'m';'k';'y'];
% 画出 Voronoi 块
vblock=0;
if vblock
C=color(dc,:);
for i=1:N^2
    plot(xn(i),yn(i),'s','MarkerFaceColor',C(i,:));
end
end

% 画出 patterns 点集和均值点，以及 x、y 轴
for i=1:c
    plot(patterns(1,i),patterns(2,i),'o','MarkerFaceColor',color(ceil(i/3),:));
    j=ceil(i/3);
    plot(m(1,j),m(2,j),'^','MarkerFaceColor',color(ceil(i/3),:));
end

plot(region(1):region(2),zeros(1,region(2)-region(1)+1),'k');
plot(zeros(1,region(4)-region(3)+1),region(3):region(4),'k');
axis(region(1:4));axis square;


function D = vr(patterns,x,y,N)
D		= zeros(N);
[r,c] = size(patterns);
% 首先在 y 轴方向上分别计算出任一 Pj 与 Xi 的坐标差值的平方
% y_dist 是一个 N*c 阶矩阵
y_dist	= (ones(N,1) * patterns(2,:) - y'*ones(1,c)).^2;
% 然后逐列计算该列 Pj 与 Xi 在 x 轴方向上的坐标差值的平方
for i = 1:N
   x_dist = ones(N,1)  * (patterns(1,:)-x(i)).^2; % x_dist 是一个 N*c 阶矩阵
   dist   = abs(x_dist + y_dist);
   [sorted_dist, indices] = min(dist'); % 求出最短距离，indices(k) 即为离 P(i,k) 最近的点 Xp 的编号，即 p = indices(k)
   D(i,:) = indices(1,:);
end