当我们对事物进行评价时,通常会使用综合评价模型来综合考虑各个方面的因素。对于数学建模而言,综合评价模型是一种系统化的方法,通过将多个评价指标进行加权和综合,得出对事物整体的评价结果。它可以应用于各种领域和场景,例如产品评价、学生综合素质评价、企业绩效评估等。

层次分析法

层次分析法(Analytic Hierarchy Process,AHP)是一种用于进行多指标综合评价的方法,它可以帮助我们确定各个指标的权重,并最终得出综合评价结果。这是最基础的模型之一,其主要用于解决评价类问题,例如:选择哪种方案最好、哪位运动员或者员工表现的更优秀。

层次分析法的基本思想是将一个复杂的问题分解成若干个层次,从而形成一个层次结构。每一层次包含若干个指标或因素,在进行评价时,我们需要对各个指标进行比较和权重计算,以确定它们在整个层次结构中的相对重要程度。假设我们要评估三种手机品牌(A、B、C)的性能,我们选择了以下指标:屏幕质量、摄像头性能和电池续航。我们需要对这三个指标进行两两比较,并根据比较结果给出权重值。

首先,你要确定一些重要的指标,比如屏幕质量、摄像头性能和电池续航。然后,你需要比较这些指标之间的相对重要程度。比如,你可能认为摄像头性能对你来说更重要,因为你经常拍照或者视频通话。接下来,你需要对每个指标进行进一步的比较。比如,你要比较不同手机的屏幕质量。你可以列出几个备选手机,并给它们打分。然后,你可以将这些分数转化为相对权重,用来表示在整个层次结构中的重要程度。最后,根据这些权重来做出决策。比如,如果一个手机在屏幕质量方面得分很高,并且屏幕质量在整个层次结构中的权重也很高,那么你可能会倾向于选择这款手机。

需要注意的是,层次分析法通过构建一个判断矩阵来进行两两比较,从而得到各个指标之间的权重值。在进行两两比较时,我们需要确保判断矩阵是一致的,即各个比较的结果相互符合。判断矩阵有时可能会出现矛盾情况。这种矛盾可能是由于主观判断的不一致或者信息不完整所导致的,这种情况很容易出现且很难被直接察觉到。以下是一个简化的例子:

假设我们的判断矩阵如下所示:

1
2
3
4
5
        屏幕质量     摄像头性能   电池续航
屏幕质量    1         3         2
摄像头性能   1/3      1         1/2
电池续航    1/2       2         1

从判断矩阵可以看出,我们认为屏幕质量比摄像头性能好3倍,而摄像头性能比电池续航好2倍,同时电池续航比屏幕质量好1/2倍。这样的比较存在矛盾,判断矩阵存在较大的一致性问题。在这种情况下,我们需要重新检查和修改判断矩阵,可能需要重新进行比较并调整权重值,以消除矛盾并提高一致性。这样做可以确保判断矩阵符合层次分析法的原则和要求,从而得到更准确和可靠的评价结果。

为了解决这个问题,我们需要进行一致性检验。

一致性检验是确保我们所构建的判断矩阵是合理和可靠的关键步骤。一致性检验的目的是验证判断矩阵中的数字是否具有内在的一致性。如果判断矩阵不够一致,那么权重计算可能会失真,从而影响最终的评价结果。在进行一致性检验时,我们使用一致性指标(Consistency Index, CI)和一致性比率(Consistency Ratio, CR)来评估判断矩阵的一致性程度。CI是通过计算判断矩阵的最大特征值与矩阵大小的差异得出的,而CR则是将CI与一个随机一致性指标进行比较。如果CR的值小于某个预定的阈值(通常为0.1),那么判断矩阵可以被认为是一致的。如果CR的值大于阈值,那么判断矩阵存在较大的一致性问题,需要重新进行比较和修改。

在通过一致性检验后,根据每个指标的权重值,将各个指标的得分进行加权求和,得到最终的综合评价结果。

完整代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
disp('请输入判断矩阵A')
A = input('A = ');

[n, n] = size(A);

% 方法1:算术平均法求权重
Sum_A = sum(A);
SUM_A = repmat(Sum_A, n, 1);
Stand_A = A ./ SUM_A;

disp('算术平均法求权重的结果为:');
disp(sum(Stand_A, 2) ./ n);

%几何平均法求权重
Product_A = prod(A, 2);
Product_n_A = Product_A .^ (1/n);
disp('几何平均法求权重的结果为:');
disp(Product_n_A ./ sum(Product_n_A));

% 特征值法求权重
[V, D] = eig(A);
Max_eig = max(max(D));
[r, c] = find(D == Max_eig, 1);
disp('特征值法求权重的结果为:');
disp(V(:, c) ./ sum(V(:, c)));

