matlab ssgs工具箱,第10章：多元分析

10.1 聚类分析

即群分析，是对多个样本(或指标)进行定量分类的一种多元统计分析方法。对样本进行分类称为Q型聚类分析，对指标进行分类称为R型聚类分析。

10.1.1 Q型聚类分析

(1)样本的相似性度量

对于定量变量，最常用的是闵式距离

绝对值距离

欧几里得距离：最常用，当坐标轴进行正交旋转时，它保持不变

切比雪夫距离

马氏距离：对一切线性变换是不变的

(2)类与类之间的相似性度量

最短距离法

最长距离法

重心法

类平均法

离差平方和法

(3)最短距离法(最近邻法)的计算步骤：

clc,clear

a=[1,0;1,1;3,2;4,3;2,5];

[m,n]=size(a);

d=zeros(m);

d=mandist(a'); %mandist求矩阵列向量组之间的两两绝对值距离

d=tril(d); %截取下三角元素

nd=nonzeros(d); %去掉d中的零元素，非零元素按列排列

nd=union([],nd) %去掉重复的非零元素

for i=1:m-1

nd_min=min(nd);

[row,col]=find(d==nd_min);tm=union(row,col); %row和col归为一类

tm=reshape(tm,1,length(tm)); %把数据tm变成行向量

fprintf('第%d次合成，平台高度为%d时的分类结果为：%s\n',...

i,nd_min,int2str(tm));

nd(nd==nd_min)=[]; %删除已经归类的元素

if length(nd)==0

break

end

end

使用Matlab工具箱

clc,clear

a=[1,0;1,1;3,2;4,3;2,5];

y=pdist(a,'cityblock'); %求a的两两行向量间的绝对值距离

yc=squareform(y) %变换成距离方阵

z=linkage(y) %产生等级聚类树

dendrogram(z) %画聚类图

T=cluster(z,'maxclust',3) %°把对象划分成3类

for i=1:3

tm=find(T==i); %求第i类的对象

tm=reshape(tm,1,length(tm)); %变成行向量

fprintf('µÚ%dÀàµÄÓÐ%s\n',i,int2str(tm)); %显示分类结果

end

(4)Matlab聚类分析的相关命令

pdist：计算矩阵中两两对象间的欧式距离

linkage：使用最短距离算法生成具有层次结构的聚类树

cluster：从连接输出中创建聚类

zsore：对数据矩阵进行标准化处理

dendrogram：画聚类树状图

clusterdata：数据分类

squareform：输出转换为方阵

10.1.2 R型聚类法

(1)变量相似性度量

相关系数

夹角余弦

(2)变量聚类法

最长距离法

最短距离法clc,clear

a=textread('ch.txt');

d=1-abs(a); %进行数据变换，把相关系数转化为距离

d=tril(d); %提出d矩阵的下三角部分

b=nonzeros(d);%去掉d中的零元素

b=b'; %化成行向量

z=linkage(b,'complete'); %按最长距离法聚类

y=cluster(z,'maxclust',2) %把变量划分成两类

ind1=find(y==1);ind1=ind1' %显示第一类对应的变量标号

ind2=find(y==2);ind2=ind2' %显示第二类对应的变量标号

h=dendrogram(z); %画聚类图

set(h,'Color','k','LineWidth',1.3) %把聚类图线的颜色改为黑色，线宽加粗

10.1.3 聚类分析案例

(1)R型聚类分析

clc, clear

a=load('gj.txt');

b=zscore(a); %数据标准化

r=corrcoef(b) %计算相关系数矩阵

%d=tril(1-r); d=nonzeros(d)'; %另外一种计算距离的方法

d=pdist(b','correlation'); %计算相关系数导出的距离

z=linkage(d,'average'); %按类平均法聚类

h=dendrogram(z); %画聚类图

set(h,'Color','k','LineWidth',1.3)

T=cluster(z,'maxclust',6) %把变量划分成6类

for i=1:6

tm=find(T==i); %求第i类的对象

tm=reshape(tm,1,length(tm)); %变成行向量

fprintf('µÚ%dÀàµÄÓÐ%s\n',i,int2str(tm)); %显示分类结果

end

(2)Q型聚类分析

clc,clear

load gj.txt

gj(:,[3:6])=[]; %删除数据矩阵的第3-6列，即使用变量1,2,7,8,9,10

gj=zscore(gj); %数据标准化

y=pdist(gj); %求对象间的欧式距离，每行是一个对象

z=linkage(y,'average'); %按类平均法聚类

h=dendrogram(z); %画聚类图

set(h,'Color','k','LineWidth',1.3)

for k=3:5

fprintf('划分成%d类的结果如下\n',k)

T=cluster(z,'maxclust',k);

for i=1:k

tm=find(T==i);

tm=reshape(tm,1,length(tm));

fprintf('第%d类的有%s\n',i,int2str(tm));

end

if k==5

break

end

fprintf('**********************************\n');

end

10.2 主成分分析

希望选用较少的变量去解释原来资料中的大部分变异，将我们手中许多相关性很高的变量转化成彼此相互独立或者不相关的变量。是一种降维方法

(1)主成分回归分析：

clc,clear

load sn.txt

[m,n]=size(sn);

x0=sn(:,[1:n-1]);y0=sn(:,n);

hg1=[ones(m,1),x0]\y0; %计算普通最下二乘法的回归系数

hg1=hg1' %变成行向量显示回归系数

fprintf('y=%f',hg1(1)); %开始显示普通最小二乘法回归结果

for i=2:n

if hg1(i)>0

fprintf('+%f*x%d',hg1(i),i-1);

else

fprintf('%f*x%d',hg1(i),i-1)

end

end

fprintf('\n')

r=corrcoef(x0) %计算相关系数矩阵

xd=zscore(x0); %对设计矩阵进行标准化处理

yd=zscore(y0); %对y0进行标准化处理

[vec1,lamda,rate]=pcacov(r) %vec1为r的特征向量，lamda为r的特征值，rate为各个主成分的贡献率

f=repmat(sign(sum(vec1)),size(vec1,1),1); %构造与vec1同维数的元素为矩阵

vec2=vec1.*f %修改特征向量的正负号，使得特征向量的所有分类和为正

contr=cumsum(rate) %计算贡献累积率，第i个分量表示前i个主成分的贡献率

df=xd*vec2; %计算所欲主成分的得分

num=input('请选项主成分的个数:') %通过累积贡献率交互式选择主成分的个数

hg21=df(:,[1:num])\yd %主成分变量的回归系数

hg22=vec2(:,1:num)*hg21 %标准化变量的回归方程系数

hg23=[mean(y0)-std(y0)*mean(x0)./std(x0)*hg22, std(y0)*hg22'./std(x0)] %计算原始变量的回归方程系数

fprintf('y=%f',hg23(1)); %开始显示主成分回归结果

for i=2:n

if hg23(i)>0

fprintf('+%f*x%d',hg23(i),i-1);

else

fprintf('%f*x%d',hg23(i),i-1);

end

end

fprintf('\n')

%下面机损两种回归分析的剩余标准差

rmse1=sqrt(sum((hg1(1)+x0*hg1(2:end)'-y0).^2)/(m-n)) %拟合了n个参数

rmse2=sqrt(sum((hg23(1)+x0*hg23(2:end)'-y0).^2)/(m-num)) %拟合了num个参数

(2)主成分分析法进行评价

对原始数据进行标准化处理

计算相关系数矩阵R

计算特征值和特征向量

选择p(p<=m)个主成分，计算综合评价值clc,clear

load gj.txt

gj=zscore(gj); %数据标准化

r=corrcoef(gj); %计算相关系数矩阵

%下面利用相关系数进行主成分分析，vec1的列为r的特征向量，即主成分系数

[vec1,lamda,rate]=pcacov(r) %lamda为r的特征值，即主成分系数

contr=cumsum(rate) %计算累积贡献率

f=repmat(sign(sum(vec1)),size(vec1,1),1);%¹构造与vec1同维数的元素为正负1的矩阵

vec2=vec1.*f %修改特征向量的正负号，是的每个特征向量的分量和为正

num=4; %num为选取的主成分的个数

df=gj*vec2(:,1:num); %计算各个主成分的得分

tf=df*rate(1:num)/100; %计算综合得分

[stf,ind]=sort(tf,'descend'); %把得分暗战从高到低的次序排列

stf=stf', ind=ind'

10.3 因子分析

10.3.1 因子载荷矩阵的估计方法

(1)主成分分析法

clc,clear

r=[1.000 0.577 0.509 0.387 0.462

0.577 1.000 0.599 0.389 0.322

0.509 0.599 1.000 0.436 0.426

0.387 0.389 0.436 1.000 0.523

0.462 0.322 0.426 0.523 1.000];

[vec1,val,rate]=pcacov(r);

f1=repmat(sign(sum(vec1)),size(vec1,1),1);

vec2=vec1.*f1;

f2=repmat(sqrt(val)',size(vec2,1),1);

a=vec2.*f2 %提出全部因子的载荷矩阵

a1=a(:,1) %提出一个因子的载荷矩阵

tcha1=diag(r-a1*a1') %计算一个因子的特殊方差

a2=a(:,[1,2])

tcha2=diag(r-a2*a2')

ccha2=r-a2*a2'-diag(tcha2) %求两个因子时的残差矩阵

con=cumsum(rate)

(2)主因子法

clc,clear

r=[1 1/5 -1/5;1/5 1 -2/5;-1/5 -2/5 1];

n=size(r,1); rt=abs(r);

rt(1:n+1:n^2)=0;

rstar=r;

rstar(1:n+1:n^2)=max(rt');

[vec1,val,rate]=pcacov(rstar)

f1=repmat(sign(sum(vec1)),size(vec1,1),1);

vec2=vec1.*f1

f2=repmat(sqrt(val)',size(vec2,1),1);

a=vec2.*f2

num=input('请选择公因子的个数：');

aa=a(:,1:num)

s1=sum(aa.^2)

s2=sum(aa.^2,2)

(3)最大似然估计法

Matlab工具箱中有factoran命令

clc,clear

r=[1 1/5 -1/5;1/5 1 -2/5;-1/5 -2/5 1];

[Lambda,Psi] = factoran(r,1,'xtype','cov') %Lambda返回的是因子载荷矩阵，Psi返回的是特殊方差

10.3.2 因子转换(正交变换)

方差最大法

四次方最大旋转

等量最大法clc,clear

r=[1 -1/3 2/3;-1/3 1 0;2/3 0 1];

[vec1,val,rate]=pcacov(r)

f1=repmat(sign(sum(vec1)),size(vec1,1),1);

vec2=vec1.*f1;

f2=repmat(sqrt(val)',size(vec2,1),1);

lambda=vec2.*f2

num=2;

[lambda2,t]=rotatefactors(lambda(:,1:num),'method', 'varimax') %对载荷矩阵进行旋转，其中lambda2位旋转载荷矩阵，t为变换的正交矩阵

10.3.3因子分析的步骤

选择分析的变量

计算所选原始变量相关的系数矩阵

提出公共因子

因子旋转

计算因子得分

10.3.4 因子分析案例

clc,clear

load ssgs.txt

n=size(ssgs,1);

x=ssgs(:,[1:4]); y=ssgs(:,5);

x=zscore(x);

r=corrcoef(x)

[vec1,val,con1]=pcacov(r)

f1=repmat(sign(sum(vec1)),size(vec1,1),1);

vec2=vec1.*f1;

f2=repmat(sqrt(val)',size(vec2,1),1);

a=vec2.*f2

num=input('请选择主因子的个数：');

am=a(:,[1:num]);

[bm,t]=rotatefactors(am,'method', 'varimax')

bt=[bm,a(:,[num+1:end])];

con2=sum(bt.^2)

check=[con1,con2'/sum(con2)*100]

rate=con2(1:num)/sum(con2)

coef=inv(r)*bm

score=x*coef

weight=rate/sum(rate)

Tscore=score*weight'

[STscore,ind]=sort(Tscore,'descend')

display=[score(ind,:)';STscore';ind']

[ccoef,p]=corrcoef([Tscore,y])

[d,dt,e,et,stats]=regress(Tscore,[ones(n,1),y]);

d,stats

10.4 判别分析

根据所研究的个体的观测指标来推断该个体所述类型。

10.4.1 距离判别

适用于连续性随机变量的判别类

10.4.2 Fisher判别

将表面上不易分类的数据透过投影得以分离

10.4.3 Bayes判别

与Bayes估计类似，利用后验概率分布

利用Matlab工具箱中的分类函数classify

clc,clear

p1=6/14;p2=8/14;

a=[24.8 24.1 26.6 23.5 25.5 27.4

-2.0 -2.4 -3.0 -1.9 -2.1 -3.1]';

b=[22.1 21.6 22.0 22.8 22.7 21.5 22.1 21.4

-0.7 -1.4 -0.8 -1.6 -1.5 -1.0 -1.2 -1.3]';

n1=6;n2=8;

train=[a;b]; %train为已知样本

group=[ones(n1,1);2*ones(n2,1)]; %已知样本类别标识

prior=[p1; p2]; %已知样本的闲言概率

sample=train; %sample一般为未知样本，这里是准备回代检验误判

[x1,y1]=classify(sample,train,group,'linear',prior) %线性分类

[x2,y2]=classify(sample,train,group,'quadratic',prior) %二次分类

%函数的第二个返回值为误判率

10.4.4 应用举例

clc,clear

a=[978898743621682

867599546342414

768537646352535];

train=a(:,[1:12])';

sample=a(:,[13:end])';

group=[ones(7,1);2*ones(5,1)]; %已知样本分类

[x1,y1]=classify(sample,train,group,'mahalanobis') %马氏距离分类

[x2,y2]=classify(sample,train,group,'linear') %线性分类

[x3,y3]=classify(sample,train,group,'quadratic') %二次分类

10.5 典型相关分析

10.5.1 原始变量和典型变量之间的相关性

(1)原始变量和典型变量之间的相关系数

(2)各组原始变量被典型变量所解释的反差

10.5.2 典型相关系数的检验

(1)计算样本的协方差矩阵

(2)整体检验

(3)部分总体典型相关系数为零的检验

10.5.3 典型相关分析案例

(1)职业满意度典型相关分析

clc,clear

load r.txt

n1=5;n2=7;num=min(n1,n2);

s1=r([1:n1],[1:n1]);

s12=r([1:n1],[n1+1:end]);

s21=s12';

s2=r([n1+1:end],[n1+1:end]);

m1=inv(s1)*s12*inv(s2)*s21;

m2=inv(s2)*s21*inv(s1)*s12;

[vec1,val1]=eig(m1);

for i=1:n1

vec1(:,i)=vec1(:,i)/sqrt(vec1(:,i)'*s1*vec1(:,i));

vec1(:,i)=vec1(:,i)/sign(sum(vec1(:,i)));

end

val1=sqrt(diag(val1));

[val1,ind1]=sort(val1,'descend');

a=vec1(:,ind1(1:num))

dcoef1=val1(1:num) %提出典型相关系数

flag=1;

xlswrite('bk.xls',a,'Sheet1','A1')

flag=n1+2; str=char(['A',int2str(flag)]);

xlswrite('bk.xls',dcoef1','Sheet1',str)

[vec2,val2]=eig(m2);

for i=1:n2

vec2(:,i)=vec2(:,i)/sqrt(vec2(:,i)'*s2*vec2(:,i));

vec2(:,i)=vec2(:,i)/sign(sum(vec2(:,i)));

end

val2=sqrt(diag(val2));

[val2,ind2]=sort(val2,'descend');

b=vec2(:,ind2(1:num))

dcoef2=val2(1:num)

flag=flag+2; str=char(['A',int2str(flag)]);

xlswrite('bk.xls',b,'Sheet1',str)

flag=flag+n2+1; str=char(['A',int2str(flag)]);

xlswrite('bk.xls',dcoef2','Sheet1',str)

x_u_r=s1*a

y_v_r=s2*b

x_v_r=s12*b

y_u_r=s21*a

flag=flag+2; str=char(['A',int2str(flag)]);

xlswrite('bk.xls',x_u_r,'Sheet1',str)

flag=flag+n1+1; str=char(['A',int2str(flag)]);

xlswrite('bk.xls',y_v_r,'Sheet1',str)

flag=flag+n2+1; str=char(['A',int2str(flag)]);

xlswrite('bk.xls',x_v_r,'Sheet1',str)

flag=flag+n1+1; str=char(['A',int2str(flag)]);

xlswrite('bk.xls',y_u_r,'Sheet1',str)

mu=sum(x_u_r.^2)/n1

mv=sum(x_v_r.^2)/n1

nu=sum(y_u_r.^2)/n2

nv=sum(y_v_r.^2)/n2

fprintf('X组的原始变量被u1~u%d解释比例为%f\n',num,sum(mu));

fprintf('Y组的原始变量被v1~v%d解释比例为%f\n',num,sum(nv));

(2)中国城市竞争力和基础设施的典型相关分析

(3)城市竞争力和基础设施关系的典型相关分析

(4)城市竞争力和基础设施关系的经济分析

10.6 对应分析

对原始数据采用适当的标度方法，把R型和Q型分析结合起来，同时得到两方面的结果

clc, clear

format long g

a=load('dy.txt');

T=sum(sum(a));

P=a/T;

r=sum(P,2), c=sum(P) %计算边缘分布

Row_prifile=a./repmat(sum(a,2),1,size(a,2)) %计算行轮廓分布阵

B=(P-r*c)./sqrt((r*c)); %计算标准化数据阵B

[u,s,v]= svd(B,'econ') %对标准化后的数据阵B作奇异值分解

w=sign(repmat(sum(v),size(v,1),1)) %修改特征向量的符号矩阵

%使得v中的每一个列向量的分量和大于0

ub=u.*w %修改特征向量的正负号

vb=v.*w

lamda=diag(s).^2 %计算惯量

ksi2square=T*(lamda) %计算卡方统计量的分解

T_ksi2square=sum(ksi2square) %计算总卡方统计量

con_rate=lamda/sum(lamda) %计算贡献率

cum_rate=cumsum(con_rate) %计算累积贡献率

beta=diag(r.^(-1/2))*ub; %求加权特征向量

G=beta*s %求行轮廓坐标

alpha=diag(c.^(-1/2))*vb;

F=alpha*s

num=size(G,1);

rang=minmax(G(:,1)'); %坐标的取值范围

delta=(rang(2)-rang(1))/(8*num); %画图的标注位置调整量

ch='LPSBEM';

hold on

for i=1:num

plot(G(i,1),G(i,2),'*','Color','k','LineWidth',1.3) %画行点散布图

text(G(i,1)+delta,G(i,2),ch(i)) %对行点进行标注

plot(F(i,1),F(i,2),'H','Color','k','LineWidth',1.3)

text(F(i,1)+delta,F(i,2),int2str(i+1972))

end

xlabel('dim1'), ylabel('dim2')

xlswrite('tt1',[diag(s),lamda,ksi2square,con_rate,cum_rate])

format

10.7 多维标度法

距离阵

欧几里得阵clc, clear

D=[0, 1, sqrt(3), 2, sqrt(3), 1, 1; zeros(1,2),1, sqrt(3), 2, sqrt(3), 1

zeros(1,3),1, sqrt(3), 2, 1;zeros(1,4), 1, sqrt(3), 1

zeros(1,5), 1, 1; zeros(1,6), 1; zeros(1,7)] %原始距离矩阵的上三角元素

d=D+D'; %构造完整的距离矩阵

%d=nonzeros(D')'; %转换成pdist函数输出格式的数据

[y,eigvals]=cmdscale(d) %求经典解，d可以为实对称矩阵或pdist函数的行向量输出

plot(y(:,1),y(:,2),'o','Color','k','LineWidth',1.3) %画出点的坐标

%下面通过求特征值求经典解

D2=D+D'; %构造对称距离

A=-D2.^2/2; %构造A矩阵

n=size(A,1);

H=eye(n)-ones(n)/n; %构造H矩阵

B=H*A*H %构造B矩阵

[vec1,val1]=eig(B); %求B矩阵的特征向量vec1和特征值val1

[val2,ind]=sort(diag(val1),'descend') %把特征按从大到小排列

vec2=vec1(:,ind) %相应的把特征向量也重新排序

vec3=orth(vec2(:,[1,2])); %构造正交特征向量

point=[vec3(:,1)*sqrt(val2(1)),vec3(:,2)*sqrt(val2(2))] %求点的坐标

hold on

plot(point(:,1),point(:,2),'D','Color','k','LineWidth',1.3) %验证得到的解和Matlab不一致

theta=0.409; %旋转的角度

T=[cos(theta),-sin(theta);sin(theta),cos(theta)];

Tpoint=point*T; %把特征向量进行一个正交变换

plot(Tpoint(:,1),Tpoint(:,2),'+','Color','k','LineWidth',1.3) %验证这样得到的解和Matlab一致

legend('Matlab命令cmdscale求得的解','按照算法求得的一个解','正交变换后得到的雨cmdscale相同的解','Location','best')