matlab的判别分析

广西某锰矿床已知两种不同锰矿石各项评价指标如下表所列。现新发现湖润锰矿床，初步

Matlab执行代码:

g1=[41.19 11.86 0.182 36.22;34.99 9.84 0.178 27.82;35.62 10.56 0.261

21.02];

g2=[23.21 5.46 0.11 21.17;25.05 6.84 0.134 27.3;19.23 6.61 0.137 26.61]; fprintf('做距离判别分析:\n')

fprintf('在两个总体的协方差矩阵相等的假设下进行判别分析:\n')

fprintf('两个样本的协方差矩阵s1,s2分别为\n')

s1=cov(g1)

s2=cov(g2)

fprintf('因为两个总体的协方差矩阵相等,所以协方差的联合估计s为:\n') [m1,n2]=size(g1);[m2,n2]=size(g2);

s=((m1-1)*s1+(m2-1)*s2)/(m1+m2-2)

fprintf('两个总体的马氏平方距离为:\n')

sn=inv(s);

u1=mean(g1);u2=mean(g2);

ucz=(u1-u2)';

dmj=(u1-u2)*sn*ucz

fprintf('该值反映了两个总体的分离程度,线性函数的判别函数为:\n')

syms x1;syms x2;syms x3;syms x4;

x=[x1;x2;x3;x4];

u1z=u1';u2z=u2';

a1=(sn*u1z)';b1=(u1*sn*u1z)/2;

a2=(sn*u2z)';b2=(u2*sn*u2z)/2;

w1=vpa((a1*x-b1),4)

w2=vpa((a2*x-b2),4)

fprintf('用回代法作出误判率p1为:\n')

fprintf('比较gwh1和gwh2大小\n')

g=[g1;g2];

[m,n]=size(g);

for i=1:m

ghdw1(i,:)=a1.*g(i,:);

ghdw2(i,:)=a2.*g(i,:);

gwh1(i)=sum(ghdw1(i,:))-b1;

gwh2(i)=sum(ghdw2(i,:))-b2;

end

gwh1

gwh2

fprintf('经比较得g1中1,2,3号判入g1;g2中1,2,3号判入g2,则误判率的回带估计为:\n')

p1=0

fprintf('用交叉估计法确认距离判别的误判率：\n')

fprintf('依次剔除g1总体中1，2，3号样本是的判别函数值x1w1，x1w2为：')

for I=1:3

xg1=g1;

xg1(I,:)=[];

xs1=cov(xg1);

x1s=(xs1+2*s2)/3;x1sn=x1s';

xu1=mean(xg1);

