Matlab 与多元统计分析胡云峰 安庆师范学院第三章习题对某地区的6名2周岁男婴的身高、胸围、上半臂进行测量。

得样本数据如表所示。

假设男婴的测量数据X （a ）（a=1，…，6）来自正态总体N 3(?,∑) 的随机样本。

根据以往的资料，该地区城市2周岁男婴的这三项的均值向量?0=（90,58,16）’，试检验该地区农村男婴与城市男婴是否有相同的均值向量。

解1．预备知识 ∑未知时均值向量的检验： H 0：?=?0 H 1：?≠?0H 0成立时122)(0,)(1)(1,)()'((1)))()'()(,1)(1)1(,)(1)P P X N n S W n n X n S X n X S X T p n n p T F P n p n pμμμμμ---∑--∑⎪⎩∴----=-----+∴--:::: 当2(,)(1)n p T F p n p p n α-≥--或者22T T α≥拒绝0H当2(,)(1)n p T F p n p p n α-<--或者22T T α<接受0H这里2(1)(, )p n T F p n p n pαα-=--2．根据预备知识用matlab 实现本例题 算样本协方差和均值程序x=[78 ;76 ;92 ;81 ;81 ;84 ]; [n,p]=size(x); i=1:1:n;xjunzhi=(1/n)*sum(x(i,:)); y=rand(p,n);for j=1:1:ny(:,j)= x(j,:)'-xjunzhi'; y=y; endA=zeros(p,p); for k=1:1:n;A=A+(y(:,k)*y(:,k)'); endxjunzhi=xjunzhi' S=((n-1)^(-1))*A 输出结果xjunzhi = S =然后u=[90;58;16];t2=n*(xjunzhi-u)'*(S^(-1))*(xjunzhi-u) f=((n-p)/(p*(n-1)))*t2 输出结果t2 = f =所以21()'()T n X S X μμ-=--=2(1)n p F T p n -=-=查表得F 3,3=< F 3,3=< 因此在a=或 a=时拒绝0H 假设相应于表再给出该地区9名2周岁女婴的三项指标的测量数据如表所示。

假设女婴的测量数据Y (a)(a=1,…，9)来自正态总体N 3(?,∑)的随机样本。

试检验2周岁男婴与女婴的均值是有无显着差异表 某地区农村2周岁女婴体格测量数据解1. 预备知识有共同未知协方差阵∑时012:H μμ= 112:H μμ≠在0H 成立的情况下且两样本独立1112)(0,)(2)(1)(1)(2,)(2))((2)))))()'()(,2)21(P X Y PX Y N n m S n S m S W n m n m n m S n m T P n m n mn m p p n ---⎧-∑⎪⎨⎪+-=-+-+-∑⎩'⎤⎤∴+--+--⎥⎥⎦⎦'⎤⎤=--⎥⎥⎦⎦⋅=--+-++--+∴X Y X Y X Y S X Y X Y S X Y :::2(,1)2)T F P n m p m +--+-:给定检验水平α，查F 分布表，使{}p F F αα>=，可确定出临界值αF ，再用样本值计算出F ，若F F α>，则否定0H ，否则接受0H 。