% 计算一致性比例CR
CI = (Max_eig - n) / (n - 1);
RI = [0 0.0001 0.52 0.89 1.12 1.26 1.36 1.41 1.46 1.49 1.52 1.54 1.56 1.58 1.59];  % RI理论上最多支持 n = 15

CR = CI / RI(n);
disp('一致性指标CI = ');
disp(CI);
disp('一致性比例CR = ');
disp(CR);

% 根据CR的值判断一致性是否可接受
if CR < 0.10
    disp('因为CR < 0.10,一致性检验通过。');
else
    disp('CR >= 0.10,该判断矩阵A需要修改。');
end

TOPSIS法(优劣解距离法)

当数据已知时,TOPSIS法是一种更常用的综合评价方法,它能充分利用原始数据的信息,更精确地反映各评价方案之间的差距。

层次分析法存在一些小问题:

  1. 评价的决策层不能太多,太多的话n会很大,所以判断矩阵和一致矩阵差异可能会很大。
  2. 权重是我们自己打分定的,如果决策层中指标的数据是已知的,应该更好地利用这些数据来使得评价的更加准确。

这些数据可能符合以下的指标:

  • 极大型(效益型)指标:越大(多)越好,比如成绩、GDP增速、企业利润
  • 极小型(成本型)指标:越小(少)越好,比如费用、坏品率、污染程度
  • 中间型指标:越接近某个值越好,比如水质量评估时的PH值
  • 区间型指标:落在某个区间最好,比如体温、水中植物性营养物量

基本过程为,先将原始数据矩阵统一指标类型。最常用的是将所有的指标转化为极大型,称为指标正向化。因为不同指标有不同的评价方式,为了方便比较,需要将不同指标的得分转化成相同类型的数值,一般会将指标值按照一定的规则进行转换,使所有指标都能够在同样的范围内进行比较。

由于不同指标的量纲不同,对评价结果的影响也不同,所以需要对所有指标的数据进行标准化处理,以消除各指标量纲的影响。标准化后,数据在不同指标上的权重是相同的,可以直接进行比较。具体数学过程略。

处理完成后,就可以找到有限方案中的最优方案和最劣方案,然后分别计算各评价对象与最优方案和最劣方案间的距离,获得各评价对象与最优方案的相对接近程度,以此作为评价优劣的依据。该方法对数据分布及样本含量没有严格限制,数据计算简单易行。

代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
%%  第一步:把数据复制到工作区,并将这个矩阵命名为X
clear;clc
load matlab.mat

%%  第二步:判断是否需要正向化
[n,m] = size(X);
disp(['共有' num2str(n) '个评价对象, ' num2str(m) '个评价指标'])
Judge = input(['这' num2str(m) '个指标是否需要经过正向化处理,需要请输入1 ,不需要输入0:  ']);
if Judge == 1
    Position = input('请输入需要正向化处理的指标所在的列,例如第2、3、6三列需要处理,那么你需要输入[2,3,6]: '); %[2,3,4]
    disp('请输入需要处理的这些列的指标类型(1:极小型, 2:中间型, 3:区间型) ')
    Type = input('例如:第2列是极小型,第3列是区间型,第6列是中间型,就输入[1,3,2]:  '); %[2,1,3]
    for i = 1 : size(Position,2)  
        X(:,Position(i)) = Positivization(X(:,Position(i)),Type(i),Position(i));
    % Positivization是我们自己定义的函数,其作用是进行正向化,其一共接收三个参数
    % 第一个参数是要正向化处理的那一列向量 X(:,Position(i))   回顾上一讲的知识,X(:,n)表示取第n列的全部元素
    % 第二个参数是对应的这一列的指标类型(1:极小型, 2:中间型, 3:区间型)
    % 第三个参数是告诉函数我们正在处理的是原始矩阵中的哪一列
    % 该函数有一个返回值,它返回正向化之后的指标,我们可以将其直接赋值给我们原始要处理的那一列向量
    end
    disp('正向化后的矩阵 X =  ')
    disp(X)
end

%% 第三步:对正向化后的矩阵进行标准化
Z = X ./ repmat(sum(X.*X) .^ 0.5, n, 1);
disp('标准化矩阵 Z = ')
disp(Z)

%% 第四步:计算与最大值的距离和最小值的距离,并算出得分
D_P = sum([(Z - repmat(max(Z),n,1)) .^ 2 ],2) .^ 0.5;   % D+ 与最大值的距离向量
D_N = sum([(Z - repmat(min(Z),n,1)) .^ 2 ],2) .^ 0.5;   % D- 与最小值的距离向量
S = D_N ./ (D_P+D_N);    % 未归一化的得分
disp('最后得分:')
stand_S = S / sum(S)
[sorted_S,index] = sort(stand_S ,'descend')

