数学建模比赛知识P2

Matlab编程Part2

前文

评价方法大体分为主观赋权和客观赋权两类。
主观赋权：采取综合咨询评分，包括层次分析法等，所需数据较少。
客观赋权：根据指标间相关关系或指标变异程度确定权数，包括Topsis等。

1.层次分析法

层次分析法(AHP)是一种解决多目标的复杂问题的定性与定量相结合的决策分析方法。

步骤

1.建立层次结构模型，将决策问题分为3个或多个层次：
最高层：目标层，通常只有一个总目标。
中间层：准则层，指标层等，表示实现总目标的中间环节。
最低层：方案层，通常有几个方案可选。

2.构造判断（成对比较）矩阵，使用一致矩阵法：
因素之间两两相互比较，每层对比的因素不宜超过9个，判断矩阵元素的标度方法由Santy的1-9标度方法给出

得到以下判断矩阵：

3.层次单排序及其一致性检验
层次单排序：
获得判断矩阵对于最大特征值的特征向量经过归一化（向量中各元素之和为1）后的矩阵$W$。
注：归一化即为每一个数除以所有数的总和
一致性检验：
(1)求出判断矩阵的最大特征根 $\lambda$ 进而求出特征向量（权向量）$\omega$。
(2)由一致性指标公式求出一致性指标($CI$)
$$CI=\frac{\lambda -n}{n-1}$$
$CI=0$有完全的一致性，$CI$越大，不一致越严重。
(3)通过查表随机一致性指标($RI$)求出一致性比率($CR$)
$$CR=\frac{CI}{RI}$$
如果$CR<0.1$则层次单排序通过一致性检验。

4.层次总排序及其一致性检验
层次总排序：
计算某一层次所有因素对于最高层相对重要性的权值，从最高层到最底层依次进行。

$B$层第$i$个因素对总目标的权值为：$B_i={\sum }_{j=1}^m{a}_j{b}_{ij}$
一致性检验：
将权值转换成判断矩阵并归一化，求出一致性指标$CI$和随机一致性指标$RI$。
设$B$层$B_1,B_2,\cdots ,B_n$对上层（$A$层）中因素$A_j(j=1,2,\cdots ,m)$的层次单排序一致性指标为$C{I}_j$，随机一致性指标为$R{I}_j$，则层次总排序的一致性比率为
$$CR=\frac{a_{1}C{I}_1+\cdots +a_{m}C{I}_m}{a_{1}R{I}_1+\cdots +a_{m}R{I}_m}$$
如果$CR<0.1$则层次总排序通过一致性检验。到此，根据最低层的层次总排序做出决策。

编程实现

这是层次单排序及一致性检验的Matlab代码，可以根据题目进行修改优化为所需代码。

disp('请输入准则层判断矩阵A(n阶)');
A=input('A=');
[n,n]=size(A);
[V,D]=eig(A);%求得特征向量和特征值
            %求出最大特征值和它所对应的特征向量
tempNum=D(1,1);
pos=1;
for h=1:n
    if D(h,h)>tempNum
        tempNum=D(h,h);
        pos=h;
    end
end    
w=abs(V(:,pos));
w=w/sum(w);
t=D(pos,pos);
disp('准则层特征向量w=');disp(w);disp('准则层最大特征根t=');disp(t);
%以下是一致性检验
CI=(t-n)/(n-1);RI=[0 0 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 1.60 1.61 1.615 1.62 1.63];
CR=CI/RI(n);
if CR<0.10
    disp('此矩阵的一致性可以接受!');
    disp('CI=');disp(CI);
    disp('CR=');disp(CR);
else disp('此矩阵的一致性验证失败，请重新进行评分!');
end

---

2.Topsis（理想解法）

Topsis是一种多指标评价方法，构造正理想解和负理想解，即最优解和最劣解，通过计算方案到理想方案的贴近度对方案进行排序。

步骤