2．根据预备知识用matlab 实现本例题 由上一题知道 xjunzhi =Sx =类似程序xjunzhi=[82;;];Sx=[ ; ; ];n=6;y=[ ; 15;78 15; ;79 ;78 ;75 ;64 ;80 ];[m,p]=size(y);i=1:1:m;yjunzhi=(1/m)*sum(y(i,:));z=rand(p,m);for j=1:1:mz(:,j)= y(j,:)'-yjunzhi';z=z;endB=zeros(p,p);for k=1:1:m;B=B+(z(:,k)*z(:,k)');endSy=((m-1)^(-1))*B;yjunzhi=yjunzhi'S=(1/(n+m-2))*((n-1)*Sx+(m-1)*Sy)得到结果yjunzhi =S =然后t=((n*m)/(n+m))*((xjunzhi-yjunzhi)')*(S^(-1))*(xjunzhi-yjunzhi) F=((n+m-p-1)/(p*(n+m-2)))*t输出结果t =F =查表得(3，11)=> (3，11)=>H假设因此在a=或a=时接受第四章习题下表列举某年级任取12名学生的5门主课的期末考试成绩，试绘制学生序号为1、2、11、12的轮廓图、雷达图。

解我们只需要数据如下1999493100100299889699971176724367781285755034371 利用matlab画轮廓图程序x=1:5;y1=[99 94 93 100 100];y2=[99 88 96 99 97];y3=[76 72 43 67 78];y4=[85 75 50 34 37];plot(x,y1,'k-o','linewidth',1);hold on;plot(x,y2,'r--*','linewidth',2);hold on;plot(x,y3,'','linewidth',2);hold onplot(x,y4,'k--o','linewidth',2);xlabel('学科');ylabel('分数');legend('1','2','11','12');set(gca,'xtick',[1 2 3 4 5])set(gca,'xticklabel',{'政治','语文','外语','数学','物理'})输出结果学科分数2 利用matlab 画雷达图此图用matlab 画起来比较复杂 首先我们修改polar 函数在命令窗口输入edit polar 结果会出现polar 函数的程序 其中我们把 % plot spokesth = (1:6)*2*pi/12;cst = cos(th); snt = sin(th); cs = [-cst; cst]; sn = [-snt; snt];line(rmax*cs,rmax*sn,'linestyle',ls,'color',tc,'linewidth',1,... 'handlevisibility','off','parent',cax) 修改为% plot spokesth = (1:3)*2*pi/6;cst = cos(th); snt = sin(th); cs = [-cst; cst]; sn = [-snt; snt];line(rmax*cs,rmax*sn,'linestyle',ls,'color',tc,'linewidth',1,... 'handlevisibility','off','parent',cax) 再将后面的所有程序中的30改为72 然后另存为work 中并命名为然后输入程序 x=[0:pi/:2*pi];y1=[99 94 93 100 100 99];y2=[99 88 96 99 97 99];y3=[76 72 43 67 78 76];y4=[85 75 50 34 37 85];mypolar(x,y1,'b');hold on;mypolar(x,y2,'m');hold on;mypolar(x,y3,'g');hold on;mypolar(x,y4,'y')legend('1','2','11','12');输出结果第五章聚类分析习题.下表给出我国历年职工人数（单位：万人），请用有序样品的fisher法聚类。

解第一步数据标准化后计算直径D程序：X=[1580 23;1881 121;2423 554;4532 662;5044 925;3303 1012;3465 1136;...3939 1264;4170 1334;4792 1424;5610 1524;6007 1644;6860 1813;...7451 2048;8019 2425];stdr=std(X);[n,m]=size(X);X=X./stdr(ones(n,1),:);[n p]=size(X);D=zeros(n,n);for i=1:1:n;for j=1:1:n;if i<jt=i:1:j;xgjunzhi=(1/(j-i+1))*sum(X(t,:));y=zeros(1,j-i+1);for s=i:1:jy(s)=(X(s,:)-xgjunzhi)*(X(s,:)-xgjunzhi)';ends=i:1:j;D(i,j)=sum(y);elseD(i,j)=0;endendendD=D'输出结果矩阵太大，所以用excel处理了一下D=我们只看下三角所有元素，其它元素理解为空第二步我们计算损失函数矩阵L程序：%设计一个把样品分为两类的程序,以及对应最后一类分割点D=D';L=zeros(n-1,n-1);alp=zeros(n-1,n-1);for m=2:n;s=zeros(1,m-1);for j=2:ms(1,j-1)=D(1,j-1)+D(j,m);endL(m-1,1)=min(s(1,1:m-1));for j=1:m-1if L(m-1,1)==s(1,j);alp(m-1,1)=j+1;endendend%分为k类for k=3:n;for m=k:ns=zeros(1,m-k+1);for j=k:m;s(1,j-k+1)=L(j-2,k-2)+D(j,m);endL(m-1,k-1)=min(s(1,1:m-k+1));for j=1:m-k+1if L(m-1,k-1)==s(1,j);alp(m-1,k-1)=j+k-1;endendendend输出结果这里由于表太大，用excel处理一下L=alp=20000000000000 33000000000000 44400000000000 44550000000000 46666000000000 46666700000000 46668880000000 46688889000000 4688101010101000000 41010101010111111110000 410101011111112121212000 4111111111313131313131300 101113131313131314141414140 1012131415151515151515151515在这里解释一下这两个矩阵行表示分为k类，k从2到15；列表示样本数m，m从2到15 我们只看下三角所有元素，其它元素理解为空，接下来我们根据结果分析如果我们要把样品分为三类，则第一个分割点为11，然后第二个分割点为6得到第一类：{1952，1954,1956,1958,1960}第二类：{1962,1964,1966,1968,1970}第三类：{1972,1974,1976,1978,1980}第六章判别分析例对全国30个省市自治区1994年影响各地区经济增长差异的制度变量x1—经济增长率，解求均值及协方差的逆的估计值程序X1=[ ; ; ;...; ; ;...20 ; ;19 ;... 16 ; ]; X2=[ ; ; ;... ; ; ;... 11 ;18 ; ;... ; ; ;... 84 ; ; ;... ];X3=[ ; ; ]; [n p]=size(X1); [m p]=size(X2); i=1:1:n;x1junzhi=(1/n)*sum(X1(i,:)); j=1:1:m;x2junzhi=(1/m)*sum(X2(j,:)); S1=cov(X1); S2=cov(X2);sigamani=(((n-1)*S1+(m-1)*S2)/(n+m-2))^(-1) x1junzhi=x1junzhi' x2junzhi=x2junzhi' 输出结果 sigamani =x1junzhi =x2junzhi =接着计算判别函数 根据111ln ''1,22g gg g f q X g μμμ--=-∑+∑=11ln 1ln0.897942716ln 2ln 0.5232527q q =≈-=≈-112342123445.86550.08960.08490.0715 1.240629.13440.08970.14430.0008 1.0591f x x x x f x x x x =-+-++=-+-++按照判别原则，若12f f >，则属于第一组，若12f f <，则属于第二组 回判 程序A=sigamani*x1junzhi; B=sigamani*x2junzhi; C=zeros(27,2); C(:,1)=[1:1:27]; for i=1:1:11f1=X1(i,:)*; f2=X1(i,:)*; if f1>f2C(i,2)=1; elseC(i,2)=2; end endfor i=1:1:16f1=X2(i,:)*; f2=X2(i,:)*; if f1>f2C(i+11,2)=1; elseC(i+11,2)=2; end end C输出结果 C =1 12 13 14 15 16 17 18 19 1 10 2 11 1 12 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 2所以误判率为1100%27⨯≈%很小，所以判别有效最后对待判样品进行判别程序D=zeros(3,2);D(:,1)=[28:1:30];for j=1:1:3f1=X3(j,:)*;f2=X3(j,:)*;if f1>f2D(j,2)=1;elseD(j,2)=2;endendD输出结果D =28 129 230 2第七章主成分分析例对全国30个省市自治区经济发展基本情况的八项指标作主成分分析，原始数据如下:解用matlab实现主成分分析第一步在matlab输入原始数据在这里由于输入数据量较大，我们可以在matlab的workspace中点击“新建变量”选项，命名为“x的变量，然后把你在excel中打好的表格中的数据直接复制粘贴到该变量中接着我们将原始数据标准化程序stdr=std(x); %求各变量的标准差[n,m]=size(x);sddata=x./stdr(ones(n,1),:) %标准化变换输出结果sddata =第二步建立指标间的相关系数矩阵R在这里标准化之后的样本数据的相关系数矩阵与样本离差阵相等所以我们接着在命令窗口输入R=cov(sddata)输出结果R =第三步求R的特征向量程序[x,B]=eig(R)输出结果x =B =00000000000000000000000000000000000000000000000000000000在这里由于输出结果数据长度太大，无法在这里显示，所以用excel对上面的矩阵B做了一点小小的处理在矩阵B中对角线上的元素对应的是R的特征值，对应的矩阵列向量为其特征向量对结果分析从上表看，前三个特征值累计贡献率已达%，这说明前三个主成分基本包含了全部指标具有的信息，为此，我们取前三个特征值，并计算出相应的特征向量:对应特征向量u1u2u3因而前三个主成为第一个主成分F1=++++在第一个主成分的表达式中第一、二、三项指标的系数较大，这三个指标起主要作用，我们可以把第一主成分看成是由国内生产总值，固定生产投资和居民消费水平所刻画的反映经济发展状况的综合指标。