聚类模型

聚类,就是将样本划分为由类似的对象组成的多个类的过程,这里的“类”是事先未知的。聚类后,我们可以更加准确的在每个类中单独使用统计模型进行估计、分析或预测;也可以探究不同类之间的相关性和主要差异。

K-means聚类

K-means聚类的算法流程:

  1. 指定需要划分的簇的个数K值(类的个数);
  2. 随机地选择K个数据对象作为初始的聚类中心(不一定要是我们的样本点);
  3. 计算其余的各个数据对象到这K个初始聚类中心的距离,把数据对象划归到距离它最近的那个中心所处在的簇类中;
  4. 调整新类并且重新计算出新类的中心;
  5. 循环步骤三和四,看中心是否收敛(不变),如果收敛或达到迭代次数则停止循环;

优点:

  1. 算法简单、快速。
  2. 对处理大数据集,该算法是相对高效率的。

缺点:

  1. 要求用户必须事先给出要生成的簇的数目K。
  2. 对初值敏感。
  3. 对于孤立点数据敏感。

K-means++算法可解决2和3这两个缺点,只对K-means算法“初始化K个聚类中心” 这一步进行了优化。k-means++算法选择初始聚类中心的基本原则是:初始的聚类中心之间的相互距离要尽可能的远。

算法描述如下:

  1. 随机选取一个样本作为第一个聚类中心;
  2. 计算每个样本与当前已有聚类中心的最短距离(即与最近一个聚类中心的距离),这个值越大,表示被选取作为聚类中心的概率较大;最后,用轮盘法(依据概率大小来进行抽选)选出下一个聚类中心;
  3. 重复步骤二,直到选出K个聚类中心。选出初始点后,就继续使用标准的K-means算法了。

这一算法的最大缺点是,簇的个数k必须事前给定。

系统(层次)聚类

系统聚类的合并算法通过计算两类数据点间的距离,对最为接近的两类数据点进行组合,并反复迭代这一过程,直到将所有数据点合成一类,并生成聚类谱系图。

系统(层次)聚类的算法流程:

  1. 将每个对象看作一类,计算两两之间的最小距离;
  2. 将距离最小的两个类合并成一个新类;
  3. 重新计算新类与所有类之间的距离;
  4. 重复二三两步,直到所有类最后合并成一类;
  5. 结束。

DBSCAN算法

聚类前不需要预先指定聚类的个数,生成的簇的个数不定(和数据有关)。该算法利用基于密度的聚类的概念,即要求聚类空间中的一定区域内所包含对象(点或其他空间对象)的数目不小于某一给定阈值。该方法能在具有噪声的空间数据库中发现任意形状的簇,可将密度足够大的相邻区域连接,能有效处理异常数据。

DBSCAN算法将数据点分为三类:

  • 核心点:在半径Eps内含有不少于MinPts数目的点
  • 边界点:在半径Eps内点的数量小于MinPts,但是落在核心点的邻域内
  • 噪音点:既不是核心点也不是边界点的点

只有两个指标,且你做出散点图后发现数据表现得很“DBSCAN”,这时候再用DBSCAN进行聚类。

其他情况下,全部使用系统聚类吧!

具体实现使用spss,详见帮助文档

主成分分析

主成分分析(Principal Component Analysis,PCA),是一种降维算法,它能将多个指标转换为少数几个主成分,这些主成分是原始变量的线性组合,且彼此之间互不相关,其能反映出原始数据的大部分信息。一般来说,当研究的问题涉及到多变量且变量之间存在很强的相关性时,我们可考虑使用主成分分析的方法来对数据进行简化。

在实际问题研究中,多变量问题是经常会遇到的。变量太多,会增加分析问题的难度与复杂性,而且在许多实际问题中,多个变量之间是具有一定的相关关系的。因此,人们会很自然地想到,能否在相关分析的基础上,用较少的新变量代替原来较多的旧变量,而且使这些较少的新变量尽可能多地保留原来变量所反映的信息?事实上,这种想法是可以实现的,主成分分析方法就是综合处理这种问题的一种强有力的工具。主成分分析是把原来多个变量划为少数几个综合指标的一种统计分析方法。从数学角度来看,这是一种降维处理技术。

降维是将高维度的数据(指标太多)保留下最重要的一些特征,去除噪声和不重要的特征,从而实现提升数据处理速度的目的。在实际的生产和应用中,降维在一定的信息损失范围内,可以为我们节省大量的时间和成本。降维也成为应用非常广泛的数据预处理方法。