x1w1(I)=sum((x1sn*xu1')'.*g1(I,:))-0.5*xu1*x1sn*xu1';

x1w2(I)=sum((x1sn*u2')'.*g1(I,:))-0.5*u2*x1sn*u2';

end

x1w1

x1w2

for I1=1:3

if(x1w1(I1)>=x1w2(I1))

zp1(I1)=1;

end

end

zg1=sum(zp1);

fprintf('依次剔除g2总体中1，2，3号样本的判别函数值x2w1,x2w2为:') for J=1:3

xg2=g2;

xg2(J,:)=[];

xs2=cov(xg2);

x2s=(2*s1+xs2)/3;x2sn=x2s';

xu2=mean(xg2);

x2w1(J)=sum((x2sn*xu2')'.*g2(J,:))-0.5*u1*x2sn*u1';

x2w2(J)=sum((x2sn*xu2')'.*g2(J,:))-0.5*xu2*x2sn*xu2';

end

x2w1

x2w2

for J1=1:3

if(x2w1(J1)

zp2(J1)=1;

end

end

zg2=sum(zp2);

fprintf('由上比较得，交叉法所得的误判率为：\n')

zp=zg1+zg2;

jwpl=(6-zp)/6

fprintf('判别新样品:\n')

yp=[26.93 12.66 0.152 30.20;25.47 10.25 0.132 33.46;27.38 10.38 0.120 31.20;28.98 10.98 0.111 31.21];

[p,q]=size(yp);

for j=1:p

w1p(j,:)=a1.*yp(j,:);

w2p(j,:)=a2.*yp(j,:);

w1pb(j)=sum(w1p(j,:))-b1;

w2pb(j)=sum(w2p(j,:))-b2;

end

w1pb

w2pb

for k=1:4

if(w1pb(k)>=w2pb(k))

fprintf('属于氧化锰矿石的样本序号是%g\n',k)

end

end

fprintf('用贝叶斯判别法分析:\n')

fprintf('\n在两个总体的协方差矩阵相等的假设下做贝叶斯判别:\n')

fprintf('\n先验概率按比例分配求得总体g1,g2的先验概率分别为:\n')

bp1=m1/(m1+m2)

bp2=m2/(m1+m2)

fprintf('两个正态总体的贝叶斯判别为:\n')

bw1=w1+log(bp1);

bw2=w2+log(bp2);

fprintf('当两个总体的协方差矩阵,误判损失相同且先验概率按比例分配时距离判别与贝叶斯判别等价\n')

fprintf('计算广义平方距离函数:')

syms bx;syms bx1;syms bx2;syms bx3;syms bx4;

bx=[bx1;bx2;bx3;bx4];

bdp1=vpa((bx-u1z)'*sn*(bx-u1z)-2*log(bp1),4)

bdp2=vpa((bx-u2z)'*sn*(bx-u2z)-2*log(bp2),4)

fprintf('后验概率pg1,pg2为:\n')

pg1=exp(-0.5*bdp1)/(exp(-0.5*bdp1)+exp(-0.5*bdp2))

pg2=exp(-0.5*bdp2)/(exp(-0.5*bdp1)+exp(-0.5*bdp2))

fprintf('此时贝叶斯判别法则为:当pg1>=pg2时,属于g1总体;当pg1

fprintf('\n贝叶斯判别的回带判别')

for t=1:m

bdg1(t)=(g(t,:)'-u1z)'*sn*(g(t,:)'-u1z)-2*log(bp1);

bdg2(t)=(g(t,:)'-u2z)'*sn*(g(t,:)'-u2z)-2*log(bp2);

p1b(t)=exp(-0.5*bdg1(t))/(exp(-0.5*bdg1(t))+exp(-0.5*bdg2(t)));

p2b(t)=exp(-0.5*bdg2(t))/(exp(-0.5*bdg1(t))+exp(-0.5*bdg2(t))); end

fprintf('回代g1,g2中六个样本,求得的后验概率为:\n')

p1b

p2b

fprintf('经比较得,误判率的回带估计bp为:\n')

bp=0

fprintf('贝叶斯判别的交叉法确认误判率:\n')

fprintf('依次踢除g1总体中1,2,3号样本,所得的广义平方距离b1d1,b1d2为:') for T=1:3

bxg1=g1;

bxg1(T,:)=[];

bju1=mean(bxg1);

b1s1=cov(bxg1);b1s=(b1s1+2*s2)/3;

bj1p1=2/5 ; bj1p2=3/5;

b1d1(T)=(g1(T,:)-bju1)*b1s'*(g1(T,:)'-bju1')-2*log(bj1p1);

b1d2(T)=(g1(T,:)-u2)*b1s'*(g1(T,:)'-u2')-2*log(bj1p2);

b1p1(T)=exp(-0.5*b1d1(T))/(exp(-0.5*b1d1(T))+exp(-0.5*b1d2(T)));

b1p2(T)=exp(-0.5*b1d2(T))/(exp(-0.5*b1d1(T))+exp(-0.5*b1d2(T))); end

b1d1

b1d2

fprintf('依次剔除g2总体中1,2,3号样本,所得的广义平方距离b2d1,b2d2为:') for T1=1:3

if(b1d1(T1)<=b1d2(T1))

b1zp(T1)=1;

end

end

for V=1:3

bxg2=g2;

bxg2(V,:)=[];

bju2=mean(bxg2);

b2s2=cov(bxg2);b2s=(2*s1+b2s2)/3;

bj2p1=3/5;bj2p2=2/5;

b2d1(V)=(g2(V,:)-u1)*b2s'*(g2(V,:)'-u1')-2*log(bj2p1);

b2d2(V)=(g2(V,:)-bju2)*b2s'*(g2(V,:)'-bju2')-2*log(bj2p2);

b2p1(V)=exp(-0.5*b2d1(V))/(exp(-0.5*b2d1(V))+exp(-0.5*b2d2(V)));

b2p2(V)=exp(-0.5*b2d2(V))/(exp(-0.5*b2d1(V))+exp(-0.5*b2d2(V))); end

b2d1

b2d2

for V1=1:3

if(b2d1(V1)>=b2d2(V1))

b2zp(V1)=1;

end

end

fprintf('由上比较贝叶斯判别时,交叉法确认的误判率为:')

byp=(6-(sum(b1zp)+sum(b2zp)))/6

fprintf('根据以上的贝叶斯判别法则,判别待判样品yp\n')

for v=1:p

ydg1(v)=(yp(v,:)'-u1z)'*sn*(yp(v,:)'-u1z)-2*log(bp1);

ydg2(v)=(yp(v,:)'-u2z)'*sn*(yp(v,:)'-u2z)-2*log(bp2);

yp1(v)=exp(-0.5*ydg1(v))/(exp(-0.5*ydg1(v))+exp(-0.5*ydg2(v)));

yp2(v)=exp(-0.5*ydg2(v))/(exp(-0.5*ydg1(v))+exp(-0.5*ydg2(v))); end

fprintf('后验概率yp1,yp2为:\n')

yp1

yp2

fprintf('比较后验概率yp1,yp2知:\n')

for w=1:p

if(yp1(w)>=yp2(w))

fprintf('属于氧化锰矿石总体的待判样品序号为:%g\n',w) end

end