模式识别实验二、贝叶斯分类参数估计matlab代码

u1=[1;2]; %类条件分布概率参数
u2=[5;3];
d1=[1.5 0;0 1];
d2=[3 0;0 3];

n1=100; %训练样本数目
n2=100;
pw1=n1/(n1+n2);
pw2=n2/(n1+n2); %先验概率
x1=mvnrnd(u1,d1,n1); %训练样本
x2=mvnrnd(u2,d2,n2);
xt=[x1;x2];
invsig1=inv(d1);
invsig2=inv(d2);
s1=abs(det(d1));
s2=abs(det(d2));

p1 = mvnpdf(x1, u1', d1);
p2 = mvnpdf(x2, u2', d2);

x1_test=mvnrnd(u1,d1,20); %测试样本
x2_test=mvnrnd(u2,d2,20);
x_test=[x1_test;x2_test];

classify=zeros((n1+n2),1); %理论值
for i=1:(n1+n2)
g1(i) = (-0.5)*(xt(i,:)-u1')*invsig1*(xt(i,:)-u1')'-log(2*pi)-0.5*log(s1)+log(pw1);
g2(i) = (-0.5)*(xt(i,:)-u2')*invsig2*(xt(i,:)-u2')'-log(2*pi)-0.5*log(s2)+log(pw2);
if(g1(i)>g2(i))
classify(i)=1;
elseif(g1(i)classify(i)=2;
end
end
ne21=0;%应该是1类，错分到2类的数目
ne12=0;%应该是2类，错分到1类的数目
for i=1:n1
if(classify(i)~=1)
ne21=ne21+1;
end
end
for i=n1+1:n1+n2
if(classify(i)~=2)
ne12=ne12+1;
end
end
pe=(ne21+ne12)/(n1+n2)
us1=mean(x1) %参数估计
us2=mean(x2)
si1=zeros(2);
for i=1:n1
si1=si1+(x1(i,:)-us1)'*(x1(i,:)-us1);
end
ds1=si1/n1
si2=zeros(2);
for i=1:n2
si2=si2+(x2(i,:)-us1)'*(x2(i,:)-us2);
end
ds2=si2/n2
invsigs1=inv(ds1);
invsigs2=inv(ds2);
ss1=abs(det(ds1));
ss2=abs(det(ds2));

x=xt(:,1);
y=xt(:,2);
[X,Y]=meshgrid(x,y);
g=zeros((n1+n2),2);
sclassify=zeros((n1+n2),1); %训练样本值
for i=1:(n1+n2)
g1(i) = (-0.5)*(xt(i,:)-us1)*invsigs1*(xt(i,:)-us1)'-log(2*pi)-0.5*log(ss1)+log(pw1);
g2(i) = (-0.5)*(xt(i,:)-us2)*invsigs2*(xt(i,:)-us2)'-log(2*pi)-0.5*log(ss2)+log(pw2);
if(g1(i)>g2(i))
sclassify(i)=1;
elseif(g1(i)sclassify(i)=2;
end
g(i,:)=[g1(i) g2(i)];
end
sne21=0;%应该是1类，错分到2类的数目
sne12=0;%应该是2类，错分到1类的数目
for i=1:n1
if(sclassify(i)~=1)
sne21=sne21+1;
end
end
for i=n1+1:n1+n2
if(sclassify(i)~=2)
sne12=sne12+1;
end
end
spe=(sne21+sne12)/(n1+n2)
tclassify=zeros((n1+n2),1); %测试理论值
for i=1:40
g1(i) = (-0.5)*(x_test(i,:)-u1')*invsig1*(x_test(i,:)-u1')'-log(2*pi)-0.5*log(s1)+log(pw1);
g2(i) = (-0.5)*(x_test(i,:)-u2')*invsig2*(x_test(i,:)-u2')'-log(2*pi)-0.5*log(s2)+log(pw2);
if(g1(i)>g2(i))
tclassify(i)=1;
elseif(g1(i)tclassify(i)=2;
end
end
tne21=0;%应该是1类，错分到2类的数目
tne12=0;%应该是2类，错分到1类的数目
for i=1:20
if(classify(i)~=1)
tne21=tne21+1;
end
end
for i=21:40
if(classify(i)~=2)
tne12=tne12+1;
end
end
tpe=(tne21+tne12)/40
test_classify=zeros((n1+n2),1); %测试样本值
for i=1:40
g1(i) = (-0.5)*(x_test(i,

:)-us1)*invsigs1*(x_test(i,:)-us1)'-log(2*pi)-0.5*log(ss1)+log(pw1);
g2(i) = (-0.5)*(x_test(i,:)-us2)*invsigs2*(x_test(i,:)-us2)'-log(2*pi)-0.5*log(ss2)+log(pw2);
if(g1(i)>g2(i))
test_classify(i)=1;
elseif(g1(i)test_classify(i)=2;
end
end
teste21=0;%应该是1类，错分到2类的数目
teste12=0;%应该是2类，错分到1类的数目
for i=1:20
if(test_classify(i)~=1)
teste21=sne21+1;
end
end
for i=21:40
if(test_classify(i)~=2)
sne12=sne12+1;
end
end
test_pe=(sne21+sne12)/40
subplot(2,2,1)
syms xx yy;
gg = (-0.5)*[xx-u1(1),yy-u1(2)]*invsig1*[xx-u1(1),yy-u1(2)]'-0.5*log(s1)+log(pw1)-((-0.5)*[xx-u2(1),yy-u2(2)]*invsig2*[xx-u2(1),yy-u2(2)]'-0.5*log(s2)+log(pw2));
ezplot(gg)
hold on
for i=1:n1+n2
if(i<=n1)
if(classify(i)~=1)
plot(x1(i,:)-1,'r*');
else
plot(x1(i,:)-1,'ro');
end
else
if(classify(i)~=2)
plot(x2(i-n1,1)+3,x2(i-n1,2),'b*'); %为了让1类和2类在图像上区分开
else
plot(x2(i-n1,:)+3,x2(i-n1,2),'bo');
end
end
end
title('理论决策面')
subplot(2,2,2)
syms xs ys;
ggs = (-0.5)*[xs-us1(1),ys-us1(2)]*invsigs1*[xs-us1(1),ys-us1(2)]'-0.5*log(ss1)+log(pw1)-((-0.5)*[xs-us2(1),ys-us2(2)]*invsigs2*[xs-us2(1),ys-us2(2)]'-0.5*log(ss2)+log(pw2));
ezplot(ggs)
hold on
for i=1:n1+n2
if(i<=n1)
if(sclassify(i)~=1)
plot(x1(i,:)-1,'r*');
else
plot(x1(i,:)-1,'ro');
end
else
if(sclassify(i)~=2)
plot(x2(i-n1,1)+3,x2(i-n1,2),'b*'); %为了让1类和2类在图像上区分开
else
plot(x2(i-n1,1)+3,x2(i-n1,2),'bo');
end
end
end
title('样本决策面')
subplot(2,2,3)
syms xx yy;
gg = (-0.5)*[xx-u1(1),yy-u1(2)]*invsig1*[xx-u1(1),yy-u1(2)]'-0.5*log(s1)+log(pw1)-((-0.5)*[xx-u2(1),yy-u2(2)]*invsig2*[xx-u2(1),yy-u2(2)]'-0.5*log(s2)+log(pw2));
ezplot(gg)
hold on
for i=1:40
if(i<=20)
if(tclassify(i)~=1)
plot(x1_test(i,:)-1,'r*');
else
plot(x1_test(i,:)-1,'ro');
end
else
if(tclassify(i)~=2)
plot(x2_test(i-20,1)+3,x2_test(i-20,2),'b*'); %为了让1类和2类在图像上区分开
else
plot(x2_test(i-20,1)+3,x2_test(i-20,2),'bo');
end
end
end
title('测试理论决策面')

subplot(2,2,4)
syms xs ys;
ggs = (-0.5)*[xs-us1(1),ys-us1(2)]*invsigs1*[xs-us1(1),ys-us1(2)]'-0.5*log(ss1)+log(pw1)-((-0.5)*[xs-us2(1),ys-us2(2)]*invsigs2*[xs-us2(1),ys-us2(2)]'-0.5*log(ss2)+log(pw2));
ezplot(ggs)
hold on
for i=1:40
if(i<=20)
if(test_classify(i)~=1)
plot(x1_test(i,:)-1,'r*');
else
plot(x1_test(i,:)-1,'ro');
end
else
if(test_classify(i)~=2)
plot(x2_test(i-20,1)+3,x2_test(i-20,2),'b*'); %为了让1类和2类在

图像上区分开
else
plot(x2_test(i-20,1)+3,x2_test(i-20,2),'bo');
end
end
end
title('测试新决策面')