主成分分析可以用于聚类和回归问题。

在聚类问题中,主成分分析可以用来降低数据的维度。通过PCA,我们可以将原始数据转换为一组新的主成分(即特征向量),这些主成分是原始数据中变化最大的方向。然后,我们可以使用聚类算法(如K均值聚类、层次聚类等)来对这些主成分进行聚类,以实现样本的分组。

在回归问题中,主成分分析可以用来处理多重共线性。多重共线性指的是自变量之间存在高度相关性的情况,在回归分析中可能导致模型不稳定或系数估计不准确。通过应用主成分分析,我们可以将相关性较高的自变量转换为一组无关的主成分,并利用这些主成分进行回归分析,以建立更稳定和准确的预测模型。

主成分分析在聚类和回归问题中的具体应用方法会有所不同。在聚类问题中,我们可以根据主成分的方差解释率选择保留的主成分数量,然后再应用聚类算法进行分组。在回归问题中,我们可以根据主成分的贡献程度选择保留的主成分数量,然后再进行回归分析。

代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
clear;clc
load data1.mat   % 主成分聚类

[n,p] = size(x);  % n是样本个数,p是指标个数

%% 第一步:对数据x标准化为X
X=zscore(x);   % matlab内置的标准化函数(x-mean(x))/std(x)

%% 第二步:计算样本协方差矩阵
R = cov(X);

%% 注意:以上两步可合并为下面一步:直接计算样本相关系数矩阵
R = corrcoef(x);
disp('样本相关系数矩阵为:')
disp(R)

%% 第三步:计算R的特征值和特征向量
% 注意:R是半正定矩阵,所以其特征值不为负数
% R同时是对称矩阵,Matlab计算对称矩阵时,会将特征值按照从小到大排列
[V,D] = eig(R);  % V 特征向量矩阵  D 特征值构成的对角矩阵

