matlab的判别分析.docx
- 文档编号:27440171
- 上传时间:2023-07-01
- 格式:DOCX
- 页数:9
- 大小:19.09KB
matlab的判别分析.docx
《matlab的判别分析.docx》由会员分享,可在线阅读,更多相关《matlab的判别分析.docx(9页珍藏版)》请在冰豆网上搜索。
matlab的判别分析
广西某锰矿床已知两种不同锰矿石各项评价指标如下表所列。
现新发现湖润锰矿床,初步定为氧化锰矿石,试用两类判别定其属性。
类别
序号
氧
化
锰
矿
石
1
2
3
41.19
34.99
35.62
11.86
9.84
10.56
0.182
0.178
0.261
36.22
27.82
21.02
碳
酸
锰
矿
石
1
2
3
23.21
25.05
19.23
5.46
6.84
6.61
0.11
0.134
0.137
21.17
27.3
26.61
湖
润
锰
石
(待判)
1
2
3
4
26.93
25.47
27.38
28.98
12.66
10.25
10.38
10.98
0.152
0.132
0.120
0.111
30.20
33.46
31.20
31.21
Matlab执行代码:
g1=[41.1911.860.18236.22;34.999.840.17827.82;35.6210.560.26121.02];
g2=[23.215.460.1121.17;25.056.840.13427.3;19.236.610.13726.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')
symsx1;symsx2;symsx3;symsx4;
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);
fori=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为:
')
forI=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
forI1=1:
3
if(x1w1(I1)>=x1w2(I1))
zp1(I1)=1;
end
end
zg1=sum(zp1);
fprintf('依次剔除g2总体中1,2,3号样本的判别函数值x2w1,x2w2为:
')
forJ=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
forJ1=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.9312.660.15230.20;25.4710.250.13233.46;27.3810.380.12031.20;28.9810.980.11131.21]; [p,q]=size(yp); forj=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 fork=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('计算广义平方距离函数: ') symsbx;symsbx1;symsbx2;symsbx3;symsbx4; 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 ! ! \n') fprintf('\n贝叶斯判别的回带判别') fort=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为: ') forT=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为: ') forT1=1: 3 if(b1d1(T1)<=b1d2(T1)) b1zp(T1)=1; end end forV=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 forV1=1: 3 if(b2d1(V1)>=b2d2(V1)) b2zp(V1)=1; end end fprintf('由上比较贝叶斯判别时,交叉法确认的误判率为: ') byp=(6-(sum(b1zp)+sum(b2zp)))/6 fprintf('根据以上的贝叶斯判别法则,判别待判样品yp\n') forv=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') forw=1: p if(yp1(w)>=yp2(w)) fprintf('属于氧化锰矿石总体的待判样品序号为: %g\n',w) end end
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- matlab 判别分析