1.求出规范决策矩阵
设多属性决策问题的决策矩阵为$A=(a_{ij}{)}_{m\times n}$，变换后的决策矩阵为$B=(b_{ij}{)}_{m\times n}$
(1)数据预处理，把多种类型的数据，如成本型，效益型，区间型等，转化效益型，即数据越大越好（或转化成本型，数据越小越好）。
线性变换：a.效益型数据变换 b.成本型转效益型
a：$b_{ij}=a_{ij}/a_j^{max}$
b：$b_{ij}=1-a_{ij}/a_j^{max}$
区间型数据转效益型数据：设最优数据区间为$[a_j^0,a_j^{\ast }]$，$a_j^1$为无法容忍的下限，$a_j^2$为无法容忍的上限
$$b_{ij}=\{\begin{array}{c}1-(a_j^0-a_{ij})/(a_j^0-a_j^1),a_j^1⩽a_{ij}<a_j^0\\ 1,a_j^0⩽a_{ij}⩽a_j^{\ast }\\ 1-(a_{ij}-a_j^{\ast })/(a_j^2-a_j^{\ast }),a_j^{\ast }<a_{ij}⩽a_j^2\\ 0,\text{其它}\end{array}$$
(2)数据规范化，消除量纲的影响。
$$b_{ij}=\frac{a_{ij}}{\sqrt{{\sum }_{i=1}^m{a}_{ij}^2}}(i=1,\cdots ,m,j=1,\cdots ,n)$$
拓展：
归一化：a.效益型归一化 b.成本型归一化
a：$b_{ij}=a_{ij}-a_j^{min}/a_j^{max}-a_j^{min}$
b：$b_{ij}=a_j^{max}-a_{ij}/a_j^{max}-a_j^{min}$
2.求出加权规范阵
设各属性的权重向量为$w=[w_1,w_2,\cdots ,w_n{]}^T$，则加权规范阵为$C=(c_{ij}{)}_{m\times n}$，其中
$$c_{ij}=w_j\cdot b_{ij}(i=1,\cdots ,m,j=1,\cdots ,n)$$
3.确定正理想解$C^{\ast }$和负理想解$C^0$
找到每组数据中的最优和最劣的值分别组合成正理想解和负理想解。

4.计算各方案的排队指标值并排序
备选方案$d_i$到正理想解的距离为
$$s_i^{\ast }=\sqrt{\sum _{j=1}^n(c_{ij}-c_j^{\ast }{)}^2}(i=1,\cdots ,m)$$
备选方案$d_i$到负理想解的距离为
$$s_i^0=\sqrt{\sum _{j=1}^n(c_{ij}-c_j^0{)}^2}(i=1,\cdots ,m)$$
排队指标值
$$f_i^{\ast }=\frac{s_i^0}{s_i^0+s_i^{\ast }}(i=1,\cdots ,m)$$
按$f_i^{\ast }$由大到小排列方案的优劣次序。

编程实现

Matlab编程示例，不包括区间型数据的转换

x=[
21584   76.7    7.3 1.01    78.3    97.5    2.0
24372   86.3    7.4 0.80    91.1    98.0    2.0
22041   81.8    7.3 0.62    91.1    97.3    3.2
21115   84.5    6.9 0.60    90.2    97.7    2.9
24633   90.3    6.9 0.25    95.5    97.9    3.6];%矩阵
[n,m]=size(x);
%将3，4，7的低优指标去倒数转化为高优指标并且把所有指标换成接近的大小
x(:,1)=x(:,1)/100;
x(:,3)=(1./x(:,3))*100;
x(:,4)=(1./x(:,4))*100;
x(:,7)=(1./x(:,7))*100;
zh=zeros(1,m);
d1=zeros(1,n); %最小值矩阵
d2=zeros(1,n); %最大值矩阵
c=zeros(1,n);  %接近程度
%数据规范化
for i=1:m
    for j=1:n
        zh(i)=zh(i)+x(j,i)^2;
    end
end
for i=1:m
    for j=1:n
       x(j,i)=x(j,i)/sqrt( zh(i));
    end
end
%计算加权
w=[0.2 0.1 0.1 0.2 0.15 0.15 0.1];
x=x.*repmat(w,n,1);
%计算距离
xx=min(x);
dd=max(x);
for i=1:n
    for j=1:m
        d1(i)=d1(i)+(x(i,j)-xx(j))^2;
    end
    d1(i)=sqrt(d1(i));
end
for i=1:n
    for j=1:m
        d2(i)=d2(i)+(x(i,j)-dd(j))^2;
    end
    d2(i)=sqrt(d2(i));
end
%计算接近程度
for i=1:n
    c(i)=d1(i)/(d2(i)+d1(i));
end
%排序结果
[sf,ind]=sort(c,'descend');
%'ascend'升序排序，'descend'降序排序。

---

3.聚类分析法

聚类是将相似的数据划分为一组或一类，且利用各数据间的距离表示相似程度。
聚类分析包括对样品的分类（Q型，主要部分）和对变量的分类（R型，一般用来降维，减少变量个数）
(1)样品相似度——距离的定义
设有n个样品的p元数据，每个样品可看作p元空间的一个点，两点间距离记为$d(x_i,x_j)$且满足：
$d(x_i,x_j)⩾0，d(x_i,x_j)=d(x_j,x_i)，d(x_i,x_j)⩽d(x_i,x_k)+d(x_k,x_j)$
常用距离：欧氏距离：$d(x_i,x_j)=\sqrt{{\sum }_{k=1}^p(x_{ik}-x_{jk}{)}^2}$
对应Matlab命令：$pdist(x);$

(2)类间距离
设$d_{ij}$表示两个样品$x_i,x_j$之间的距离，$G_p,G_q$表示两个类，两类之间的距离通常用两类中样品之间的最短距离表示，即$D_{pq}=\min d_{ij},i\in G_p,j\in G_q$

(3)变量相似度——相似系数的定义
设$C_{a,b}$表示变量之间的相似系数，则满足：
$|C_{ab}|⩽1，C_{aa}=1，C_{ab}=C_{ba}$
相似系数常用相关系数和夹角余弦表示。

谱系聚类法

步骤
1.定义距离并计算
2.取距离最小的两个样品合为一个类，并计算新的距离，循环操作，直到合为一个大类
3.画出聚类图
4.决定类的个数和分类
编程实现
可以一步一步的运行程序，方便观察

x=[7.9 39.77 8.49 12.94 19.27 11.05 2.04 13.29
7.68 50.37 11.35 13.3 19.25 14.59 2.75 14.87
9.42 27.93 8.2 8.14 16.17 9.42 1.55 9.76
9.16 27.98 9.01 9.32 15.99 9.1 1.82 11.35
10.06 28.64 10.52 10.05 16.18 8.39 1.96 10.81];
y=pdist(x);%计算欧氏距离
z=linkage(y);%计算类间距离，并循环
H=dendrogram(z);%画聚类图
T=cluster(z,3);%把样品分为3类

K-平均聚类算法(k-means)

步骤
1.随机设定确定个数的聚类中心
2.样品划分到最近的类中
3.求每个类中样品的平均数并更新聚类中心，并循环进行2，3步
4.画出聚类结果
编程实现

function [Idx, Center] = K_means(X, xstart)
% K-means聚类
% Idx是数据点属于哪个类的标记，Center是每个类的中心位置
% X是全部二维数据点，xstart是类的初始中心位置

len = length(X);        %X中的数据点个数
Idx = zeros(len, 1);    %每个数据点的Id，即属于哪个类

C1 = xstart(1,:);       %第1类的中心位置
C2 = xstart(2,:);       %第2类的中心位置
C3 = xstart(3,:);       %第3类的中心位置

for i_for = 1:100
    %为避免循环运行时间过长，通常设置一个循环次数
    %或相邻两次聚类中心位置调整幅度小于某阈值则停止
    
    %更新数据点属于哪个类
    for i = 1:len
        x_temp = X(i,:);    %提取出单个数据点
        d1 = norm(x_temp - C1);    %与第1个类的距离
        d2 = norm(x_temp - C2);    %与第2个类的距离
        d3 = norm(x_temp - C3);    %与第3个类的距离
        d = [d1;d2;d3];
        [~, id] = min(d);   %离哪个类最近则属于那个类
        Idx(i) = id;
    end
    
    %更新类的中心位置
    L1 = X(Idx == 1,:);     %属于第1类的数据点
    L2 = X(Idx == 2,:);     %属于第2类的数据点
    L3 = X(Idx == 3,:);     %属于第3类的数据点
    C1 = mean(L1);      %更新第1类的中心位置
    C2 = mean(L2);      %更新第2类的中心位置
    C3 = mean(L3);      %更新第3类的中心位置
end
Center = [C1; C2; C3];  %类的中心位置

%演示数据
%% 1 random sample
%随机生成三组数据
a = rand(30,2) * 2;
b = rand(30,2) * 5;
c = rand(30,2) * 10;
figure(1);
subplot(2,2,1); 
plot(a(:,1), a(:,2), 'r.'); hold on
plot(b(:,1), b(:,2), 'g*');
plot(c(:,1), c(:,2), 'bx'); hold off
grid on;
title('raw data');

%% 2 K-means cluster
X = [a; b; c];  %需要聚类的数据点
xstart = [2 2; 5 5; 8 8];  %初始聚类中心
subplot(2,2,2);
plot(X(:,1), X(:,2), 'kx'); hold on
plot(xstart(:,1), xstart(:,2), 'r*'); hold off
grid on;
title('raw data center');

[Idx, Center] = K_means(X, xstart);
subplot(2,2,4);
plot(X(Idx==1,1), X(Idx==1,2), 'kx'); hold on
plot(X(Idx==2,1), X(Idx==2,2), 'gx');
plot(X(Idx==3,1), X(Idx==3,2), 'bx');
plot(Center(:,1), Center(:,2), 'r*'); hold off
grid on;
title('K-means cluster result');

disp('xstart = ');
disp(xstart);
disp('Center = ');
disp(Center);