matlab与多元统计讲解

Matlab与多元统计 胡云峰  安庆师范学院 第三章习题 3.1对某地区的6名2周岁男婴的身高、胸围、上半臂进行测量。得样本数据如3.1所示。假设男婴的测量数据X（a）（a=1，…，6）来自正态总体N3(,∑)    的随机样本。根据以往的资料，该地区城市2周岁男婴的这三项的均值向量0=（90,58,16）’，试检验该地区农村男婴与城市男婴是否有相同的均值向量。 表3.1  某地区农村2周岁男婴的体格测量数据 男婴 身高（X1）cm 胸围身高(X2)cm 上半臂围身高（X3）cm 1 78 60.6 16.5 2 76 58.1 12.5 3 92 63.2 14.5 4 81 59 14 5 81 60.8 15.5 6 84 59.5 14         解 1．预备知识  ∑未知时均值向量的检验： H0：=0     H1：≠0 H0成立时  当 或者 拒绝 当 或者 接受 这里 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:n y(:,j)= x(j,:)'-xjunzhi'; y=y; end A=zeros(p,p); for k=1:1:n; A=A+(y(:,k)*y(:,k)'); end xjunzhi=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 所以 =420.4447 =84.0889 查表得F3,3(0.05)=9.28<84.0889      F3,3(0.01)=29.5<84.0889 因此在a=0.05或 a=0.01时拒绝 假设 3.2 相应于表3.1再给出该地区9名2周岁女婴的三项指标的测量数据如表3.2所示。假设女婴的测量数据Y(a)(a=1,…，9)来自正态总体N3(,∑)的随机样本。试检验2周岁男婴与女婴的均值是有无显著差异 表3.2 某地区农村2周岁女婴体格测量数据 女婴 身高（X1）cm 胸围身高(X2)cm 上半臂围身高（X3）cm 1 80 58.4 14 2 75 59.2 15 3 78 60.3 15 4 75 57.4 13 5 79 59.5 14 6 78 58.1 14.5 7 75 58 12.5 8 64 55.5 11 9 80 59.2 12.5         解 1. 预备知识 有共同未知协方差阵 时 在 成立的情况下且两样本独立 给定检验水平 ，查 分布表，使 ，可确定出临界值 ，再用样本值计算出 ，若 ，则否定 ，否则接受 。 2．根据预备知识用matlab实现本例题 由上一题知道 xjunzhi = 82.0000 60.2000 14.5000 Sx = 31.6000    8.0400    0.5000 8.0400    3.1720    1.3100 0.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:m z(:,j)= y(j,:)'-yjunzhi'; z=z; end B=zeros(p,p); for k=1:1:m; B=B+(z(:,k)*z(:,k)'); end Sy=((m-1)^(-1))*B; yjunzhi=yjunzhi' S=(1/(n+m-2))*((n-1)*Sx+(m-1)*Sy) 得到结果yjunzhi = 76.0000 58.4000 13.5000 S = 27.2308    6.5615    2.8462 6.5615    2.4323    1.4000 2.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.3117 F =1.4982 查表得F0.05(3，11)=3.59>1.4982  F0.01(3，11)=6.22>1.4982 因此在a=0.05或 a=0.01时接受 假设 第四章习题 4.1 下表列举某年级任取12名学生的5门主课的期末考试成绩，试绘制学生序号为1、2、11、12的轮廓图、雷达图。 表4.1 学生学习成绩 序号 政治 语文 外语 数学 物理 1 99 94 93 100 100 2 99 88 96 99 97 3 100 98 81 96 100 4 93 88 88 99 96 5 100 91 72 96 78 6 90 78 82 75 97 7 75 73 88 97 89 8 93 84 83 68 88 9 87 73 60 76 84 10 95 82 90 62 39 11 76 72 43 67 78 12 85 75 50 34 37             解  我们只需要数据如下 1 99 94 93 100 100 2 99 88 96 99 97 11 76 72 43 67 78 12 85 75 50 34 37             1 利用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 on plot(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 spokes th = (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 spokes th = (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法聚类。 年份 全民所有制 集体所有制 1952 1580 23 1954 1881 121 1956 2423 554 1958 4532 662 1960 5044 925 1962 3303 1012 1964 3465 1136 1966 3939 1264 1968 4170 1334 1970 4792 1424 1972 5610 1524 1974 6007 1644 1976 6860 1813 1978 7451 2048 1980 8019 2425       解 第一步数据标准化后计算直径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); 继续阅读