计量地理学

for z=1:length(y)
newy(z)=y(length(y)+1-z);
end
fprintf('%g\n',newy)
rate=y/sum(y);
fprintf('\n 贡献率：\n')
newrate=newy/sum(newy)
sumrate=0;
newi=[];
for k=length(y):-1:1
主成分 zi 的贡献率为
累计贡献率为
λi
p
∑ λk
k =1
(i = 1,2,L, p)
i
∑ λk
k =1
p
∑ λk
k =1
(i = 1,2,L, p)
一般取累计贡献率达 85%~95%的特征值 λ1, λ2 ,L, λm 所对应的第 1、第 2，…， 第 m（m≤p）个主成分。
④ 计算主成分载荷。 其计算公式为
福 州 583.13 0.273 3 2 152 288 6 555 351 8 851 7 190 467 524 5 030 220 69.59 680 607
厦 门 128.99 0.486 5 333 374 5 751 124 3 728 2 570 418 758 2 108 331 46.93 657 484
z 22
L
z
2m
⎥ ⎥
⎢M M
M⎥
⎢ ⎣
z
n1
zn2
L
z
nm
⎥ ⎦
2.程序结构及函数作用
（4）
在软件 Matlab 中实现主成分分析可以采取两种方式实现：一是通过编程来 实现；二是直接调用 Matlab 中自带程序实现。下面主要介绍利用 Matlab 的矩阵 计算功能编程实现主成分分析。
程序结构
子函数
城乡居 民年底 储蓄余额 /万元
在岗职 工人数 /万人
在岗职工 工资总额
/万元
北 京 1 249.90 0.597 8 1 843 427 19 999 706 20 323 45 562 2 790 863 26 806 646 410.80 5 773 301
天 津 910.17 0.580 9 1 501 136 22 645 502 3 259 26 317 1 128 073 11 301 931 202.68 2 254 343
fprintf('特征向量(vec)及特征值(val)：\n') [vec,val]=eig(std) %求特征值(val)及特征向量(vec)
newval=diag(val) ;
[y,i]=sort(newval) ;
%对特征根进行排序，y 为排序结果，i 为索引
fprintf('特征根排序：\n')
石 家 庄 875.40 0.233 2 2 918 680 6 885 768 2 929 1 911 352 348 7 095 875 95.60 758 877
太 原 299.92 0.656 3 236 038 2 737 750 1 937 11 895 203 277 3 943 100 88.65 654 023
南 昌 424.20 0.398 8 688 289 2 305 881 3 674 3 189 167 714 2 640 460 62.08 479 ,555
济 南 557.63 0.408 5 1 486 302 6 285 882 5 915 11 775 460 690 4 126 970 83.31 756 696
for i=1:a
for j=1:b
std(i,j)= vector(i,j)/cwsum(j);
end
end
cwfac.m
%cwfac.m
function result=cwfac(vector); fprintf('相关系数矩阵:\n')
std=CORRCOEF(vector) %计算相关系数矩阵
大 连 545.31 0.494 6 1 879 739 8 426 385 10 780 19 187 709 227 7 556 796 94.15 965 922
长 春 691.23 0.406 8 1 853 210 5 966 343 4 810 9 532 357 096 4 803 744 102.63 884 447
sumrate=sumrate+rate(k);
newi(length(y)+1-k)=i(k);
if sumrate>0.85 break;
end end
%记下累积贡献率大于 85%的特征值的序号放入 newi 中
4
《计量地理学》（徐建华，高等教育出版社，2006）配套实习指导
fprintf('主成分数：%g\n\n',length(newi)); fprintf('主成分载荷：\n')
广 州 685.00 0.621 4 1 600 738 23 348 139 22 007 23 854 1 761 499 20 401 811 182.81 3 047 594
深 圳 119.85 0.793 1 299 662 20 368 295 8 754 4 274 1 847 908 9 519 900 91.26 1 890 338
lij = p(zi , x j ) = λi eij (i, j = 1,2,L, p)
（3）
得到各主成分的载荷以后，还可以按照教材中（3.5.2）式进一步计算，得到 各主成分的得分
2
《计量地理学》（徐建华，高等教育出版社，2006）配套实习指导
⎡ z11 z12 L z1m ⎤
Z
=
⎢ ⎢
z
Leabharlann Baidu
21
for p=1:length(newi)
for q=1:length(y)
result(q,p)=sqrt(newval(newi(p)))*vec(q,newi(p));
end
end
%计算载荷
disp(result)
cwscore.m
%cwscore.m,计算得分 function score=cwscore(vector1,vector2); sco=vector1*vector2; csum=sum(sco,2); [newcsum,i]=sort(-1*csum); [newi,j]=sort(i); fprintf('计算得分：\n') score=[sco,csum,j]
λi (i = 1,2,L, p) ，并使其按大小顺序排列，即 λ1 ≥ λ2 ≥ L ≥ λp ≥ 0 ；然后分别求出
p
∑ 对应于特征值 λi 的特征向量 ei (i = 1,2,L, p) 。这里要求 ei =1，即 ei2j = 1，其中 j =1
eij 表示向量 ei 的第 j 个分量。
③ 计算主成分贡献率及累计贡献率。
k =1
n
n
∑ ∑ (xki − xi )2 (xkj − x j )2
k =1
k =1
（2）
因为R是实对称矩阵（即rij=rji），所以只需计算上三角元素或下三角元素即可。 ② 计算特征值与特征向量
1
《计量地理学》（徐建华，高等教育出版社，2006）配套实习指导
首 先 解 特 征 方 程 λI − R = 0 ， 通 常 用 雅 可 比 法 （ Jacobi ） 求 出 特 征 值
武 汉 740.20 0.586 9 1 211 291 7 506 085 9 793 15 442 604 658 5 748 055 149.20 1 314 766
长 沙 582.47 0.310 7 1 146 367 3 098 179 8 706 5 718 323 660 3 461 244 69.57 596 986
%得分矩阵：sco 为各主成分得分；csum 为综合得分；j 为排序结果
cwprint.m
%cwprint.m function print=cwprint(filename,a,b);
%filename 为文本文件名，a 为矩阵行数(样本数)，b 为矩阵列数(变量指标数) fid=fopen(filename,'r') vector=fscanf(fid,'%g',[a b]); fprintf('标准化结果如下：\n') v1=cwstd(vector) result=cwfac(v1); cwscore(v1,result);
4.程序测试
原始数据
中国大陆 35 个大城市某年的 10 项社会经济统计指标数据见下表。
5
《计量地理学》（徐建华，高等教育出版社，2006）配套实习指导
城市 名称
年底总 人口 /万人
非农业 人口比 例/%
农业 总产值 /万元
工业 总产值 /万元
客运 总量 /万人
货运 总量 /104 t
地方财政 预算内收 入/万元
呼和浩特 207.78 0.441 2 365 343 816 452 2 351 2 623 105 783 1 396 588 42.11 309 337
沈 阳 677.08 0.629 9 1 295 418 5 826 733 7 782 15 412 567 919 9 016 998 135.45 1 152 811
Cwprint.m
主函数
Cwstd.m
Cwfac.m
Cwscore.m
函数作用
Cwstd.m——用总和标准化法标准化矩阵。 Cwfac.m——计算相关系数矩阵；计算特征值和特征向量；对主成分进行排 序；计算各特征值贡献率；挑选主成分（累计贡献率大于 85%），输出主成分个 数；计算主成分载荷。 Cwscore.m——计算各主成分得分、综合得分并排序。 Cwprint.m——读入数据文件；调用以上 3 个函数并输出结果。
哈 尔 滨 927.09 0.462 7 2 663 855 4 186 123 6 720 7 520 481 443 6 450 020 172.79 1 309 151
上 海 1 313.12 0.738 4 2 069 019 54 529 098 6 406 44 485 4 318 500 25 971 200 336.84 5 605 445
南 京 537.44 0.534 1 989 199 13 072 737 14 269 11 193 664 299 5 680 472 113.81 1 357 861
杭 州 616.05 0.355 6 1 414 737 12 000 796 17 883 11 684 449 593 7 425 967 96.90 1 180 947
3
《计量地理学》（徐建华，高等教育出版社，2006）配套实习指导
3.源程序
cwstd.m
%cwstd.m,用总和标准化法标准化矩阵
function std=cwstd(vector) cwsum=sum(vector,1); [a,b]=size(vector);
%对列求和 %矩阵大小,a 为行数,b 为列数
青 岛 702.97 0.369 3 2 382 320 11 492 036 13 408 17 038 658 435 4 978 045 103.52 961 704
郑 州 615.36 0.342 4 677 425 5 287 601 10 433 6 768 387 252 5 135 338 84.66 696 848
《计量地理学》（徐建华，高等教育出版社，2006）配套实习指导
§10 利用 Matlab 编程实现主成分分析
1.概述
Matlab 语言是当今国际上科学界(尤其是自动控制领域)最具影响力、也是最 有活力的软件。它起源于矩阵运算，并已经发展成一种高度集成的计算机语言。 它提供了强大的科学运算、灵活的程序设计流程、高质量的图形可视化与界面设 计、与其他程序和语言的便捷接口的功能。Matlab 语言在各国高校与研究单位起 着重大的作用。主成分分析是把原来多个变量划为少数几个综合指标的一种统计 分析方法，从数学角度来看，这是一种降维处理技术。
宁 波 538.41 0.254 7 1 428 235 10 622 866 22 215 10 298 501 723 5 246 350 62.15 824 034
合 肥 429.95 0.318 4 628 764 2 514 125 4 893 1 517 233 628 1 622 931 47.27 369 577
主成分分析计算步骤：
① 计算相关系数矩阵
⎡ r11 r12 L r1p ⎤
R
=
⎢ ⎢
r21
⎢M
r22 M
L
r2
p
⎥ ⎥
M⎥
⎢⎢⎣rp1
rp2
L
rpp
⎥ ⎥⎦
（1）
在（1）式中，rij（i，j=1，2，…，p）为原变量的xi与xj之间的相关系数， 其计算公式为
rij =
n
∑ (xki − xi )(xkj − x j )