matlab与多元统计分析

m a t l a b与多元统计分

析

Company Document number：WTUT-WT88Y-W8BBGB-BWYTT-19998

Matlab 与多元统计分析

胡云峰 安庆师范学院

第三章习题

对某地区的6名2周岁男婴的身高、胸围、上半臂进行测量。得样本数据如表所示。假设男婴的测量数据X （a ）（a=1，…，6）来自正态总体N 3(,∑)

的

随机样本。根据以往的资料，该地区城市2周岁男婴的这三项的均值向量0=

（90,58,16）’，试检验该地区农村男婴与城市男婴是否有相同的均值向量。

表 某地区农村2周岁男婴的体格测量数据

解

1．预备知识 ∑未知时均值向量的检验： H 0：=0 H 1：≠0

H 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: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 = 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. 预备知识

有共同未知协方差阵∑时 在0H 成立的情况下且两样本独立

给定检验水平α，查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: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 =

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的轮廓图、雷达图。

表学生学习成绩

解我们只需要数据如下

1999493100100

29988969997

117672436778

128575503437

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,'','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中并命名为

然后输入程序

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

t=i:1:j;

xgjunzhi=(1/(j-i+1))*sum(X(t,:));

y=zeros(1,j-i+1);

for s=i:1:j

y(s)=(X(s,:)-xgjunzhi)*(X(s,:)-xgjunzhi)';

end

s=i:1:j;

D(i,j)=sum(y);

else

D(i,j)=0;

end

end

end

D=D'

输出结果矩阵太大，所以用excel处理了一下

D=

000000000000000 00000000000000 0000000000000

00000000000

0000000000

000000000

00000000

0000000

000000

00000

0000

000

00 280我们只看下三角所有元素，其它元素理解为空

第二步我们计算损失函数矩阵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:m

s(1,j-1)=D(1,j-1)+D(j,m);

end

L(m-1,1)=min(s(1,1:m-1));

for j=1:m-1

if L(m-1,1)==s(1,j);

alp(m-1,1)=j+1;

end

end

end

%分为k类

for k=3:n;

for m=k:n

s=zeros(1,m-k+1);

for j=k:m;

s(1,j-k+1)=L(j-2,k-2)+D(j,m);

end

L(m-1,k-1)=min(s(1,1:m-k+1));

for j=1:m-k+1

if L(m-1,k-1)==s(1,j);

alp(m-1,k-1)=j+k-1;

end

end

end

end

输出结果这里由于表太大，用excel处理一下

L=

00000000000000 0000000000000 000000000000

00000000000

0000000000

000000000

00000000

0000000

000000

00000

0000

000

00 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—经济增长率，x2—非国有化水平，x3—开放度，x4—市场化程度作贝叶斯判别分析。

解求均值及协方差的逆的估计值

程序

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 =

接着计算判别函数

根据111

ln ''1,22

g g g

g f q X g μμμ--=-∑+∑=

按照判别原则，若12f f >，则属于第一组，若12f f <，则属于第二组 回判 程序

A=sigamani*x1junzhi; B=sigamani*x2junzhi; C=zeros(27,2); C(:,1)=[1:1:27]; for i=1:1:11 f1=X1(i,:)*; f2=X1(i,:)*; if f1>f2 C(i,2)=1; else

C(i,2)=2; end end

for i=1:1:16 f1=X2(i,:)*; f2=X2(i,:)*; if f1>f2

C(i+11,2)=1; else

C(i+11,2)=2; end end C 输出结果 C =

1 1

2 1

3 1

4 1

5 1

7 1

8 1

9 1

10 2

11 1

12 2

13 2

14 2

15 2

16 2

17 2

18 2

19 2

20 2

21 2

22 2

23 2

24 2

25 2

26 2

27 2

所以误判率为1

100%

27

?≈%很小，所以判别有效

最后对待判样品进行判别程序

D=zeros(3,2);

D(:,1)=[28:1:30];

for j=1:1:3

f1=X3(j,:)*;

f2=X3(j,:)*;

if f1>f2

D(j,2)=1;

else

D(j,2)=2;

end

end

D

输出结果

D =

28 1

30 2

第七章主成分分析

例对全国30个省市自治区经济发展基本情况的八项指标作主成分分析，原始数据如下:

表 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 =

0000000 0000000

0000000

0000000

0000000

0000000

0000000

0000000

在这里由于输出结果数据长度太大，无法在这里显示，所以用excel对上面的矩阵B做了一点小小的处理

在矩阵B中对角线上的元素对应的是R的特征值，对应的矩阵列向量为其特征向量

对结果分析如下表

从上表看，前三个特征值累计贡献率已达%，这说明前三个主成分基本包含了全部指标具有的信息，为此，我们取前三个特征值，并计算出相应的特征向量: 对应特征向量

u1u2u3

因而前三个主成为

第一个主成分

F1=++++在第一个主成分的表达式中第一、二、三项指标的系数较大，这三个指标起主要作用，我们可以把第一主成分看成是由国内生产总值，固定生产投资和居民消费水平所刻画的反映经济发展状况的综合指标。

在第二主成分中的表达式中第四、五、七、八项指标的系数较大·尤其四、五、七，我们把这四项指标看成一个综合指标。

在第三个主成分中，第六项指数影响最大，远远超出其它指标的影响，可单独看成是

居民消费价格指数的影响。

第八章因子分析

习题下表列出邓阜仙岩石的部分化学成分，试用此组数据作因子分析。

解利用matlab做因子分析

第一步将原始数据标准化

x=[ ;

;

;

;

;

;

];

stdr=std(x); %求各变量的标准差

[n,m]=size(x);

sddata=x./stdr(ones(n,1),:) %标准化变换