3.正态分布概率模型下的最小错误率贝叶斯决策MATLAB程序代码

一、题目：

【题目】:已知三个类别分别为W1:[0,0]T，[2,1]T,[1,0]T; W2:[-1,1]T，[2,0]T,[-2,-1]T; W3:[0,-2]T，[0,-1]T,[1,-2]T。以协方差矩阵相等和不相等两种情况,基于最小错误率的Bayes分类器进行仿真,画出三类的分界线,并判断样本[-2,2]T属于哪一类?

 

二、Matlab程序代码：

%【题目】:已知三个类别分别为W1:[0,0]T，[2,1]T,[1，0]T; W2:[-1，1]T，[2,0]T,[-2,-1]T; W3:[0,-2]T，[0,-1]T,[1，-2]T。
%以协方差矩阵相等和不相等两种情况,基于最小错误率的Bayes分类器进行仿真,画出三类的分界线,并判断样本[-2,2]T属于哪一类?
 
%基于最小错误率的Bayes的MATLAB程序
%题中所给各类样本数据及观察x
clc
clear
w1=[0 0;2 1;1 0];w2=[-1,1;-2,0;-2,-1];w3=[0,-2;0,-1;1,-2];
x=[-2,2]';  %'表示向量的转置,题目中要求判断的样本x
syms x1 x2% Create symbolic scalar variables x1 and .创建符号标量变量x1和x2,syms: Symbolic Scalar
Xr=[x1,x2];Xc=[x1;x2];%Xr里面是逗号，Xc里面是分号

%算出各类状态的先验概率,
samplenum= length(w1)+length(w2)+length(w3);% Length of largest array dimension最大数组维度的长度length(w1)= length(w2)= length(w3)=3
pw1=length(w1)/samplenum;% pw1 = 0.3333 求先验概率
pw2=length(w2)/samplenum;% pw2 = 0.3333 求先验概率
pw3=length(w3)/samplenum;% pw3 = 0.3333 求先验概率
%算出各类的均值,协方差,及协方差的逆矩阵. mean求平均值，var求方差，std求标准差，cov求协方差，corr与corrcoef求相关系数
u1=mean(w1)';u2=mean(w2)';u3=mean(w3)';%求均值u1=（（0+2+1）/3  （0+1+0）/3）'=(1.0000 0.3333)’
D1=cov(w1);D2=cov(w2);D3=cov(w3);%求协方差
invD1=inv(D1);invD2=inv(D2);invD3=inv(D3);% 求协方差的逆矩阵

%在三类协方差不相等时,分别计算观察x的gi(x) i=1,2,3;决定观察x属于的类别.
g1=-1/2*(x-u1)'*invD1*(x-u1)-1/2*log(det(D1))+log(pw1);%各类的协方差阵不相等时的判别函数式，书中公式（2-91）
g2=-1/2*(x-u2)'*invD2*(x-u2)-1/2*log(det(D2))+log(pw2); %各类的协方差阵不相等时的判别函数式，书中公式（2-91）
g3=-1/2*(x-u3)'*invD3*(x-u3)-1/2*log(det(D3))+log(pw3); %各类的协方差阵不相等时的判别函数式，书中公式（2-91）
g=[g1,g2,g3]; % 定义一个向量g, g = [-64.5228  -12.5228  -20.5228]
w=find(g==max(g)); % 确定向量A中最大值所在位置
% max(g)表示取g的最大值，maxIndex=g==max(g)表示向量g中每个元素与max(g)的值进行比较，如果相等，则向量maxIndex对应位置1；否则，置0。
% 执行以上程序后，向量maxIndex=[0,1,0]
%find(x)表示查找非零元素的线性索引，即非零元素所出的位置，例如X = [1 0 2; 0 1 1; 0 0 4]   k = find(X)   k =（1 5 7 8 9）＇
fprintf('协方差不相等时,x=[-2;2]  belongs to w%d\n',w);%将多个数值和字面文本输出到屏幕。%d整数型，\n换行