%% 第四步:计算主成分贡献率和累计贡献率
lambda = diag(D);  % diag函数用于得到一个矩阵的主对角线元素值(返回的是列向量)
lambda = lambda(end:-1:1);  % 因为lambda向量是从小大到排序的,我们将其调个头
contribution_rate = lambda / sum(lambda);  % 计算贡献率
cum_contribution_rate = cumsum(lambda)/ sum(lambda);   % 计算累计贡献率  cumsum是求累加值的函数
disp('特征值为:')
disp(lambda')  % 转置为行向量,方便展示
disp('贡献率为:')
disp(contribution_rate')
disp('累计贡献率为:')
disp(cum_contribution_rate')
disp('与特征值对应的特征向量矩阵为:')
% 注意:这里的特征向量要和特征值一一对应,之前特征值相当于颠倒过来了,因此特征向量的各列需要颠倒过来
%  rot90函数可以使一个矩阵逆时针旋转90度,然后再转置,就可以实现将矩阵的列颠倒的效果
V=rot90(V)';
disp(V)

%% 计算我们所需要的主成分的值
m =input('请输入需要保存的主成分的个数:  ');
F = zeros(n,m);  %初始化保存主成分的矩阵(每一列是一个主成分)
for i = 1:m
    ai = V(:,i)';   % 将第i个特征向量取出,并转置为行向量
    Ai = repmat(ai,n,1);   % 将这个行向量重复n次,构成一个n*p的矩阵
    F(:, i) = sum(Ai .* X, 2);  % 注意,对标准化的数据求了权重后要计算每一行的和
end

%% (1)主成分聚类 :将主成分指标所在的F矩阵复制到Excel表格,然后再用Spss进行聚类
% 在Excel第一行输入指标名称(F1,F2, ..., Fm)
% 双击Matlab工作区的F,进入变量编辑中,然后复制里面的数据到Excel表格
% 导出数据之后,我们后续的分析就可以在Spss中进行。

%%(2)主成分回归:将x使用主成分得到主成分指标,并将y标准化,接着导出到Excel,然后再使用Stata回归
% Y = zscore(y);  % 一定要将y进行标准化哦~
% 在Excel第一行输入指标名称(Y,F1, F2, ..., Fm)
% 分别双击Matlab工作区的Y和F,进入变量编辑中,然后复制里面的数据到Excel表格
% 导出数据之后,我们后续的分析就可以在Stata中进行。

灰色关联分析

数理统计中的回归分析、方差分析、主成分分析等都是用来进行系统分析的方法。这些方法都有下述不足之处:

  1. 要求有大量数据,数据量少就难以找出统计规律。
  2. 要求样本服从某个典型的概率分布,要求各因素数据与系统特征数据之间呈线性关系且各因素之间彼此无关,这种要求往往难以满足。
  3. 可能出现量化结果与定性分析结果不符的现象,导致系统的关系和规律遭到歪曲和颠倒。

对一个抽象的系统或现象进行分析,首先要选准反映系统行为特征的数据序列,称为找系统行为的映射量,用映射量来间接地表征系统行为。例如,用国民平均接受教育的年数来反映教育发达程度,用刑事案件的发案率来反映社会治安面貌和社会秩序,用医院挂号次数来反映国民的健康水平等。有了系统行为特征数据和相关因素的数据,即可作出各个序列的图形,从直观上进行分析。

灰色关联分析的基本思想是根据序列曲线几何形状的相似程度来判断其联系是否紧密。曲线越接近,相应序列之间的关联度就越大,反之就越小。

这不是传统的数理统计方法,但是可以用。

分析代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
%% 灰色关联分析用于系统分析
clear;clc
load gdp.mat  % 导入数据 一个6*4的矩阵

Mean = mean(gdp);  % 求出每一列的均值以供后续的数据预处理
gdp = gdp ./ repmat(Mean,size(gdp,1),1);  %size(gdp,1)=6, repmat(Mean,6,1)可以将矩阵进行复制,复制为和gdp同等大小,然后使用点除(对应元素相除),这
disp('预处理后的矩阵为:'); disp(gdp)
Y = gdp(:,1);  % 母序列
X = gdp(:,2:end); % 子序列
absX0_Xi = abs(X - repmat(Y,1,size(X,2)))  % 计算|X0-Xi|矩阵(在这里我们把X0定义为了Y)
a = min(min(absX0_Xi))    % 计算两级最小差a
b = max(max(absX0_Xi))  % 计算两级最大差b
rho = 0.5; % 分辨系数取0.5
gamma = (a+rho*b) ./ (absX0_Xi  + rho*b)  % 计算子序列中各个指标与母序列的关联系数
disp('子序列中各个指标的灰色关联度分别为:')
disp(mean(gamma))

评价代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
%% 灰色关联分析用于综合评价模型例
clear;clc
load data_water_quality.mat

%%  判断是否需要正向化
[n,m] = size(X);
disp(['共有' num2str(n) '个评价对象, ' num2str(m) '个评价指标'])
Judge = input(['这' num2str(m) '个指标是否需要经过正向化处理,需要请输入1 ,不需要输入0:  ']);   %1

if Judge == 1
    Position = input('请输入需要正向化处理的指标所在的列,例如第2、3、6三列需要处理,那么你需要输入[2,3,6]: '); %[2,3,4]
    disp('请输入需要处理的这些列的指  标类型(1:极小型, 2:中间型, 3:区间型) ')
    Type = input('例如:第2列是极小型,第3列是区间型,第6列是中间型,就输入[1,3,2]:  '); %[2,1,3]
    % 注意,Position和Type是两个同维度的行向量
    for i = 1 : size(Position,2)  %这里需要对这些列分别处理,因此我们需要知道一共要处理的次数,即循环的次数
        X(:,Position(i)) = Positivization(X(:,Position(i)),Type(i),Position(i));
    % Positivization是我们自己定义的函数,其作用是进行正向化,其一共接收三个参数
    end
    disp('正向化后的矩阵 X =  ')
    disp(X)
end

%% 对正向化后的矩阵进行预处理
Mean = mean(X);  % 求出每一列的均值以供后续的数据预处理
Z = X ./ repmat(Mean,size(X,1),1);  
disp('预处理后的矩阵为:'); disp(Z)

%% 构造母序列和子序列
Y = max(Z,[],2);  % 母序列为虚拟的,用每一行的最大值构成的列向量表示母序列
X = Z; % 子序列就是预处理后的数据矩阵

%% 计算得分
absX0_Xi = abs(X - repmat(Y,1,size(X,2)))  % 计算|X0-Xi|矩阵
a = min(min(absX0_Xi))    % 计算两级最小差a
b = max(max(absX0_Xi))  % 计算两级最大差b
rho = 0.5; % 分辨系数取0.5
gamma = (a+rho*b) ./ (absX0_Xi  + rho*b)  % 计算子序列中各个指标与母序列的关联系数
weight = mean(gamma) / sum(mean(gamma));  % 利用子序列中各个指标的灰色关联度计算权重
score = sum(X .* repmat(weight,size(X,1),1),2);   % 未归一化的得分
stand_S = score / sum(score);   % 归一化后的得分
[sorted_S,index] = sort(stand_S ,'descend') % 进行排序