做网站用的大图,免费的个人简历模板下载,wordpress 媒体库,濮阳网站设计公司【数学建模】《实战数学建模#xff1a;例题与讲解》第十三讲-相关分析#xff08;含Matlab代码#xff09; 基本概念典型相关分析综合评价模型对应分析因子分析聚类分析 习题10.41. 题目要求2.解题过程3.程序 习题10.51. 题目要求2.解题过程3.程序 习题10.6#xff08;1例题与讲解》第十三讲-相关分析含Matlab代码 基本概念典型相关分析综合评价模型对应分析因子分析聚类分析 习题10.41. 题目要求2.解题过程3.程序 习题10.51. 题目要求2.解题过程3.程序 习题10.611. 题目要求2.解题过程——对应分析3.程序4.结果 习题10.621. 题目要求2.解题过程——R型因子分析3.程序4.结果 习题10.631. 题目要求2.解题过程——聚类分析3.程序4.结果 本系列侧重于例题实战与讲解希望能够在例题中理解相应技巧。文章开头相关基础知识只是进行简单回顾读者可以搭配课本或其他博客了解相应章节然后进入本文正文例题实战效果更佳。
如果这篇文章对你有帮助欢迎点赞与收藏~
基本概念
变量类型在进行相关分析之前首先要确定所研究的变量类型。这些变量可以是连续的如身高、体重或者离散的如性别、婚姻状况。相关系数相关分析的核心是计算相关系数这是一个度量值表明两个变量之间的关系有多紧密。最常用的相关系数是皮尔逊相关系数Pearson correlation coefficient用于度量两个连续变量之间的线性关系。线性与非线性关系皮尔逊相关系数主要用于评估线性关系。对于非线性关系可能需要使用其他类型的相关系数如斯皮尔曼等级相关系数Spearman’s rank correlation coefficient。方向相关系数的值范围通常在-1到1之间。一个正的相关系数意味着一个变量增加时另一个变量也增加负的相关系数则意味着一个变量增加时另一个变量减少。统计显著性仅仅计算出相关系数是不够的还需要评估这种相关性是否具有统计显著性。通常这通过进行假设检验如t检验来实现。因果关系值得注意的是即使两个变量之间存在强相关性也不能自动推断出因果关系。相关性只能揭示变量之间的关联而不是因果。
本文中的相关知识
典型相关分析
计算两组变量表示为 X 组和 Y 组的典型变量这些典型变量是通过线性组合构成目的是最大化两组变量间的相关性。 计算原始变量与典型变量之间的相关系数这展示了原始变量与新构造的典型变量之间的关联强度。 分析典型变量之间的典型相关系数这反映了两组变量间的整体关联强度。
综合评价模型
用于评估多个指标的综合影响。在这个例子中使用了无量纲化方法来处理不同单位或量级的数据使得比较和评估成为可能。 使用欧氏距离和绝对值距离来评估不同湖泊水质的富营养化等级。
对应分析
这是一种用于探索分类数据的多元统计技术它可以揭示分类变量之间的关系。在这个例子中通过对应分析可以看出不同地区和指标之间的关系。
因子分析
用于数据简化和变量归约。在这个例子中通过因子分析将多个变量简化为少数几个因子从而使得解释变得更加简洁。
聚类分析
包括 R 型和 Q 型聚类分别用于变量和样本点的分类。通过聚类分析可以将变量或样本点分为不同的组从而发现数据的内在结构。
习题10.4
1. 题目要求 2.解题过程
解 根据计算我们得到了 X \mathbf{X} X 组的典型变量为 u 1 0.7689 x 1 0.2721 x 2 u_10.7689x_10.2721x_2 u10.7689x10.2721x2 u 2 − 1.4787 x 1 1.6443 x 2 u_2-1.4787x_11.6443x_2 u2−1.4787x11.6443x2 我们也得到了原始变量与 X \mathbf{X} X 组典型变量之间的相关系数如表 1、原始变量与 Y \mathbf{Y} Y 组典型变量之间的相关系数如表 2以及两组典型变量之间的典型相关系数如表 3。 表 1原始变量与 X 组典型变量之间的相关系数 x 1 x_{1} x1 x 2 x_{2} x2 y 1 y_{1} y1 y 2 y_{2} y2 y 3 y_{3} y3 u 1 u_{1} u10.9865860.8872150.2897050.6757040.35394 u 2 u_{2} u2-0.163240.4613560.158162-0.020580.05631表 2原始变量与 Y 组典型变量之间的相关系数 x 1 x_{1} x1 x 2 x_{2} x2 y 1 y_{1} y1 y 2 y_{2} y2 y 3 y_{3} y3 v 1 v_{1} v10.678720.6103580.4211150.9822030.514487 v 2 v_{2} v2-0.03050.0862120.846397-0.110140.301341表 3两组典型变量之间的典型相关系数 10.687920.1869从这些表中我们可以看到两个表示外出活动特性的变量 x 1 x_1 x1 和 x 2 x_2 x2 与 u 1 u_1 u1 的相关系数较大说明 u 1 u_1 u1 可以被视为形容外出活动特性的指标。在第一对典型变量中 v 1 v_1 v1 与 y 2 y_2 y2 代表家庭年收入的相关系数较大因此我们可以认为 v 1 v_1 v1 主要代表了家庭的年收入。此外 u 1 u_1 u1 和 v 1 v_1 v1 之间的典型相关系数为 0.6879说明这两个典型变量在一定程度上是相关的。 最后对于原始变量的解释度 u 1 u_1 u1 和 v 1 v_1 v1 解释的本组原始变量的比率分别为 0.8803 和 0.4689。这表示 u 1 u_1 u1 和 v 1 v_1 v1 在一定程度上能够解释原始变量的变异。而且 X \mathbf{X} X 组的原始变量被 u 1 u_1 u1 和 u 2 u_2 u2 完全100%解释了而 Y \mathbf{Y} Y 组的原始变量被 v 1 v_1 v1 和 v 2 v_2 v2 解释了 74.2%说明典型变量确实能够在很大程度上表现出原始变量的特性。
3.程序
求解的MATLAB程序如下
clc, clear;% 加载数据
correlation_matrix load(data104.txt);% X组和Y组变量的数量
num_X 2;
num_Y 3;
num_min min(num_X, num_Y);% 提取X与XY与YY与X的相关系数
correlation_XX correlation_matrix([1:num_X], [1:num_X]);
correlation_YX correlation_matrix([1:num_X], [num_X1:end]);
correlation_XY correlation_YX;
correlation_YY correlation_matrix([num_X1:end], [num_X1:end]);% 计算矩阵M1和M2
matrix_M1 inv(correlation_XX) * correlation_YX * inv(correlation_YY) * correlation_XY;
matrix_M2 inv(correlation_YY) * correlation_XY * inv(correlation_XX) * correlation_YX;% 求M1的特征向量和特征值
[eigVec_M1, eigVal_M1] eig(matrix_M1);
for i 1:num_X% 特征向量归一化满足as1a1特征向量乘±1保证所有分量和为正eigVec_M1(:, i) eigVec_M1(:, i) / sqrt(eigVec_M1(:, i) * correlation_XX * eigVec_M1(:, i));eigVec_M1(:, i) eigVec_M1(:, i) / sign(sum(eigVec_M1(:, i)));
end
% 计算特征值的平方根按照从大到小排列
eigVal_M1 sqrt(diag(eigVal_M1));
[eigVal_M1, ind1] sort(eigVal_M1, descend);% 取出X组的系数阵和典型相关系数
coef_X eigVec_M1(:, ind1(1:num_min));
canCorrelation_X eigVal_M1(1:num_min);% 求M2的特征向量和特征值
[eigVec_M2, eigVal_M2] eig(matrix_M2);
for i 1:num_Y% 特征向量归一化满足bs2b1特征向量乘±1保证所有分量和为正eigVec_M2(:, i) eigVec_M2(:, i) / sqrt(eigVec_M2(:, i) * correlation_YY * eigVec_M2(:, i));eigVec_M2(:, i) eigVec_M2(:, i) / sign(sum(eigVec_M2(:, i)));
end
% 计算特征值的平方根按照从大到小排列
eigVal_M2 sqrt(diag(eigVal_M2));
[eigVal_M2, ind2] sort(eigVal_M2, descend);% 取出Y组的系数阵和典型相关系数
coef_Y eigVec_M2(:, ind2(1:num_min));
canCorrelation_Y eigVal_M2(1:num_min);% 计算XuYvXvYu的相关系数
correlation_XU correlation_XX * coef_X;
correlation_YV correlation_YY * coef_Y;
correlation_XV correlation_YX * coef_Y;
correlation_YU correlation_XY * coef_X;% 计算X组、Y组原始变量被ui、vi解释的方差比例
ratio_XU sum(correlation_XU.^2) / num_X;
ratio_XV sum(correlation_XV.^2) / num_X;
ratio_YU sum(correlation_YU.^2) / num_Y;
ratio_YV sum(correlation_YV.^2) / num_Y;% 输出结果
fprintf(X组的原始变量被u1~u%d解释的比例为%f\n, num_min, sum(ratio_XU));
fprintf(Y组的原始变量被v1~v%d解释的比例为%f\n, num_min, sum(ratio_YV));习题10.5
1. 题目要求 2.解题过程
解
在进行综合评价之前首先要对评价的指标进行分析。通常的评价指标分为效益型、成本型和固定型指标。效益型指标是指那些数值越大影响力越大的统计指标(也称为正向型指标);成本型指标是指数值越小越好的指标(亦称为逆向型指标);而固定型指标是指数值越接近某个常数越好的指标(又称为适度型指标)。如果各评价指标的属性不一致,则在进行综合评估时容易发生偏差,必须先对各评价指标统一属性。 建立无量纲化实测数据矩阵和评价标准矩阵。实测数据矩阵为 A \mathbf{A} A 和等级标准矩阵为 B \mathbf{B} B然后建立无量纲化实测数据矩阵 C \mathbf{C} C 和无量纲化等级标准矩阵 D \mathbf{D} D。其中 c i j c_{ij} cij 和 d k t d_{kt} dkt 的计算方法如下 c i j c_{ij} cij 和 d k t d_{kt} dkt 的计算方式 c i j { a i j / max i a i j , if j ≠ 3 min i a i j / a i j , if j 3 c_{ij} \begin{cases} a_{ij} / \max_i a_{ij}, \text{if } j \neq 3 \\ \min_i a_{ij} / a_{ij}, \text{if } j 3 \end{cases} cij{aij/maxiaij,miniaij/aij,if j3if j3 d k t { b k t / max i b k t , if k ≠ 3 min i b k t / b k t , if k 3 d_{kt} \begin{cases} b_{kt} / \max_i b_{kt}, \text{if } k \neq 3 \\ \min_i b_{kt} / b_{kt}, \text{if } k 3 \end{cases} dkt{bkt/maxibkt,minibkt/bkt,if k3if k3 计算 D \mathbf{D} D 的各行向量的均值 μ i \mu_i μi、标准差 s i s_i si 和变异系数 w i w_i wi以及权向量 w \boldsymbol{w} w。计算公式如下 μ i \mu_i μi s i s_i si 和 w i w_i wi 的计算方式 μ i 0.2 ∑ j 1 5 d i j s i ∑ j 1 5 ( d i j − μ i ) 2 4 w i s i / μ i \mu_i 0.2 \sum_{j1}^{5} d_{ij} s_i \sqrt{\frac{\sum_{j1}^{5} (d_{ij}-\mu_i)^2}{4}} w_i s_i / \mu_i μi0.2j1∑5dijsi4∑j15(dij−μi)2 wisi/μi w \boldsymbol{w} w 的计算方式 w [ 0.277 , 0.2447 , 0.2347 , 0.2442 ] \boldsymbol{w} [0.277,0.2447,0.2347,0.2442] w[0.277,0.2447,0.2347,0.2442] 建立综合评价模型计算 C \mathbf{C} C 中各个行向量到 D \mathbf{D} D 中各个列向量的欧氏距离 x i j x_{ij} xij 和绝对值距离 y i j y_{ij} yij。若 x i k min 1 ≤ j ≤ 5 x i j x_{ik}\min_{1\leq j\leq5}{x_{ij}} xikmin1≤j≤5xij 或 y i k min 1 ≤ j ≤ 5 y i j y_{ik}\min_{1\leq j\leq5}{y_{ij}} yikmin1≤j≤5yij 则表明第 i i i 个湖泊属于第 k k k 级。 然后我们可以得到以下欧氏距离判别表和绝对值距离判别表。 x i j x_{ij} xij 和 y i j y_{ij} yij 的计算方式 x i j ∑ k 1 4 ( c i k − d k j ) 2 y i j ∑ k 1 4 ∣ c i k − d k j ∣ x_{ij} \sqrt{\sum_{k1}^{4} (c_{ik}-d_{kj})^2}\\ y_{ij} \sum_{k1}^{4} |c_{ik}-d_{kj}| xijk1∑4(cik−dkj)2 yijk1∑4∣cik−dkj∣
然后我们可以得到以下欧氏距离判别表和绝对值距离判别表。
欧氏距离判别表
湖泊 x i 1 x_{i1} xi1 x i 2 x_{i2} xi2 x i 3 x_{i3} xi3 x i 4 x_{i4} xi4 x i 5 x_{i5} xi5级别杭州西湖1.84721.83121.73741.37690.28815武汉东湖1.59591.57981.48591.12710.50345青海湖0.21850.20450.13670.33831.79173巢湖1.32011.30381.20820.83920.95914滇池1.07931.06500.98670.73281.34504
绝对值距离判别表
湖泊 y i 1 y_{i1} yi1 y i 2 y_{i2} yi2 y i 3 y_{i3} yi3 y i 4 y_{i4} yi4 y i 5 y_{i5} yi5级别杭州西湖3.66313.63033.43742.67830.32315武汉东湖3.14363.11082.91782.15870.84275青海湖0.40620.37340.21100.57873.58003巢湖2.40712.37432.18141.42231.57914滇池1.67011.63741.44441.06602.31614
从上面的计算可知尽管欧几里得距离与绝对值距离意义不同但是对各湖泊水质的富营养化的评价等级是一样的表明此处给出的方法具有稳定性。
3.程序
求解的MATLAB程序如下
% 清除环境变量
clc; clear;% 初始化实测数据矩阵
rawData [130,10.3,0.35,2.76;105,10.7,0.4,2.0;20,1.4,4.5,0.22;30,6.26,0.25,1.67;20,10.13,0.5,0.23];% 初始化等级标准矩阵
standardMatrix [1,4,23,110,660;0.09,0.36,1.8,7.1,27.1;37,12.,2.4,0.55,0.17;0.02,0.06,0.31,1.2,4.6];% 进行无量纲化处理
normalizedRawData rawData ./ repmat(max(rawData), size(rawData, 1), 1);
normalizedRawData(:, 3) min(rawData(:, 3)) ./ rawData(:, 3);normalizedStandardMatrix standardMatrix ./ repmat(max(standardMatrix, [], 2), 1, size(standardMatrix, 2));
normalizedStandardMatrix(3, :) min(standardMatrix(3, :)) ./ standardMatrix(3, :);% 计算均值和标准差
average mean(standardMatrix, 2);
standardDeviation std(standardMatrix, [], 2);% 计算权重
weights standardDeviation ./ average;
weights weights / sum(weights);% 计算欧式距离
euclideanDistances dist(normalizedRawData, normalizedStandardMatrix);% 计算绝对值距离
absoluteDistances mandist(normalizedRawData, normalizedStandardMatrix);习题10.61
1. 题目要求 2.解题过程——对应分析
解
分别用 i 1 , … , 16 i1,\dots,16 i1,…,16 表示以 a i j a_{ij} aij 表示第 i 个地区第 j 个指标变量 x j x_j xj 的取值。
记 A ( a i j ) 16 × 6 \mathbf{A}(a_{ij})_{16\times6} A(aij)16×6 记 a i ⋅ ∑ j 1 6 a i j , a ⋅ j ∑ i 1 16 a i j a_{i\cdot}\sum_{j1}^{6}a_{ij},a_{\cdot j}\sum_{i1}^{16}a_{ij} ai⋅j1∑6aij,a⋅ji1∑16aij 首先把 A \mathbf{A} A 化为规格化的“概率”矩阵 P \mathbf{P} P 记 P ( p i j ) 16 × 6 \mathbf{P}(p_{ij})_{16\times6} P(pij)16×6 其中 p i j a i j / T p_{ij}a_{ij}/\mathbf{T} pijaij/T T ∑ i 1 16 ∑ j 1 6 a i j \mathbf{T}\sum_{i1}^{16}\sum_{j1}^{6}a_{ij} T∑i116∑j16aij。再对数据进行对应变换令 B ( b i j ) 16 × 6 \mathbf{B}(b_{ij})_{16\times6} B(bij)16×6 其中 b i j p i j − p i ⋅ p ⋅ j p i ⋅ p ⋅ j a i j − a i ⋅ a ⋅ j / T a i ⋅ a ⋅ j , i 1 , 2 , … , 16 , j 1 , 2 , … , 6. b_{ij}\frac{p_{ij}-p_{i\cdot}p_{\cdot j}}{\sqrt{p_{i\cdot}p_{\cdot j}}} \frac{a_{ij}-a_{i\cdot}a_{\cdot j}/\mathbf{T}}{\sqrt{a_{i\cdot}a_{\cdot j}}}, i1,2,\dots, 16,j1,2,\dots,6. bijpi⋅p⋅j pij−pi⋅p⋅jai⋅a⋅j aij−ai⋅a⋅j/T,i1,2,…,16,j1,2,…,6. 对 B \mathbf{B} B 进行奇异值分解 B U Λ V T \mathbf{B}\mathbf{U}\varLambda \mathbf{V}^{\mathbf{T}} BUΛVT其中 U \mathbf{U} U 为 16 × 16 16\times 16 16×16 正交矩阵 V \mathbf{V} V为 6 × 6 6\times6 6×6 正交矩阵 Λ [ Λ m 0 0 0 ] \varLambda \begin{bmatrix} \varLambda_m 0\\ 00 \end{bmatrix} Λ[Λm000]这里 $ \varLambda_mdiag(d_1,\dots,d_m) $ 其中 d i ( i 1 , 2 , … , m ) d_i(i1,2,\dots,m) di(i1,2,…,m) 为 B \mathbf{B} B 的奇异值。
记 U [ U 1 ⋮ U 2 ] , V [ V 1 ⋮ V 2 ] \mathbf{U}[\mathbf{U_1} \vdots \mathbf{U_2}],\mathbf{V}[\mathbf{V_1} \vdots \mathbf{V_2}] U[U1⋮U2],V[V1⋮V2] 其中 U 1 \mathbf{U_1} U1 为 16 × m 16\times m 16×m 的列正交矩阵 V 1 \mathbf{V_1} V1 为 6 × m 6\times m 6×m 的列正交矩阵则 B \mathbf{B} B 的奇异值分解式等价于 B U 1 Λ V 1 T \mathbf{B}\mathbf{U_1}\varLambda\mathbf{V_1}^T BU1ΛV1T。
记 D r d i a g ( p 1 ⋅ , p 2 ⋅ , … , p 16 ⋅ ) , D c d i a g ( p ⋅ 1 , p ⋅ 2 , … , p ⋅ 6 ) \mathbf{D_r}diag(p_{1\cdot},p_{2\cdot},\dots,p_{16\cdot}),\mathbf{D_c}diag(p_{\cdot1},p_{\cdot2},\dots,p_{\cdot6}) Drdiag(p1⋅,p2⋅,…,p16⋅),Dcdiag(p⋅1,p⋅2,…,p⋅6) 其中 p i ⋅ ∑ j 1 6 p i j p_{i\cdot}\sum_{j1}^{6}p_{ij} pi⋅∑j16pij p ⋅ j ∑ i 1 16 p i j p_{\cdot j}\sum_{i1}^{16}p_{ij} p⋅j∑i116pij。则列轮廓的坐标为 F D c − 1 / 2 V 1 Λ m \mathbf{F}\mathbf{D}_{c}^{-1/2}\mathbf{V_1}\varLambda_m FDc−1/2V1Λm 行轮廓的坐标为 G D r − 1 / 2 V 1 Λ m \mathbf{G}\mathbf{D}_{r}^{-1/2}\mathbf{V_1}\varLambda_m GDr−1/2V1Λm 。最后通过贡献率的比较确定需要截取的维数形成对应分析图。
计算 B T B \mathbf{B^T}\mathbf{B} BTB 的特征值惯量表示相应维数对各类别的解释量最大维数 m min 16 − 1 , 6 − 1 m\min{16-1,6-1} mmin16−1,6−1 本例最多可以产生5个维数。从下表可看出第一维数的解释量达 77.4% 前两个维数的解释度已达92.1%。
维数奇异值惯量贡献率累积贡献率10.1898930.0360590.7737640.77376420.0828310.0068610.1472240.92098830.0471380.0022220.0476180.96866940.031130.0009690.0207950.98946450.0221590.0004910.0105361
行坐标
北京天津河北山西内蒙古辽宁第一维-0.07905-0.06783-0.263540.4577660.07715-0.13567第二维-0.03540.138818-0.10045-0.057150.156316-0.08455
吉林黑龙江上海江苏浙江安徽第一维-0.27126-0.197570.3868090.0869550.079122-0.14212第二维-0.000740.045985-0.07833-0.04222-0.01969-0.14225
福建江西山东河南第一维-0.17469-0.188590.069823-0.1462第二维-0.11317-0.15270.1003180.032858
列坐标
x1x2x3x4x5x6第一维-0.07905-0.06783-0.263540.4577660.07715-0.13567第二维-0.03540.138818-0.10045-0.057150.156316-0.08455
在下图中给出16个地区和6个指标在相同坐标系上绘制的散布图。 从图中可以看出地区和指标点可以分为两类
第一类包括指标点 x 4 , x 5 x_4,x_5 x4,x5 地区点为北京、天津、河北、上海、江苏、浙江、山东
第二类包括指标点 x 1 , x 2 , x 3 , x 6 x_1,x_2,x_3,x_6 x1,x2,x3,x6 地区点为其余地区。
3.程序
求解的MATLAB程序如下
clc, clear% 原始数据其中包含了16个地区的6项指标
originalData [190.33, 43.77, 9.73, 60.54, 49.01, 9.04; ...135.20, 36.40, 10.47, 44.16, 36.49, 3.94; ...95.21, 22.83, 9.30, 22.44, 22.81, 2.80; ...104.78, 25.11, 6.40, 9.89, 18.17, 3.25; ...128.41, 27.63, 8.94, 12.58, 23.99, 3.27; ...145.68, 32.83, 17.79, 27.29, 39.09, 3.47; ...159.37, 33.38, 18.37, 11.81, 25.29, 5.22; ...116.22, 29.57, 13.24, 13.76, 21.75, 6.04; ...221.11, 38.64, 12.53, 115.65, 50.82, 5.89; ...144.98, 29.12, 11.67, 42.60, 27.30, 5.74; ...169.92, 32.75, 12.72, 47.12, 34.35, 5.00; ...153.11, 23.09, 15.62, 23.54, 18.18, 6.39; ...144.92, 21.26, 16.96, 19.52, 21.75, 6.73; ...140.54, 21.50, 17.64, 19.19, 15.97, 4.94; ...115.84, 30.26, 12.20, 33.61, 33.77, 3.85; ...101.18, 23.26, 8.46, 20.20, 20.50, 4.30];% 计算原始数据的总和
totalSum sum(sum(originalData));% 计算原始数据的比例
ratioData originalData / totalSum;% 计算行和列的比例
rowRatio sum(ratioData, 2);
columnRatio sum(ratioData);% 计算行剖面数据
Row_prifile originalData ./ repmat(sum(originalData, 2), 1, size(originalData, 2));% 计算BB为对应分析的基础矩阵
B (ratioData - rowRatio * columnRatio) ./ sqrt(rowRatio*columnRatio);% 对B矩阵进行奇异值分解得到USV矩阵
[U, S, V] svd(B, econ);% 对V矩阵列求和并根据其符号调整权重
W1 sign(repmat(sum(V), size(V, 1), 1));% 对U矩阵列求和并根据其符号调整权重
W2 sign(repmat(sum(V), size(U, 1), 1));% 应用权重调整V和U矩阵
V_adjusted V .* W1;
U_adjusted U .* W2;% 计算lambda特征值的平方
lambda diag(S).^2;% 计算卡方统计量
chi2Square totalSum * (lambda);% 计算总的卡方统计量
totalChi2Square sum(chi2Square);% 计算贡献率
contributionRate lambda / sum(lambda);% 计算累计贡献率
cumulativeRate cumsum(contributionRate);% 计算行轮廓坐标
beta diag(rowRatio.^(-1 / 2)) * U_adjusted;
G beta * S% 计算列轮廓坐标
alpha diag(columnRatio.^(-1 / 2)) * V_adjusted;
F alpha * S% 计算样本点的个数
numOfSample size(G, 1);% 计算坐标的取值范围
range minmax(G(:, [1, 2]));% 画图略4.结果 详见上文分析过程地区和指标点可以分为两类
第一类包括指标点 x 4 , x 5 x_4,x_5 x4,x5 地区点为北京、天津、河北、上海、江苏、浙江、山东
第二类包括指标点 x 1 , x 2 , x 3 , x 6 x_1,x_2,x_3,x_6 x1,x2,x3,x6 地区点为其余地区。
习题10.62
1. 题目要求
同上。
2.解题过程——R型因子分析
解
对数据进行标准化处理。
计算变量间的相关系数得出相关矩阵 R \mathbf{R} R 然后计算初等载荷矩阵 Λ 1 \mathbf{\Lambda_1} Λ1 。
计算得到特征根与各因子的贡献如下表所示。
valuex1x2x3x4x5x6特征值3.558421.3162520.6082390.3733830.1071780.036527贡献率59.3070121.9375410.137326.2230521.7862950.608786累积贡献率59.3070181.2445491.3818797.6049299.39121100
选择 m ( m ≤ 6 ) m(m\leq6) m(m≤6) 个主因子构造因子模型 { x ~ 1 α 11 F ~ 1 α 12 F ~ 2 α 13 F ~ 3 ⋮ x ~ 7 α 71 F ~ 1 α 72 F ~ 2 α 73 F ~ 3 \begin{align*} \begin{cases} \widetilde{x}_1\alpha_{11}\widetilde{F}_1\alpha_{12}\widetilde{F}_2\alpha_{13}\widetilde{F}_3 \\ \vdots\\ \widetilde{x}_7\alpha_{71}\widetilde{F}_1\alpha_{72}\widetilde{F}_2\alpha_{73}\widetilde{F}_3 \end{cases}\end{align*} ⎩ ⎨ ⎧x 1α11F 1α12F 2α13F 3⋮x 7α71F 1α72F 2α73F 3 求得因子载荷等估计。
最终通过表格可以看出得到了3个因子第一个因子是穿住用因子第二个因子是燃料因子第3个因子是文化因子。
第(1)问中得到 x 4 , x 5 x_4,x_5 x4,x5 是一类变量这里得到 x 2 , x 4 , x 5 x_2,x_4,x_5 x2,x4,x5 是一类变量略有差异。
3.程序
求解的MATLAB程序如下
clc, clear % 原始数据包含了16个地区的6项指标
originalData [190.33, 43.77, 9.73, 60.54, 49.01, 9.04; ...135.20, 36.40, 10.47, 44.16, 36.49, 3.94; ...95.21, 22.83, 9.30, 22.44, 22.81, 2.80; ...104.78, 25.11, 6.40, 9.89, 18.17, 3.25; ...128.41, 27.63, 8.94, 12.58, 23.99, 3.27; ...145.68, 32.83, 17.79, 27.29, 39.09, 3.47; ...159.37, 33.38, 18.37, 11.81, 25.29, 5.22; ...116.22, 29.57, 13.24, 13.76, 21.75, 6.04; ...221.11, 38.64, 12.53, 115.65, 50.82, 5.89; ...144.98, 29.12, 11.67, 42.60, 27.30, 5.74; ...169.92, 32.75, 12.72, 47.12, 34.35, 5.00; ...153.11, 23.09, 15.62, 23.54, 18.18, 6.39; ...144.92, 21.26, 16.96, 19.52, 21.75, 6.73; ...140.54, 21.50, 17.64, 19.19, 15.97, 4.94; ...115.84, 30.26, 12.20, 33.61, 33.77, 3.85; ...101.18, 23.26, 8.46, 20.20, 20.50, 4.30];% 数据标准化
standardizedData zscore(originalData);% 计算相关性矩阵
correlationMatrix corrcoef(standardizedData);% 使用主成分分析方法对相关性矩阵进行处理得到特征向量、特征值和贡献率
[eigenvectors, eigenvalues, contribution] pcacov(correlationMatrix);% 计算累积贡献率
cumulativeContribution cumsum(contribution);% 根据特征向量的符号进行调整
adjustedSigns repmat(sign(sum(eigenvectors)), size(eigenvectors, 1), 1);
adjustedEigenvectors eigenvectors .* adjustedSigns;% 根据特征值进行缩放
scaledFactors repmat(sqrt(eigenvalues), size(adjustedEigenvectors, 1), 1);
scaledEigenvectors adjustedEigenvectors .* scaledFactors;% 计算贡献率
contribution1 sum(scaledEigenvectors.^2);% 选择的因子数量
factorNum 3;% 根据选择的因子数量得到对应的因子
selectedFactors scaledEigenvectors(:, [1:factorNum]);% 使用方差最大法进行因子旋转
[rotatedFactors, factorMatrix] rotatefactors(selectedFactors, method, varimax)% 合并旋转后的因子和其他因子
mergedFactors [rotatedFactors, scaledEigenvectors(:, [factorNum 1:end])];% 计算因子载荷量
factorLoads sum(rotatedFactors.^2, 2)% 计算贡献率
contribution2 sum(mergedFactors.^2)% 计算每个因子的贡献率
contributionRate contribution2(1:factorNum) / sum(contribution2);% 计算因子得分系数矩阵
factorScoreCoefficients inv(correlationMatrix) * rotatedFactors;4.结果 求得因子载荷等估计如下表所示。 可以看出得到了3个因子第一个因子是穿住用因子第二个因子是燃料因子第3个因子是文化因子。
第(1)问中得到 x 4 , x 5 x_4,x_5 x4,x5 是一类变量这里得到 x 2 , x 4 , x 5 x_2,x_4,x_5 x2,x4,x5 是一类变量略有差异。
习题10.63
1. 题目要求
同上。
2.解题过程——聚类分析
解
首先计算变量间的相关系数。用两变量 x j x_j xj与 x k x_k xk的相关系数作为它们的相似性度量,即 x j x_j xj与 x k x_k xk的相似系数为 r j k ∑ i 1 16 ( a i j − μ j ) ( a i k − μ k ) [ ∑ i 1 16 ( a i j − μ j ) 2 ∑ i 1 16 ( a i k − μ k ) 2 ] 1 / 2 , j , k 1 , … , 6. r_{jk} \frac{\sum_{i1}^{16}(a_{ij}-\mu_{j})(a_{ik}-\mu_k)} {[\sum_{i1}^{16}(a_{ij}-\mu_{j})^2\sum_{i1}^{16}(a_{ik}-\mu_k)^2]^{1/2}},j,k1,\dots,6. rjk[∑i116(aij−μj)2∑i116(aik−μk)2]1/2∑i116(aij−μj)(aik−μk),j,k1,…,6. 然后计算6个变量两两之间的距离构造距离矩阵。
接着使用最短距离法来测量类与类之间的距离,记类 G p G_{p} Gp和 G q G_{q} Gq之间的距离: D ( G p , G q ) min i ∈ G p , k ∈ G q d i k . D(G_p,G_q) \min_{i\in G_p,k\in G_q }{d_{ik}}. D(Gp,Gq)i∈Gp,k∈Gqmindik. 变量聚类的结果是变量 x 3 x_3 x3自成一类其他变量为一类。画出的变量聚类图如下图所示。 最后进行样本点聚类的 Q \mathbf{Q} Q型聚类分析。
计算16个样本点之间的两两马氏距离。
类与类间相似性度量。
画聚类图并对样本点进行分类。
样本点的聚类结果如下图所示。通过聚类图可以把地区分成4类北京自成一类吉林自成一类上海自成一类其他地区为一类。 3.程序
求解的MATLAB程序如下R型聚类
clc, clear % 原始数据包含了16个地区的6项指标
originalData [190.33, 43.77, 9.73, 60.54, 49.01, 9.04; ...135.20, 36.40, 10.47, 44.16, 36.49, 3.94; ...95.21, 22.83, 9.30, 22.44, 22.81, 2.80; ...104.78, 25.11, 6.40, 9.89, 18.17, 3.25; ...128.41, 27.63, 8.94, 12.58, 23.99, 3.27; ...145.68, 32.83, 17.79, 27.29, 39.09, 3.47; ...159.37, 33.38, 18.37, 11.81, 25.29, 5.22; ...116.22, 29.57, 13.24, 13.76, 21.75, 6.04; ...221.11, 38.64, 12.53, 115.65, 50.82, 5.89; ...144.98, 29.12, 11.67, 42.60, 27.30, 5.74; ...169.92, 32.75, 12.72, 47.12, 34.35, 5.00; ...153.11, 23.09, 15.62, 23.54, 18.18, 6.39; ...144.92, 21.26, 16.96, 19.52, 21.75, 6.73; ...140.54, 21.50, 17.64, 19.19, 15.97, 4.94; ...115.84, 30.26, 12.20, 33.61, 33.77, 3.85; ...101.18, 23.26, 8.46, 20.20, 20.50, 4.30];% 计算相关性矩阵
correlationMatrix corrcoef(originalData);% 计算距离矩阵
distanceMatrix 1 - abs(correlationMatrix);
distanceMatrix tril(distanceMatrix);% 将距离矩阵转为一维向量
distanceVector nonzeros(distanceMatrix);
distanceVector distanceVector;% 使用层次聚类算法进行聚类
linkageCluster linkage(distanceVector);% 选择最大聚类数量为2得到每个样本的类别
clusterLabels cluster(linkageCluster, maxclust, 2);% 找到属于第一类的样本
index1 find(clusterLabels 1);
index1 index1% 找到属于第二类的样本
index2 find(clusterLabels 2);
index2 index2% 生成树状图
h dendrogram(linkageCluster);% 设置树状图的颜色和线条宽度
set(h, Color, k, LineWidth, 1.3)求解的MATLAB程序如下Q型聚类
clc, clear % 原始数据包含了16个地区的6项指标
originalData [190.33, 43.77, 9.73, 60.54, 49.01, 9.04; ...135.20, 36.40, 10.47, 44.16, 36.49, 3.94; ...95.21, 22.83, 9.30, 22.44, 22.81, 2.80; ...104.78, 25.11, 6.40, 9.89, 18.17, 3.25; ...128.41, 27.63, 8.94, 12.58, 23.99, 3.27; ...145.68, 32.83, 17.79, 27.29, 39.09, 3.47; ...159.37, 33.38, 18.37, 11.81, 25.29, 5.22; ...116.22, 29.57, 13.24, 13.76, 21.75, 6.04; ...221.11, 38.64, 12.53, 115.65, 50.82, 5.89; ...144.98, 29.12, 11.67, 42.60, 27.30, 5.74; ...169.92, 32.75, 12.72, 47.12, 34.35, 5.00; ...153.11, 23.09, 15.62, 23.54, 18.18, 6.39; ...144.92, 21.26, 16.96, 19.52, 21.75, 6.73; ...140.54, 21.50, 17.64, 19.19, 15.97, 4.94; ...115.84, 30.26, 12.20, 33.61, 33.77, 3.85; ...101.18, 23.26, 8.46, 20.20, 20.50, 4.30];% 计算原始数据的协方差
covarianceMatrix cov(originalData);% 初始化距离矩阵
distanceMatrix zeros(size(originalData, 1));% 计算Q型距离
for j 1:15for i j1:16distanceMatrix(i,j) sqrt((originalData(i,:) - originalData(j,:)) * inv(covarianceMatrix) * (originalData(i,:) - originalData(j,:)));end
end% 将距离矩阵转为一维向量
distanceVector nonzeros(distanceMatrix);
distanceVector distanceVector;% 使用层次聚类算法进行聚类
linkageCluster linkage(distanceVector);% 生成树状图
dendro dendrogram(linkageCluster);% 设置树状图的颜色和线条宽度
set(dendro,Color,k,LineWidth,1.3)4.结果 变量 x 3 x_3 x3自成一类其他变量为一类。
北京自成一类吉林自成一类上海自成一类其他地区为一类。
如果这篇文章对你有帮助欢迎点赞与收藏~