%求出三类相邻区域的分界线方程.
W1=-0.5*invD1;W2=-0.5*invD2;W3=-0.5*invD3;%判别函数的Wi,书中公式（2-92）
ww1=invD1*u1;ww2=invD2*u2;ww3=invD3*u3; %判别函数的wi,书中公式（2-93）
w10=-1/2*u1'*invD1*u1-1/2*log(det(D1))+log(pw1); %判别函数的wi0,书中公式（2-94）
w20=-1/2*u2'*invD2*u2-1/2*log(det(D2))+log(pw2); %判别函数的wi0,书中公式（2-94）
w30=-1/2*u3'*invD3*u3-1/2*log(det(D3))+log(pw3); %判别函数的wi0,书中公式（2-94）

%W1与W2的分界线.
rear=w10-w20;%决策方程中的项wi0-wj0，定义rear为wi0-wj0
g12=Xr*round(W1-W2)*Xc+(ww1-ww2)'*Xc+rear; %决策方程，书中公式（2-95），g12 = 4*x1^2 + 22*x1 - 4*x2^2 - 12*x2 + 16。round四舍五入为最近的小数或整数
g12=simplify(g12);% 代数简化g12 =4*x1^2 + 22*x1 - 4*x2^2 - 12*x2 + 16
%Algebraic simplification 代数简化
%syms x a b c
%S = simplify(sin(x)^2 + cos(x)^2)  运行程序后 S=1
%S = simplify(exp(c*log(sqrt(a+b)))) 运行程序后 S =(a + b)^(c/2)

%W1与W3的分界线.
rear=w10-w30; %决策方程中的项wi0-wj0，定义rear为wi0-wj0
g13=Xr*round(W1-W3)*Xc+(ww1-ww3)'*Xc+rear; %决策方程，书中公式（2-95），g13 = 4*x1 + 4*x2 + 4*x1*x2 + x2*(4*x1 - 4*x2) + 4。round四舍五入为最近的小数或整数
g13=simplify(g13);% 代数简化g13 =4*x1 + 4*x2 + 8*x1*x2 - 4*x2^2 + 4

%W2与W3的分界线.
rear=w20-w30; %决策方程中的项wi0-wj0，定义rear为wi0-wj0
g23=Xr*round(W2-W3)*Xc+(ww2-ww3)'*Xc+rear; %决策方程，书中公式（2-95），g23 =16*x2 - 18*x1 + 4*x1*x2 - x1*(4*x1 - 4*x2) - 12 %g23 =16*x2 - 18*x1 + 8*x1*x2 - 4*x1^2 - 12。round四舍五入为最近的小数或整数
g23=simplify(g23);%代数简化g23 =16*x2 - 18*x1 + 8*x1*x2 - 4*x1^2 - 12

%绘制各类数据点及分界线.
figure('name','三类协方差不相等时')%创建图窗窗口，Figure 1:三类协方差不相等时
subplot(2,2,1)%subplot(m,n,p)将当前图窗划分为 m×n 网格，并在 p 指定的位置创建坐标区。subplot(2,2,1)2*2=4个图的第1个图
plot(w1(:,1),w1(:,2),'b+',w2(:,1),w2(:,2),'r*',w3(:,1),w3(:,2),'k.')%w1(:,1)取w1的第一列，'b+'黑色加号，'r*'红色星号，'k.'黑色点
hold on%添加新绘图时保留当前绘图
set(ezplot(g12),'color','red')%设置g12线条的颜色为红色,%ezplot(fun)函数绘图
set(ezplot(g13),'color','blue')%设置g13线条的颜色为蓝色
set(ezplot(g23),'color','magenta')%设置g23线条的颜色为紫色
legend('第一类','第二类','第三类','一二类分界线','一三类分界线','二三类分界线')
title('各类间混合分界线')%第一个小图的标题


subplot(2,2,2)%subplot(m,n,p)将当前图窗划分为 m×n 网格，并在 p 指定的位置创建坐标区。subplot(2,2,2)2*2=4个图的第2个图
plot(w1(:,1),w1(:,2),'b+',w2(:,1),w2(:,2),'r*')%w1(:,1)取w1的第一列，'b+'黑色加号，'r*'红色星号
hold on%添加新绘图时保留当前绘图
set(ezplot(g12),'color','red')%设置g12线条的颜色为红色
legend('第一类','第二类','一二类分界线')%legend(label1,...,labelN)在坐标区上添加图例，'b+'黑色加号为第一类，'r*'红色星号为第二类，

