判别分析.docx
- 文档编号:30319704
- 上传时间:2023-08-13
- 格式:DOCX
- 页数:17
- 大小:17.83KB
判别分析.docx
《判别分析.docx》由会员分享,可在线阅读,更多相关《判别分析.docx(17页珍藏版)》请在冰豆网上搜索。
判别分析
§5.1
例5.1.1
a=[13.852.797.849.6;22.314.6712.3147.8;28.824.6316.1862.15;15.293.547.543.2;28.794.916.1258.1]
b=[2.181.061.2220.6;3.850.84.0647.1;11.403.50;3.662.422.1415.1;12.105.680]
两组间的平方距离(马氏距离),即一组的平均向量到另一组的平均向量的距离。
s=1/(10-2)*(cov(a,1)+cov(b,1))*5
(mean(a)-mean(b))*inv(s)*(mean(a)-mean(b))'
ans=
.028*********
检验两均值向量是否相等
组内离差阵:
(cov(a,1)+cov(b,1))*5
总的离差阵:
cov([a;b],1)*10
Wilks统计量为:
det((cov(a,1)+cov(b,1))*5)/det(cov([a;b],1)*10)
ans=
.0795********
F=((10-2)-4+1)/4*(1-ans)/ans
F=14.46435971681518
p=1-fcdf(ans,4,(10-2)-4+1)
p=
0.00589110677703
c=[8.853.385.1726.1;28.62.41.2127;20.76.77.630.2;7.92.44.333.2;3.193.321.439.9;12.45.14.4324.6;16.83.42.3131.3;152.75.0264]
d=zeros(8,1);
fori=1:
8;
d(i)=(c(i,:
)-0.5*(mean(a)+mean(b)))*inv(s)*(mean(a)-mean(b))'
end
d=
1.0e+002*
-0.06411600273506
1.90455683456686
0.45812985284833
.024*********
-0.131********092
0.22903978009930
0.59717160679228
0.48691733202362
find(d>0)
ans=
2
3
6
7
8
这几个为含钾盐
find(d<0)
ans=
1
4
5
symsx1x2x3x4
([x1x2x3x4]-0.5*(mean(a)+mean(b)))*inv(s)*(mean(a)-mean(b))'
Ans=
17088622476233187053/3602879701896396800*x1-1336113180613332608267383/36028797018963968000000+75513315568051364907/18014398509481984000*x2-4835311064286396129/562949953421312000*x3+26138143998347875787/36028797018963968000*x4
与书上的结果不一样。
§5.2
a=[2281342011;2451341040;2001671227;17015078;1001672014]
b=[225125714;130100612;15011776;1201331026;160100510]
c=[185115519;17012564;16514253;135108212;10011772]
d21=(mean(b)-mean(a))*inv(cov(a))*(mean(b)-mean(a))'
d21=
22.1219
d12=(mean(a)-mean(b))*inv(cov(b))*(mean(a)-mean(b))'
d12=
486.0310
D21=log(det(cov(a)))+d21
D21=
43.0647
D12=log(det(cov(b)))+d12
D12=
498.2681
D11=log(det(cov(a)))
D11=
20.9428
abc=[a;b;c]
a1=zeros(15,1);
fori=1:
15
a1(i)=(abc(i,:
)-mean(a))*inv(cov(a))*(abc(i,:
)-mean(a))'+log(det(cov(a)))
end
a2=zeros(15,1);
fori=1:
15
a2(i)=(abc(i,:
)-mean(b))*inv(cov(b))*(abc(i,:
)-mean(b))'+log(det(cov(b)))
end
a3=zeros(15,1);
fori=1:
15
a3(i)=(abc(i,:
)-mean(c))*inv(cov(c))*(abc(i,:
)-mean(c))'+log(det(cov(c)))
end
a123=[a1a2a3]
a123=
1.0e+003*
0.02416.59930.1087
0.02410.12480.1070
0.02410.03010.1651
0.02410.49000.0314
0.02411.90620.2650
0.02610.01540.0241
0.06960.01540.0197
0.04190.01540.0180
0.04040.01540.1223
0.05870.01540.0270
0.03890.14960.0185
0.03320.10890.0185
0.02790.82860.0185
0.06200.91940.0185
0.05620.04520.0185
找每一行最小值。
sa=zeros(15,1);
fori=1:
15
sa(i)=find(a123(i,:
)==min(a123(i,:
)))
end
sa=
1
1
1
1
1
2
2
2
2
2
3
3
3
3
3
§5.3
例5.3.2
a=[2281342011;2451341040;2001671227;17015078;1001672014]
b=[225125714;130100612;15011776;1201331026;160100510]
c=[185115519;17012564;16514253;135108212;10011772]
A=5*(cov(a,1)+cov(b,1)+cov(c,1))%组内离差阵A
B=15*cov([a;b;c],1)-A%组间离差阵B
[ld]=eig(inv(A)*B)
l1=l(:
1)
a123=(repmat([a;b;c]*l1,1,3)-repmat([mean(a)*l1mean(b)*l1mean(c)*l1],15,1))*inv(sqrt(diag([l1'*cov(a)*l1l1'*cov(b)*l1l1'*cov(c)*l1])))
sa=zeros(15,1);
fori=1:
15
sa(i)=find(abs(a123(i,:
))==min(abs(a123(i,:
))))
end
sa=
1
1
1
2
1
2
2
3
2
3
2
2
2
2
3
判对率为:
sum((sa-[ones(5,1);2*ones(5,1);3*ones(5,1)])==0)/15
ans=0.5333
§5.4判别效果的检验与各变量判别能力的检验
二、k个总体判别效果的检验(k>2)
a=[2281342011;2451341040;2001671227;17015078;1001672014]
b=[225125714;130100612;15011776;1201331026;160100510]
c=[185115519;17012564;16514253;135108212;10011772]%三总体判别效果的检验,
A=5*(cov(a,1)+cov(b,1)+cov(c,1))%组内离差阵A
T=15*cov([a;b;c],1)%总的离差阵T
wi=det(A)/det(T)
p=4;n=15-3;pf=1-fcdf((n-p+1)*(1-sqrt(wi))/(p*sqrt(wi)),2*p,2*(n-p+1))%pf=0.0370,拒绝原假设。
ab=(mean(a)-mean(b))*inv((cov(a)+cov(b)+cov(c))*4/(15-3))*(mean(a)-mean(b))'%a与b的马氏距离
ac=(mean(a)-mean(c))*inv((cov(a)+cov(b)+cov(c))*4/(15-3))*(mean(a)-mean(c))'%a与c的马氏距离
bc=(mean(b)-mean(c))*inv((cov(a)+cov(b)+cov(c))*4/(15-3))*(mean(b)-mean(c))'%c与b的马氏距离
n=15;k=3;m=4;
pab=1-fcdf((n-k-m+1)*5*5*ab/(m*(n-k)*(5+5)),m,n-k-m+1);
pac=1-fcdf((n-k-m+1)*5*5*ac/(m*(n-k)*(5+5)),m,n-k-m+1);
pbc=1-fcdf((n-k-m+1)*5*5*bc/(m*(n-k)*(5+5)),m,n-k-m+1);
[pabpbcpac]
ans=
0.02440.72830.0144
结论与p199相同,即第一类与二三类之间有显著差异,二三类之间没有显著差异。
三、各变量判别能力的检验
a=[2281342011;2451341040;2001671227;17015078;1001672014]
b=[225125714;130100612;15011776;1201331026;160100510]
c=[185115519;17012564;16514253;135108212;10011772]
A=5*(cov(a,1)+cov(b,1)+cov(c,1))%组内离差阵A
T=15*cov([a;b;c],1)%总的离差阵T
A1=zeros(4,4);
fori=1:
4,j=1:
4
A1(i,j)=A(i,j)-A(1,j)*A(i,1)/A(1,1)
end
A1(i,1)=-A(i,1)/A(1,1);A1(1,j)=A(1,j)/A(1,1);A1(1,1)=1/A(1,1);
A2=zeros(4,4);
fori=1:
4,j=1:
4
A2(i,j)=A1(i,j)-A1(2,j)*A1(i,2)/A1(2,2)
end
A2(i,2)=-A1(i,2)/A1(2,2);A2(2,j)=A1(2,j)/A1(2,2);A2(2,2)=1/A1(2,2);
A3=zeros(4,4);
fori=1:
4,j=1:
4
A3(i,j)=A2(i,j)-A2(3,j)*A2(i,3)/A2(3,3)
end
A3(i,3)=-A2(i,3)/A2(3,3);A3(3,j)=A2(3,j)/A2(3,3);A3(3,3)=1/A2(3,3);
T1=zeros(4,4);
fori=1:
4,j=1:
4
T1(i,j)=T(i,j)-T(1,j)*T(i,1)/T(1,1)
end
T1(i,1)=-T(i,1)/T(1,1);T1(1,j)=T(1,j)/T(1,1);T1(1,1)=1/T(1,1);
T2=zeros(4,4);
fori=1:
4,j=1:
4
T2(i,j)=T1(i,j)-T1(2,j)*T1(i,2)/T1(2,2)
end
T2(i,2)=-T1(i,2)/T1(2,2);T2(2,j)=T1(2,j)/T1(2,2);T2(2,2)=1/T1(2,2);
T3=zeros(4,4);
fori=1:
4,j=1:
4
T3(i,j)=T2(i,j)-T2(3,j)*T2(i,3)/T2(3,3)
end
T3(i,3)=-T2(i,3)/T2(3,3);T3(3,j)=T2(3,j)/T2(3,3);T3(3,3)=1/T2(3,3);
U4=A3(4,4)/T3(4,4);
n=15;k=3;r=3;p4=1-fcdf((n-k-r)*(1-U4)/((k-1)*U4),k-1,n-k-r);
p4=0.5768,去掉x4
先去掉x4,再看是否要去掉x1
a=[2281342011;2451341040;2001671227;17015078;1001672014]
b=[225125714;130100612;15011776;1201331026;160100510]
c=[185115519;17012564;16514253;135108212;10011772]
a1=a(:
[1:
3]);b1=b(:
[1:
3]);c1=c(:
[1:
3])
AA=5*(cov(a1,1)+cov(b1,1)+cov(c1,1))%组内离差阵AA
TT=15*cov([a1;b1;c1],1)%总的离差阵TT
AA1=zeros(3,3);
fori=1:
3,j=1:
3
AA1(i,j)=AA(i,j)-AA(2,j)*AA(i,2)/AA(2,2)
end
AA1(i,2)=-AA(i,2)/AA(2,2);AA1(2,j)=AA(2,j)/AA(2,2);AA1(2,2)=1/AA(2,2);
AA2=zeros(3,3);
fori=1:
3,j=1:
3
AA2(i,j)=AA1(i,j)-AA1(3,j)*AA1(i,3)/AA1(3,3)
end
AA2(i,3)=-AA1(i,3)/AA1(3,3);AA2(3,j)=AA1(3,j)/AA1(3,3);AA2(3,3)=1/AA1(3,3);
TT1=zeros(3,3);
fori=1:
3,j=1:
3
TT1(i,j)=TT(i,j)-TT(2,j)*TT(i,2)/TT(2,2)
end
TT1(i,2)=-TT(i,2)/TT(2,2);TT1(2,j)=TT(2,j)/TT(2,2);TT1(2,2)=1/TT(2,2);
TT2=zeros(3,3);
fori=1:
3,j=1:
3
TT2(i,j)=TT1(i,j)-TT1(3,j)*TT1(i,3)/TT1(3,3)
end
TT2(i,3)=-TT1(i,3)/TT1(3,3);TT2(3,j)=TT1(3,j)/TT1(3,3);TT2(3,3)=1/TT1(3,3);
U1=AA2(1,1)/TT2(1,1);
n=15;k=3;r=2;p1=1-fcdf((n-k-r)*(1-U1)/((k-1)*U1),k-1,n-k-r)
p1=
0.2979
去掉x1
§5.5逐步判别
例5.5.1
a=[2281342011;2451341040;2001671227;17015078;1001672014]
b=[225125714;130100612;15011776;1201331026;160100510]
c=[185115519;17012564;16514253;135108212;10011772]%看能否引入变量
A=5*(cov(a,1)+cov(b,1)+cov(c,1))%组内离差阵A
T=15*cov([a;b;c],1)%总的离差阵T
U1=min(diag(A./T));n=15;k=3;
P1=1-fcdf((1-U1)*(n-k)/(U1*(k-1)),k-1,n-k)
P1=
0.0060
find(diag(A./T)==min(diag(A./T)))
ans=
2
引入x2变量,并对A、T作相应的消去变换
A1=zeros(4,4);
fori=1:
4,j=1:
4
A1(i,j)=A(i,j)-A(2,j)*A(i,2)/A(2,2)
end
A1(i,2)=-A(i,2)/A(2,2);A1(2,j)=A(2,j)/A(2,2);A1(2,2)=1/A(2,2);
T1=zeros(4,4);
fori=1:
4,j=1:
4
T1(i,j)=T(i,j)-T(2,j)*T(i,2)/T(2,2)
end
T1(i,2)=-T(i,2)/T(2,2);T1(2,j)=T(2,j)/T(2,2);T1(2,2)=1/T(2,2);
AT=diag(A1./T1)
AT=
0.8221
2.3457
0.6917
0.8175
看是否引入x3
n=15;k=3;r=1;P3=1-fcdf((1-AT(3))*(n-k-r)/(AT(3)*(k-1)),k-1,n-k-r)
P3=
0.1317
在显著性水平0.15下,拒绝原假设。
引入x3
对A1、T1做消去变换
T2=zeros(4,4);
fori=1:
4,j=1:
4
T2(i,j)=T1(i,j)-T1(3,j)*T1(i,3)/T1(3,3)
end
T2(i,3)=-T1(i,3)/T1(3,3);T2(3,j)=T1(3,j)/T1(3,3);T2(3,3)=1/T1(3,3);
A2=zeros(4,4);
fori=1:
4,j=1:
4
A2(i,j)=A1(i,j)-A1(3,j)*A1(i,3)/A1(3,3)
end
A2(i,3)=-A1(i,3)/A1(3,3);A2(3,j)=A1(3,j)/A1(3,3);A2(3,3)=1/A1(3,3)
看能否去掉变量P207
diag(A1./T1)
ans=
0.8221
2.3457
0.6917
0.8175
U32=0.6917,即对x2进行消去变换,
diag(T2./A2)
ans=
1.2741
0.4263
0.6917
1.2513
T3=zeros(4,4);
fori=1:
4,j=1:
4
T3(i,j)=T(i,j)-T(3,j)*T(i,3)/T(3,3)
end
T3(i,3)=-T(i,3)/T(3,3);T3(3,j)=T(3,j)/T(3,3);T3(3,3)=1/T(3,3);
A3=zeros(4,4);
fori=1:
4,j=1:
4
A3(i,j)=A(i,j)-A(3,j)*A(i,3)/A(3,3)
end
A3(i,3)=-A(i,3)/A(3,3);A3(3,j)=A(3,j)/A(3,3);A3(3,3)=1/A(3,3)
diag(A3./T3)
ans=
0.8079
0.6666
2.2607
0.7816
T4=zeros(4,4);
fori=1:
4,j=1:
4
T4(i,j)=T3(i,j)-T3(2,j)*T3(i,2)/T3(2,2)
end
T4(i,2)=-T3(i,2)/T3(2,2);T4(2,j)=T3(2,j)/T3(2,2);T4(2,2)=1/T3(2,2);
A4=zeros(4,4);
fori=1:
4,j=1:
4
A4(i,j)=A3(i,j)-A3(2,j)*A3(i,2)/A3(2,2)
end
A4(i,2)=-A3(i,2)/A3(2,2);A4(2,j)=A3(2,j)/A3(2,2);A4(2,2)=1/A3(2,2);
diag(T4./A4)
ans=
1.2741
0.6666
0.4423
1.2513
取两个红色值最大者
U=diag(T2./A2)
P=1-fcdf((1-U(3))*(n-r+1-k)/(U(3)*(k-1)),k-1,n-k-r+1)
P=
0.1095
0.15下,拒绝原假设。
不去掉x2、x3
看能否引入新变量
AT=diag(A2./T2)
n=15;k=3;r=2;P3=1-fcdf((1-AT(3))*(n-k-r)/(AT(3)*(k-1)),k-1,n-k-r)
P3=
1
不再接受变量,即最终入选x2、x3两变量
用距离法判断,并假定协方差相等。
a1=a(:
[23]);b1=b(:
[23]);c1=c(:
[23]);
S=4*(cov(a1)+cov(b1)+cov(c1))/(15-3);
abc=[a1;b1;c1]
z1=zeros(15,1);
fori=1:
15
z1(i)=(abc(i,:
)-mean(a1))*inv(S)*(abc(i,:
)-mean(a1))'
end
z2=zeros(15,1);
fori=1:
15
z2(i)=(abc(i,:
)-mean(b1))*inv(S)*(abc(i,:
)-mean(b1))'
end
z3=zeros(15,1);
fori=1:
15
z3(i)=(abc(i,:
)-mean(c1))*inv(S)*(abc(i,:
)-mean(c1))'
end
z123=[z1z2z3]
sa=zeros(15,1);
fori=1:
15
sa(i)=find(z123(i,:
)==min(z123(i,:
)))
end
sa=
1
1
1
1
1
3
2
2
2
2
3
3
3
3
2
与书上p211答案一模一样。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 判别分析
![提示](https://static.bdocx.com/images/bang_tan.gif)