matlab的判别分析
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
广西某锰矿床已知两种不同锰矿石各项评价指标如下表所列。现新发现湖润锰矿床,初步
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)