subplot(2,2,3)%subplot(m,n,p)将当前图窗划分为 m×n 网格，并在 p 指定的位置创建坐标区。subplot(2,2,2)2*2=4个图的第3个图
plot(w1(:,1),w1(:,2),'b+',w3(:,1),w3(:,2),'k.')%w1(:,1)取w1的第一列，'b+'黑色加号，'k.'黑色点
hold on%添加新绘图时保留当前绘图
set(ezplot(g13),'color','blue')%设置g13线条的颜色为蓝色
legend('第一类','第三类','一三类分界线')%legend(label1,...,labelN)在坐标区上添加图例，'b+'黑色加号为第一类，'k.'黑色点为第三类

subplot(2,2,4)%subplot(m,n,p)将当前图窗划分为 m×n 网格，并在 p 指定的位置创建坐标区。subplot(2,2,2)2*2=4个图的第4个图
plot(w2(:,1),w2(:,2),'r*',w3(:,1),w3(:,2),'k.')%w1(:,1)取w1的第一列，'r*'红色星号，'k.'黑色点
hold on%添加新绘图时保留当前绘图
%ezplot(g23)%ezplot(fun)函数绘图
set(ezplot(g23),'color','magenta')%设置g23线条的颜色为紫色
legend('第二类','第三类','二三类分界线')%legend(label1,...,labelN)在坐标区上添加图例，'b+'黑色加号为第一类，'k.'黑色点为第三类


%在三类协方差相等时,分别计算观察x的gi(x) i=1,2,3;决定观察x属于的类别.
D=D1+D2+D3;invD=inv(D);%协方差的逆
g1=(x-u1)'*invD*(x-u1);%各类的协方差阵相等时的判别函数，书中公式（2-82），此时决策规则：为了对样本x进行分类，只要计算出x到每类的均值ui的马氏（Mahalanobis）距离平方γ^2，最后把x归于γ^2最小的类别
g2=(x-u2)'*invD*(x-u2);%各类的协方差阵相等时的判别函数，书中公式（2-82），此时决策规则：为了对样本x进行分类，只要计算出x到每类的均值ui的马氏（Mahalanobis）距离平方γ^2，最后把x归于γ^2最小的类别
g3=(x-u3)'*invD*(x-u3);%各类的协方差阵相等时的判别函数，书中公式（2-82），此时决策规则：为了对样本x进行分类，只要计算出x到每类的均值ui的马氏（Mahalanobis）距离平方γ^2，最后把x归于γ^2最小的类别
g=[g1,g2,g3]; % 定义一个向量g, g =[13.4222    3.8222   21.9556]
w=find(g==min(g)); % 确定向量A中最小值所在位置，注意：各类的协方差阵不相等时w=find(g==max(g))找最大
% min(g)表示取g的最大值，minIndex=g==min(g)表示向量g中每个元素与min(g)的值进行比较，如果相等，则向量maxIndex对应位置1；否则，置0。
% 执行以上程序后，向量minIndex=[0,1,0]
%find(x)表示查找非零元素的线性索引，即非零元素所出的位置，例如X = [1 0 2; 0 1 1; 0 0 4]   k = find(X)   k =（1 5 7 8 9）＇
fprintf('协方差相等时,x=[-2;2] belongs to w%d',w);%将多个数值和字面文本输出到屏幕。%d整数型，\n换行

