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

得样本数据如表3.1所示。

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

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

表3.1 某地区农村2周岁男婴的体格测量数据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 60.6 16.5;76 58.1 12.5;92 63.2 14.5;81 59.0 14.0;81 60.8 15.5;84 59.5 14.0]; [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 =82.0000 60.2000 14.5000 S =31.6000 8.0400 0.5000 8.0400 3.1720 1.3100 0.5000 1.3100 1.900 然后u=[90;58;16];t2=n*(xjunzhi-u)'*(S^(-1))*(xjunzhi-u) f=((n-p)/(p*(n-1)))*t2 输出结果t2 = 420.4447 f =84.0889所以21()'()T n X S X μμ-=--=420.44472(1)n p F T p n -=-=84.0889查表得F 3,3(0.05)=9.28<84.0889 F 3,3(0.01)=29.5<84.0889 因此在a=0.05或 a=0.01时拒绝0H 假设3.2 相应于表3.1再给出该地区9名2周岁女婴的三项指标的测量数据如表3.2所示。

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

试检验2周岁男婴与女婴的均值是有无显著差异表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 =82.000060.200014.5000Sx =31.6000 8.0400 0.50008.0400 3.1720 1.31000.5000 1.3100 1.900类似程序xjunzhi=[82;60.2;14.5];Sx=[31.6 8.04 0.5;8.04 3.1720 1.3100;0.5 1.31 1.9];n=6;y=[80.0 58.4 14.0;75.0 59.2 15;78 60.3 15;75.0 57.4 13.0;79 59.5 14.0;78 58.1 14.5;75 58.0 12.5;64 55.5 11.0;80 59.2 12.5];[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 =76.000058.400013.5000S =27.2308 6.5615 2.84626.5615 2.4323 1.40002.8462 1.4000 1.8462然后t=((n*m)/(n+m))*((xjunzhi-yjunzhi)')*(S^(-1))*(xjunzhi-yjunzhi)F=((n+m-p-1)/(p*(n+m-2)))*t输出结果t =5.3117F =1.4982查表得F0.05(3，11)=3.59>1.4982 F0.01(3，11)=6.22>1.4982H假设因此在a=0.05或a=0.01时接受第四章习题4.1 下表列举某年级任取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,'b-.p','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中并命名为mypolar.m然后输入程序x=[0:pi/2.5: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');输出结果第五章聚类分析习题5.3.下表给出我国历年职工人数（单位：万人），请用有序样品的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=000000000000000 0.022567000000000000000.448980.2457800000000000002.0632 1.39810.600240000000000003.9256 2.651 1.18020.11098000000000004.5022 3.0091 1.42380.569530.4086200000000005.179 3.4353 1.66480.825760.538310.020440000000006.0823 4.021 1.976 1.0230.633430.127810.047757000000007.0311 4.6502 2.3255 1.23130.7550.263410.112750.01245600000008.3322 5.5762 2.9094 1.6045 1.05310.606190.338810.131220.06003200000010.3127.1034 4.0117 2.4126 1.7772 1.37930.923140.526640.315410.0994*******12.6968.9972 5.4422 3.5114 2.7548 2.3553 1.669 1.04570.654960.256320.03671000016.29111.9987.8688 5.5038 4.5686 4.1193 3.1032 2.1468 1.47070.771220.308580.1276200021.11716.12811.3218.42987.2316 6.6487 5.2116 3.8312 2.7793 1.68770.88810.460160.10709002822.16716.52812.97811.38610.5468.5596 6.627 5.0716 3.4539 2.1748 1.34430.598320.199510我们只看下三角所有元素，其它元素理解为空第二步我们计算损失函数矩阵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=000000000000000.02256700000000000000.448980.0225670000000000000.559960.133550.022567000000000001.01850.559960.133550.02256700000000001.27470.58040.153990.0430070.020440000000001.4720.687770.261360.150380.0430070.02044000000001.68030.823370.396960.166440.0554640.0328970.01245600000002.0535 1.16620.711620.285210.166440.0554640.0328970.0124560000002.8616 1.77970.922770.496360.265840.154860.0554640.0328970.012456000003.9604 1.9366 1.07970.653280.321920.203150.0921740.0554640.0328970.0124600005.9528 2.3621 1.4747 1.02020.593790.321920.203150.0921740.0554640.03290.0124560008.7188 2.9416 2.0437 1.18680.760370.429010.310240.199270.0921740.055460.0328970.01245600 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}第六章判别分析例6.6对全国30个省市自治区1994年影响各地区经济增长差异的制度变量x1—经济增长率，x2—非国有化水平，x3—开放度，x4—市场化程度作贝叶斯判别分析。