%求出三类相邻区域的分界线方程.
W12=invD*(u1-u2);W13=invD*(u1-u3);W23=invD*(u2-u3);%决策方程w，公式(2-87)
x120=1/2*(u1+u2)-(log(pw1/pw2)*(u1-u2))/((u1-u2)'*invD*(u1-u2));%决策方程x0，公式(2-88)
x130=1/2*(u1+u3)-(log(pw1/pw3)*(u1-u3))/((u1-u3)'*invD*(u1-u3))%决策方程x0，公式(2-88)
x230=1/2*(u2+u3)-(log(pw2/pw3)*(u2-u3))/((u2-u3)'*invD*(u2-u3))%决策方程x0，公式(2-88)

%W1与W2的分界线.
g12=W12'*(Xc-x120);%判别函数的Wi,书中公式（2-86）w^T(x-x0)=0,g12 =2*x1 - (4*x2)/5 + 4/5
g12=simplify(g12);%代数化简,g12 =2*x1 - (4*x2)/5 + 4/5

%W1与W3的分界线.
g13=W13'*(Xc-x130);%判别函数的Wi,书中公式（2-86）w^T(x-x0)=0,g13 =(4*x2)/3 - (4*x1)/15 + 16/15
g13=simplify(g13);%代数化简,g13 =(4*x2)/3 - (4*x1)/15 + 16/15

%W2与W3的分界线.
g23=W23'*(Xc-x230);%判别函数的Wi,书中公式（2-86）书中公式（2-86）w^T(x-x0)=0,g23 =(32*x2)/15 - (34*x1)/15 + 4/15
g23=simplify(g23);%代数化简,g23 =(32*x2)/15 - (34*x1)/15 + 4/15

%绘制各类数据点及分界线.
figure('name','三类协方差相等时')%创建图窗窗口，Figure 2:三类协方差不相等时
subplot(2,2,1)%subplot(m,n,p)将当前图窗划分为 m×n 网格，并在 p 指定的位置创建坐标区。subplot(2,2,1)2*2=4个图的第1个图
plot(w1(:,1),w1(:,2),'b+',w2(:,1),w2(:,2),'r*',w3(:,1),w3(:,2),'k.')%w1(:,1)取w1的第一列，'b+'黑色加号，'r*'红色星号
hold on%添加新绘图时保留当前绘图
set(ezplot(g12),'color','red')%设置g12线条的颜色为红色,%ezplot(fun)函数绘图
set(ezplot(g13),'color','blue')%设置g13线条的颜色为蓝色
set(ezplot(g23),'color','magenta')%设置g23线条的颜色为紫色
legend('第一类','第二类','第三类','一二类分界线','一三类分界线','二三类分界线')
title('各类间混合分界线')%第一个小图的标题

subplot(2,2,2)%subplot(m,n,p)将当前图窗划分为 m×n 网格，并在 p 指定的位置创建坐标区。subplot(2,2,1)2*2=4个图的第2个图
plot(w1(:,1),w1(:,2),'b+',w2(:,1),w2(:,2),'r*')%w1(:,1)取w1的第一列，'b+'黑色加号，'r*'红色星号
hold on%添加新绘图时保留当前绘图
set(ezplot(g12),'color','red')%设置g12线条的颜色为红色
legend('第一类','第二类','一二类分界线')

subplot(2,2,3)%subplot(m,n,p)将当前图窗划分为 m×n 网格，并在 p 指定的位置创建坐标区。subplot(2,2,1)2*2=4个图的第3个图
plot(w1(:,1),w1(:,2),'b+',w3(:,1),w3(:,2),'k.')%w1(:,1)取w1的第一列，'b+'黑色加号，'k.'黑色点
hold on%添加新绘图时保留当前绘图
set(ezplot(g13),'color','blue')%设置g13线条的颜色为蓝色
legend('第一类','第三类','一三类分界线')

subplot(2,2,4)%subplot(m,n,p)将当前图窗划分为 m×n 网格，并在 p 指定的位置创建坐标区。subplot(2,2,1)2*2=4个图的第4个图
plot(w2(:,1),w2(:,2),'r*',w3(:,1),w3(:,2),'k.')%w1(:,1)取w1的第一列，'r*'红色星号，'k.'黑色点
hold on%添加新绘图时保留当前绘图
set(ezplot(g23),'color','magenta')%设置g23线条的颜色为紫色
legend('第二类','第三类','二三类分界线')

三、运行结果：

三类协方差不相等时的分类

三类协方差相等时的分类