8. 层次分析模型

数学建模

张瑞
中国科学技术大学数学科学学院

综合评价模型

综合评价是使用比较系统的、规范的方法对于多个指标、多个单位同时进行评价的方法，也叫多指标综合评价方法。

  • 综合评价的应用领域和范围非常广泛。比如，环境监测综合评价、药物临床试验综合评价、社会治安综合评价，综合经济效益评价、房地产市场景气程度综合评价等等

评价方法分为两类：

  • 基于主观赋权的方法。如，综合指数法，层次分析法
  • 基于客观赋权，根据各指标间相关关系或各指标值变异程度来确定权数。如，主成份分析、因子分析法、理想解法（TOPSIS）

层次分析模型

例 1. 想去旅游。在黄山、桂林、北戴河中，综合考虑景色、费用、居住条件、饮食、旅途等因素，选一个为目的地。

  • 作比较判断时人的主观选择起相当大的作用,各因素的重要性难以量化。
  • 人们在进行社会的、经济的以及科学管理领域问题的系统分析中,面临的常常是一个由相互关联、相互制约的众多因素构成的复杂而往往缺少定量数据的系统。
  • 层次分析法(Analytic Hierarchy Process,简称 AHP)是对一些较为复杂、较为模糊的问题作出决策的简易方法,它特别适用于那些难于完全定量分析的问题。
  • 它是美国运筹学家 T. L. Saaty 教授于上世纪 70 年代初期提出的一种简便、灵活而又实用的多准则决策方法。
  • AHP是一种定性与定量相结合的、系统化、层次化的分析方法

运用层次分析法建模，大体上可按下面四个步骤进行：

  1. 建立递阶层次结构模型；
  2. 构造出各层次中的所有判断矩阵；
  3. 层次单排序及一致性检验；
  4. 层次总排序及一致性检验。

构建层次结构模型

我们把将 决策的目标，考虑的因素，还有决策对象按照它们之间的相互关系分为最高层，中间层和最低层。

  • 最高层：是我们决策的目的和解决的问题，为_目标层_O
  • 最低层：决策的备选方案，为_方案层_P
  • 中间层：考虑的因素和决策的准则，为_准则层_C

每层有若干元素, 各层元素间的关系用相连的直线表示。

mm8-ex-trip

在如何选择旅游地这一事例中，准则层里面还可以加，比如说：安全系数，交通便利等等因素。

接下来我们就要确定各层次各因素的权重，比如说我们认为，这次旅行的安全是最重要的，那么我们可以给他这个权重赋值10，接下来景色也很重要，将它赋值8，最不重要的比如说是距离，将它赋值为1。

但这样感觉，特别的随意，而且常常这个模型不能够推广，也不容易被别人接受，主观性太强。

  • 在确定影响某因素的诸因子在该因素中所占的比重时,遇到的主要困难是这些比重常常不易定量化。
  • 此外,当影响某因素的因子较多时,直接考虑各因子对该因素有多大程度的影响时,常常会因考虑不周全、顾此失彼而使决策者提出与他实际认为的重要性程度不相一致的数据,甚至有可能提出一组隐含矛盾的数据

所以Santy先生提出了一致矩阵法。

构造判断矩阵或对比矩阵

一致矩阵法的核心如下：

  1. 不把所有的因素放在一起比较，只是两两之间相互比较，全部比较结果用矩阵 $A = ( a_{ij} )_{n\times n}$ 表示,称 $A$ 为成对比较判断矩阵(简称判断矩阵)。
  2. 是对比的时候采用相对的尺度，尽量减少因为性质不同而造成的各个因素之间相互比较的困难，以提高准确度。
  3. 一般地作$\frac{n(n-1)}2$次两两判断是必要的。
    • 可以提供更多的信息, 通过各种不同角度的反复比较, 从而导出一个合理的排序。
    • 把所有元素都和某个元素比较的弊病在于, 任何一个判断的失误均可导致不合理的排序, 而个别判断的失误对于难以定量的系统往往是难以避免的。

我们可以参考成对比较矩阵的标度表。

标度 含义
1 两个因素具有相同的重要性
3 一个比另一个稍微重要
5 一个比另一个明显重要
7 一个比另一个强烈重要
9 一个比另一个极端重要
2,4,6,8 相邻判断的中值
  • 从心理学观点来看,分级太多会超越人们的判断能力,既增加了作判断的难度,又容易因此而提供虚假数据。
  • Saaty 等人还用实验方法比较了在各种不同标度下人们判断结果的正确性,实验结果也表明,采用 1~9 标度最为合适。

mm8-trip-pairwise

  • 矩阵元素的值$a_{ij}=C_i:C_j$，是准则两两对比的相对重要性。
  • 在这个矩阵里面C1比C1是1，他自己和自己比肯定都是一样重要的。
  • C1:C2是1:2。那么C2(费用)比C1(景色)稍微重要一些，而且这个重要程度要少于标度表中的 “稍微重要”。

定义 1.
满足$a_{ij}>0$，且$a_{ij}a_{ji}=1$的矩阵$A=(a_{ij})$，称为正互反矩阵。

判断矩阵是正互反矩阵

层次单排序及一致性检验

定理 1.
正互反矩阵 $A$ 的最大特征根 $\lambda_{max}$ 必为正实数,其对应特征向量的所有分量均为正实数。 $A$ 的其余特征值的模均严格小于 $\lambda_{max}$。

判断矩阵 $A$ 对应于最大特征值 $\lambda_{max}$ 的特征向量 $W$ ,经归一化后即为同一层次相应因素对于上一层次某因素相对重要性的排序权值,这一过程称为层次单排序。

定义 2.
满足$a_{ik}a_{kj}=a_{ij}$的正互反矩阵称为一致矩阵。

定理 2.
若 $A$ 为一致矩阵,则

  1. $A$ 必为正互反矩阵。
  2. $A$ 的转置矩阵 $A^T$ 也是一致矩阵。
  3. $A$ 的任意两行成比例,比例因子大于零, 从而 $rank ( A ) = 1$ (同样, $A$ 的任意两列也成比例)。

若 $A$ 为一致矩阵,则

  1. $A$ 的最大特征值 $\lambda_{max} = n$ ,其中 $n$ 为矩阵 $A$ 的阶。 $A$ 的其余特征根均为零。
  2. 若 $A$ 的最大特征值 $\lambda_{max}$ 对应的特征向量为 $W = ( w_1 , \cdots , w_n )^T$ ,则 $a_{ij} = \frac{w_i}{w_j}$
  • 这里的$W$就是想要知道的权重，
  • 所以通过求比较矩阵的最大特征值所对应的特征向量，就可以获得不同因素的权重，
  • 归一化一下（每个权重除以权重和作为自己的值，最终总和为1）就更便于使用了。

构造的成对比较判断矩阵虽能减少其它因素的干扰, 较客观地反映出一对因子影响力的差别。但综合全部比较结果时,其中难免包含一定程度的非一致性。

mm8-trip-pairwise

如： $C2:C1=2$, $C1:C3=4$，则$C2:C3$应该是$8$，但实际为$7$

一致性检验: 对A确定不一致的允许范围

定理 3.
$n$ 阶正互反矩阵 $A$ 为一致矩阵当且仅当其最大特征根 $\lambda_{max}=n$, 且当正互反矩阵 $A$ 非一致时, 必有$\lambda_{max}>n$

  • 根据定理，可以由 $\lambda_{max}$ 是否等于 $n$ 来检验判断矩阵 $A$ 是否为一致矩阵。
  • $\lambda_{max}$ 比 $n$大得越多, $A$ 的非一致性程度也就越严重，$\lambda_{max}$对应的标准化特征向量也就越不能真实地反映出因素影响的比重

对判断矩阵的一致性检验的步骤如下:

  1. 计算$CI=\frac{\lambda_{max}-n}{n-1}$
  2. 查找相应的平均随机一致性指标(Random consistency Index) RI。对$n$由1到15，有RI为 
    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]
  3. 计算一致性比例 $CR=\frac{CI}{RI}$。当$CR<0.10$时，认为判断矩阵的一致性是可以接受的,否则应对判断矩阵作适当修正。

“选择旅游地”中准则层对目标的权向量及一致性检验

  1. 最大特征根$\lambda=5.073$
  2. 一致性指标: $CI=\frac{\lambda-5}{5-4}=0.018021102142554035$
  3. 随机一致性指标 $RI=1.12$ (查表)
  4. 一致性比率$CR=\frac{CI}{RI}=\frac{0.018}{1.12}=0.016$，$<0.1$，通过一致性检验
  5. 相应特征向量$w^{(2)}$ 
    相应特征向量:  [0.46582183+0.j 0.84086331+0.j 0.09509743+0.j 0.17329948+0.j
0.19204866+0.j]
import numpy as np
# 判断矩阵
aa=np.array([[1.0,1.0/2.0, 4.0, 3.0, 3.0],
             [2.0,1.0,7.0,5.0,5.0],
             [1/4.0, 1/7.0, 1, 1/2.0, 1/3.0],
             [1/3.0, 1/5.0, 2.0, 1.0, 1.0],
             [1/3.0, 1/5.0, 3.0, 1.0, 1.0]
            ])
na=aa.shape[0] # 矩阵大小
eigs, eigvecs=np.linalg.eig(aa) # 特征根
rho=np.max(np.abs(eigs)) # 最大特征值
rhoind=np.argmax(np.abs(eigs))

# 一致性指标
ci = (rho-na)/(na-1)
# 随机性指标 (Random consistancy Index)
ri=np.array([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])
# 一致性比率， <0.1 表示通过一致性检验
cr=ci/ri[na-1]
if cr<0.1:
  print("矩阵大小: ",na)
  print("最大特征值: ",rho)
  print("相应特征向量: ",eigvecs[:,rhoind])
  print("一致性指标: ",ci)
  print("通过一致性检验(cr<0.1)", cr)
else:
  print("--未通过--一致性检验(cr>0.1)", cr)

层次总排序及一致性检验

同样求第3层(方案)对第2层每一元素(准则)的权向量

mm8-ex-trip-criteria-pairwise

不同的矩阵之间确定的这些系数，我们通过查阅资料和自己主观判断来填写完成。

同样计算每一个矩阵的最大特征值，相应的特征向量，CR值等

k 1 2 3 4 5
$w_k^{(3)}$ [0.89021421 0.41320083 0.19179084] [0.11283181 0.32546327 0.93879851] [0.6882472 0.6882472 0.22941573] [0.92551596 0.28029569 0.25466553] [0.23570226 0.23570226 0.94280904]
$\lambda_k$ 3.005 3.002 2.9999 3.009 3
$CI_k$ 0.003 0.001 0 0.005 0
$RI_k$ 0.58 0.58 0.58 0.58 0.58
$CR_k$ 0.005 0.001 -8e-16 0.009 0

记矩阵$B=(w_1^{(3)},\cdots,w_k^{(3)})$，则最终方案层对目标的组合权向量(层次总排序权重)为

\[v=B\cdot w^{(2)} \]
[0.78066562+0.j 0.62544003+0.j 1.12575705+0.j]

对层次总排序也需作一致性检验,检验仍象层次总排序那样由高层到低层逐层进 行。这是因为虽然各层次均已经过层次单排序的一致性检验,各成对比较判断矩阵都 已具有较为满意的一致性。但当综合考察时,各层次的非一致性仍有可能积累起来, 引起最终分析结果较严重的非一致性。

记$w^{(2)}=(a_1,\cdots,a_m)$，则总排序随机一致性比例为

\[CR=\frac{\sum_{k=1}^m CI_k a_k}{\sum_{k=1}^m RI_ka_k} \]
总排序随机一致性比例:  (0.0029760853941554807+0j)

# 方案层对准则层的判断矩阵
b1=np.array([[1.0, 2, 5],[1/2.0, 1, 2.0],[1/5.0, 1/2.0, 1]])
b2=np.array([[1, 1/3.0, 1/8.0],[3, 1, 1/3.0],[8,3,1]])
b3=np.array([[1,1,3],[1,1,3],[1/3.0,1/3.0,1]])
b4=np.array([[1,3,4],[1/3.0,1,1],[1/4.0,1,1]])
b5=np.array([[1,1,1/4.0],[1,1,1/4.0],[4,4,1]])


bb=np.array([b1,b2,b3,b4,b5])
#print(bb[0])
#print(b1)
evals_b=[]
evecs_b=[]

for bi in bb:
    nb=bi.shape[0] # size
    eigs, eigvecs=np.linalg.eig(bi) # 特征根
    rho=np.max(np.abs(eigs)) # 最大特征值
    rhoind=np.argmax(np.abs(eigs))
    evec=eigvecs[:,rhoind]  # 对应的特征向量
    if evec[0]<0.0:
        evec=-evec
    #print(np.sum(np.power(wa,2)))
    #evals_b.append(rho)
    #print(evec)
    evecs_b.append(evec.copy())
    # 一致性指标
    ci = (rho-nb)/(nb-1)
    # 一致性比率， <0.1 表示通过一致性检验
    cr=ci/ri[nb-1]
    evals_b.append((rho,ci,cr,ri[nb-1]))

    #print(bi)

#print(evals_b)
ci_b=np.array(evals_b)
mat_b=np.array(evecs_b).T
print(mat_b)  # 各个特征向量
print(ci_b)  # 各个特征值，CI，CR，RI

final=np.dot(mat_b, wa)
print("组合权向量: ",final)

# 总排序随机一致性比例
cr_f=np.dot(ci_b[:,1],wa)/np.dot(ci_b[:,3],wa)
print("总排序随机一致性比例: ",cr_f)

应用

  • 应用领域: 经济计划和管理, 能源政策和分配, 人才选拔和评价, 生产决策, 交通运输, 科研选题, 产业结构, 教育, 医疗, 环境, 军事等。
  • 处理问题类型: 决策、评价、分析、预测等。
  • 建立层次分析结构模型是关键一步,要有主要决策层参与。
  • 构造成对比较阵是数量依据, 应由经验丰富、判断力强的专家给出。

例 2. 国家综合实力分析

mm8-ex-state

例 3. 挑选合适的工作。经双方恳谈,已有三个单位表示愿意录用某毕业生。该生根据已有信息建立了一个层次结构模型

mm8-ex-job

mm8-ex-job-c mm8-ex-job-p

mm8-ex-job-f

例 4. 渡河方案

mm8-ex-river

mm8-ex-river-dis

例 5. 科技成果的综合评价

mm8-ex-tech

不完全结构分析

例 6. 评价教师贡献的层次结构

mm8-ex-teacher

$P_1$, $P_2$只作教学，$P_4$只作科研，$P_3$兼作教学、科研

上一层第一元素与下层所有元素相关联的结构称为完全层次结构。否则为不完全层次结构

  • 如果元素不相关，则权向量的相应位置设置为0。

    设第2层对第1层权向量$w^{(2)}=(w^{(2)}_1,w^{(2)}_2)^T$已定，第3层对第2层权向量 $w_1^{(3)}=(w^{(3)}_{11},w^{(3)}_{12},w^{(3)}_{13},0)^T$ $w_2^{(3)}=(0,0,w^{(3)}_{23},w^{(3)}_{24})^T$已得

  • 如何得到第3层对第1层权向量？

    特例： $C_1$与$C_2$同等重要，则$w^{(2)}=(\frac12,\frac12)^T$。$P_1$,$\cdots$,$P_4$的能力相同，则$w^{(3)}_1=(\frac13,\frac13,\frac13)^T$, $w^{(3)}_2=(0,0,\frac12,\frac12)^T$。则公正的评价应该是：$P_1:P_2:P_3:P_4=1:1:2:1$
  • 不考虑支配元素数目的影响
    \[w^{(3)}=(w^{(3)}_1,w^{(3)}_2)w^{(2)}=(\frac16,\frac16,\frac5{12},\frac14)^T \]
  • 支配元素越多权重越大。用支配元素数目$n_1$,$n_2$对$w^{(2)}$加权修正
    \[\tilde w^{(2)}=\frac{1}{n_1w^{(2)}_1+n_2w^{(2)}_2}(n_1w^{(2)}_1, n_2w^{(2)}_2)^T \]
    然后用$\tilde w^{(2)}$计算最后的权重
    \[w^{(3)}=(w^{(3)}_1,w^{(3)}_2)\tilde w^{(2)}=(\frac15,\frac15,\frac2{5},\frac15)^T \]
  • 支配元素越多权重越小

TOPSIS

例 7. 为了客观地评价我国研究生教育的实际状况和各研究生院的教学质量，国务院学位委员会办公室组织过一次研究生院的评估。为了取得经验，先选5所研究生院，收集有关数据资料进行了试评估，下表是所给出的部分数据：

mm8-ex-graduated

  • TOPSIS 法是一种常用的组内综合评价方法

C.L.Hwang 和 K.Yoon 于1981年首次提出 TOPSIS (Technique for Order Preference by Similarity to an Ideal Solution)。

  • TOPSIS 法是一种常用的组内综合评价方法，能充分利用原始数据的信息，其结果能精确地反映各评价方案之间的差距。
  • 基本过程为基于归一化后的原始数据矩阵，采用余弦法找出有限方案中的最优方案和最劣方案，然后分别计算各评价对象与最优方案和最劣方案间的距离，获得各评价对象与最优方案的相对接近程度，以此作为评价优劣的依据。
  • 该方法对数据分布及样本含量没有严格限制，数据计算简单易行。

通俗的例子：小明数学考试 134 分，要怎么知道他的成绩是好还是不好呢？

  • 基于分布的评价方法会观察小明的分数位于班级分数的哪个水平（如前 5%、前 10%），但这种评价方法只能给出一个方向的情况。如班上成绩除了最高分外，其余都是 134 分，那么小明的成绩就是并列的倒数第一，但是正向评价给出的结果是前 5%。
  • 而 TOPSIS 就是找出班上最高分（假设是 147 分）、最低分（假设是 69 分），然后计算小明的分数和这两个分数之间的差距，从而得到自己分数好坏的一个客观评价。距离最高分越近，那么评价情况越好，距离最低分越近，那么评价情况越糟。

属性值的规范化

数据的属性值具有多种类型，包括效益型、成本型、区间型。

  • 效益型： 属性值越大越好。比如人均专著，当然是越多越好。
  • 成本型： 属性值越小越好。比如逾期毕业率，越低越好。
  • 区间型： 属性值在某个区间最佳。比如师生比最优区间为$[5,6]$，过大或者过小都不好。

属性值的规范化包括：

  1. 对各种类型的数据，变为性能越优的方案，变换后属性值越大，即指标属性同向化
  2. 排除量纲的选用对决策或评估结果的影响，即无量纲化
  3. 把属性值表中的数据归一化，即数据变换到$[0,1]$区间上

常用的属性规范化方法

  1. 线性变换。变换后的最差属性值不一定为$0$，但最优值为1
    • 若$x$为效益型，则$x'=\frac{x}{x_{max}}$，其中$x_{max}$为最大值
    • 若$x$为成本型，则$x'=1-\frac{x}{x_{max}}$
  2. 标准0-1变换。$x_{min}$为最小值，
    • 若$x$为效益型
      \[x'=\frac{x-x_{min}}{x_{max}-x_{min}} \]
    • 若$x$为成长型
      \[x'=\frac{x_{max}-x}{x_{max}-x_{min}} \]
  1. 成本型的变换。若$x>0$，则取$x'=\frac1x$。或者取$x'=x_{max}-x$。
  2. 区间属性的变换。先设置最优区间$[a,b]$，然后给出数据的容忍区间$[a',b']$，取
    \[x'=\begin{cases} 1-\frac{a-x}{a-a'} , & a'\leq x\leq a \\ 1, & a\leq x\leq b \\ 1-\frac{x-b}{b'-b}, & b\leq x\leq b' \\ 0, & otherwise \end{cases} \]
  1. 向量规范化。无论是成本型，还是效益型，使用相同的公式
    \[x'_i=\frac{x_i}{\sqrt{\sum_{k=1}^nx_k^2}} \]
  2. 标准化处理。可以消除变量的量纲化效应，
    \[x'_i=\frac{x_i-\bar x}{s} \]
    其中
    \[\bar x=\frac1n\sum_{k=1}^nx_k \ , \ \ \ s=\sqrt{\frac1{n-1}\sum_{k=1}^n(x_k-\bar x)^2} \]

python scipy.stats.zscore 可以实现标准化处理数据

对于原始决策矩阵(原始数据)，先做同向化处理（性能越优的，值越大），然后做规范化

  • 同向化。对于“生师比”，取最优区间为$[5,6]$，容忍区间为$[2,12]$。对于“逾期毕业率”，取$\frac1x$。 得到 
    === 数据同向化
   人均专著       生师比   科研经费     逾期毕业率
院校A   0.1  1.000000   5000  0.212766
院校B   0.2  1.000000   6000  0.178571
院校C   0.4  0.833333   7000  0.149254
院校D   0.9  0.333333  10000  0.434783
院校E   1.2  0.000000    400  0.555556
  • 规范化。 
           人均专著       生师比      科研经费     逾期毕业率
院校A  0.063758  0.597022  0.344901  0.275343
院校B  0.127515  0.597022  0.413882  0.231092
院校C  0.255031  0.497519  0.482862  0.193151
院校D  0.573819  0.199007  0.689803  0.562658
院校E  0.765092  0.000000  0.027592  0.718952
# 数据
data = pd.DataFrame(
    {'人均专著': [0.1, 0.2, 0.4, 0.9, 1.2], '生师比': [5, 6, 7, 10, 2], 
     '科研经费': [5000, 6000, 7000, 10000, 400],
     '逾期毕业率': [4.7, 5.6, 6.7, 2.3, 1.8]}, index=['院校' + i for i in list('ABCDE')])
# 数据同向化
data['生师比'] = dataDirection_3(data['生师比'], 5, 6, 2, 12)   # 师生比数据为区间型指标
data['逾期毕业率'] = dataDirection_1(data['逾期毕业率']) # 逾期毕业率为成本型指标
print("=== 数据同向化")
print(data)
# 数据规范化
for column in data.columns:
    data[column]=data[column]/np.sqrt(np.sum(np.power(data[column],2)))
    print(column)
print(data)

构造加权规范阵

设由决策人给定各属性的权重向量为$\vec w=(w_1,\cdots,w_n)^T$，则加权规范矩阵$C=(c_{ij})$为

\[c_ij=w_j\cdot b_{ij} \]

其中$B=(b_{ij})$是前面得到的规范决策矩阵

取权为$w=(0.2,0.3,0.4,0.1)^T$，则 

== 加权规范决策矩阵
         人均专著       生师比      科研经费     逾期毕业率
院校A  0.012752  0.179107  0.137961  0.027534
院校B  0.025503  0.179107  0.165553  0.023109
院校C  0.051006  0.149256  0.193145  0.019315
院校D  0.114764  0.059702  0.275921  0.056266
院校E  0.153018  0.000000  0.011037  0.071895
# 加权规范矩阵
w=np.array([0.2,0.3,0.4,0.1])
i=0
for column in data.columns:
    data[column]=data[column]*w[i]
    i=i+1
print("== 加权规范决策矩阵")
print(data)

确定正负理想解

在每一列中分别选取最大值、最小值构成最优（正理想解$C^*$）最劣（负理想解$C^0$）方案

得到结果 

== 正理想解（最优方案）
人均专著     0.153018
生师比      0.179107
科研经费     0.275921
逾期毕业率    0.071895
dtype: float64
== 负理想解（最劣方案）
人均专著     0.012752
生师比      0.000000
科研经费     0.011037
逾期毕业率    0.019315

综合评价指数

使用Eulerian距离。方案$d_i$到正理想解的距离

\[s^*_i=\sqrt{\sum_{k=1}^n(c_{ij}-c^*_j)^2} \]

到负理想解的距离

\[s^0_i=\sqrt{\sum_{k=1}^n(c_{ij}-c^0_j)^2} \]

综合评价指数由上面的距离计算

\[f_i=\frac{s_i^0}{s_i^0+s_i^*} \]

计算出的距离 

== 到正、负理想的距离
         人均专著       生师比      科研经费     逾期毕业率       正距离       负距离
院校A  0.012752  0.179107  0.137961  0.027534  0.201682  0.219673
院校B  0.025503  0.179107  0.165553  0.023109  0.175560  0.236921
院校C  0.051006  0.149256  0.193145  0.019315  0.144617  0.238545
院校D  0.114764  0.059702  0.275921  0.056266  0.126353  0.292404
院校E  0.153018  0.000000  0.011037  0.071895  0.319754  0.149798

最终的排序 

== 综合评价
     人均专著  生师比   科研经费  逾期毕业率    综合评价指数  评价排序
院校A   0.1    5   5000    4.7  0.478651   4.0
院校B   0.2    6   6000    5.6  0.425621   3.0
院校C   0.4    7   7000    6.7  0.377431   2.0
院校D   0.9   10  10000    2.3  0.301734   1.0
院校E   1.2    2    400    1.8  0.680977   5.0
# 到正、负理想解的距离
norm_max=[]
norm_min=[]
#f_star=[] # 综合评价指数
for index,row in data.iterrows():
    norm_max.append(np.linalg.norm(row-c_max))
    norm_min.append(np.linalg.norm(row-c_min))
    #print(index,np.linalg.norm(row-c_max))
data['正距离']=np.array(norm_max)
data['负距离']=np.array(norm_min)
print("== 到正、负理想的距离")
print(data)

data0=origin_data.copy()
data0['综合评价指数']=data['正距离']/(data['正距离']+data['负距离'])
data0['评价排序']=data0['综合评价指数'].rank(ascending=1,method='first')
print("== 综合评价")
print(data0)
  • 分析中使用的各个指标的权重，对TOPSIS的影响很大。
  • 权重除了主观方法外，也可以用熵值法。方法根据各指标所含信息有序程度的差异性来客观确定指标的权重，仅依赖于数据本身的离散程度。
  • 权重还可以用AHP（层次分析）方法来得到。

RSR秩和比综合评价

RSR法(Rank-sum ratio) 是我国田凤调教授于1988年提出的，集古典参数统计与近代非参数统计各自优点于一体的统计分析方法。RSR法现在广泛地应用于医疗卫生、科技、经济等邻域的多指标综合评价、统计预测预报、鉴别分类、统计质量控制等方面。

一般过程是将效益型指标从小到大排序进行排名、成本型指标从大到小排序进行排名，再计算秩和比，最后统计回归、分档排序。通过秩转换，获得无量纲统计量RSR；在此基础上，运用参数统计分析的概念与方法，研究RSR的分布；以RSR值对评价对象的优劣直接排序或分档排序，从而对评价对象做出综合评价。

优点：以非参数法为基础，对指标的选择无特殊要求，适用于各种评价对象，由于计算时使用的数值是秩次，可以消除异常值的干扰

缺点：排序的主要依据是利用原始数据的秩次，最终算得的 RSR 值反映的是综合秩次的差距，而与原始数据的顺位间的差距程度大小无关，这样在指标转化为秩次是会失去一些原始数据的信息，如原始数据的大小差别等。当 RSR 值实际上不满足正态分布时，分档归类的结果与实际情况会有偏差，且只能回答分级程度是否有差别，不能进一步回答具体的差别情况。

例 8. 请对某省10个地区孕产妇保健工作进行综合评价

地区编码产前检查率(%)孕妇死亡率(1/100000)围产儿死亡率(1/1000)
A 99.54 60.27 16.15
B 96.52 59.6720.10
C 99.3643.9115.60
D92.8358.9917.04
E 91.7135.4015.01
F95.3544.7113.93
G96.0949.8117.43
H99.2731.6913.89
I94.7622.9119.87
J84.8081.4923.63

编秩

RSR方法的第一步，列出原始数据表并编秩(Rank)。得到秩矩阵，记为 $R=(r_{ij})_{m\times n}$。

  • 整数次秩。将每个指标各评价对象按大小排序，序号就是秩。效益型指标从小到大编秩，成本型指标从大到小编秩，同一指标数据相同者编平均秩。
  • 非整数次秩。用类似线性插值的方式对指标值时行编秩
    • 对于效益型指标：
      \[r_j=1+(n-1)\frac{x_j-\min(x_1,\cdots,x_n)}{\max(x_1,\cdots,x_n)-\min(x_1,\cdots,x_n)} \]
    • 对于成本型指标：
      \[r_j=1+(n-1)\frac{\max(x_1,\cdots,x_n)-x_j}{\max(x_1,\cdots,x_n)-\min(x_1,\cdots,x_n)} \]
   X1：产前检查率  R1：产前检查率  X2：孕妇死亡率  R2：孕妇死亡率  X3：围产儿死亡率  R3：围产儿死亡率
A     99.54      10.0     60.27       9.0      16.15        5.0
B     96.52       7.0     59.67       8.0      20.10        9.0
C     99.36       9.0     43.91       4.0      15.60        4.0
D     92.83       3.0     58.99       7.0      17.04        6.0
E     91.71       2.0     35.40       3.0      15.01        3.0
F     95.35       5.0     44.71       5.0      13.93        2.0
G     96.09       6.0     49.81       6.0      17.43        7.0
H     99.27       8.0     31.69       2.0      13.89        1.0
I     94.76       4.0     22.91       1.0      19.87        8.0
J     84.80       1.0     81.49      10.0      23.63       10.0
data = pd.DataFrame({'产前检查率': [99.54, 96.52, 99.36, 92.83, 91.71, 95.35, 96.09, 99.27, 94.76, 84.80],
                     '孕妇死亡率': [60.27, 59.67, 43.91, 58.99, 35.40, 44.71, 49.81, 31.69, 22.91, 81.49],
                     '围产儿死亡率': [16.15, 20.10, 15.60, 17.04, 15.01, 13.93, 17.43, 13.89, 19.87, 23.63]},
                    index=list('ABCDEFGHIJ'), columns=['产前检查率', '孕妇死亡率', '围产儿死亡率'])

# 处理成本型数据
data['孕妇死亡率']=-data['孕妇死亡率']
data['围产儿死亡率']=-data['围产儿死亡率']

# 对原始数据编秩
for i, X in enumerate(data.columns):
            #print(i,X,type(X),type(i))
            aa="X{index}:{name}".format(index=i,name=X)
            Result["X{index}：{name}".format(index=i+1, name=X)] = data.iloc[:, i]
            Result['R{index}：{name}'.format(index=i+1, name=X)] = data.iloc[:, i].rank(method="dense")

计算秩和比并排序

当各个指标的权重不一样时，计算秩和比(RSR)

\[WSRS_i=\frac1n\sum_{j=1}^nw_jr_{ij} , i=1,2,\cdots,n \]

若各个指标权重相同时，则$w_j=\frac1m$，则秩和比为

\[SRS_i=\frac1{mn}\sum_{j=1}^nr_{ij} , i=1,2,\cdots,n \]
==== 秩和比(RSR):
   X1：产前检查率  R1：产前检查率  X2：孕妇死亡率  ...  R3：围产儿死亡率       RSR  RSR_Rank
A     99.54      10.0    -60.27  ...        6.0  0.600000       4.5
B     96.52       7.0    -59.67  ...        2.0  0.400000       8.5
C     99.36       9.0    -43.91  ...        7.0  0.766667       2.0
D     92.83       3.0    -58.99  ...        5.0  0.400000       8.5
E     91.71       2.0    -35.40  ...        8.0  0.600000       4.5
F     95.35       5.0    -44.71  ...        9.0  0.666667       3.0
G     96.09       6.0    -49.81  ...        4.0  0.500000       7.0
H     99.27       8.0    -31.69  ...       10.0  0.900000       1.0
I     94.76       4.0    -22.91  ...        3.0  0.566667       6.0
J     84.80       1.0    -81.49  ...        1.0  0.100000      10.0
    # 计算秩和比
    weight = 1 / m if weight is None else pd.np.array(weight) / sum(weight)
    Result['RSR'] = (Result.iloc[:, 1::2] * weight).sum(axis=1) / n
    Result['RSR_Rank'] = Result['RSR'].rank(ascending=False)
    print("==== 秩和比(RSR):")
    print(Result)

计算概率单位

按从小到大的顺序编制RSR（或WRSR）频率分布表，

  1. 列出各组频数$f_i$，
  2. 计算各组累积频数$cf_i$，
  3. 计算累积频率$p_i=\frac{cf_i}n$，最后一项记为 $1-\frac1{4n}$ 进行修正。
  4. 将$p_i$转换为$Probit_i$ ($Probit_i$为标准正态分布的$p_i$分位数加5) 
           f  Σ f  \bar{R} f  \bar{R}/n*100%    Probit
0.100000  1    1        1.0           0.100  3.718448
0.400000  2    3        2.5           0.250  4.325510
0.500000  1    4        4.0           0.400  4.746653
0.566667  1    5        5.0           0.500  5.000000
0.600000  2    7        6.5           0.650  5.385320
0.666667  1    8        8.0           0.800  5.841621
0.766667  1    9        9.0           0.900  6.281552
0.900000  1   10       10.0           0.975  6.959964
    # 绘制 RSR 分布表
    RSR = Result['RSR']
    RSR_RANK_DICT = dict(zip(RSR.values, RSR.rank().values))
    Distribution = pd.DataFrame(index=sorted(RSR.unique()))
    Distribution['f'] = RSR.value_counts().sort_index()
    Distribution['Σ f'] = Distribution['f'].cumsum()
    Distribution[r'\bar{R} f'] = [RSR_RANK_DICT[i] for i in Distribution.index]
    Distribution[r'\bar{R}/n*100%'] = Distribution[r'\bar{R} f'] / n
    Distribution.iat[-1, -1] = 1 - 1 / (4 * n)
    Distribution['Probit'] = 5 - norm.isf(Distribution.iloc[:, -1])

计算直线回归方程

以累积频率所对应的概率单位 $Probit_i$ 为自变量，以 $RSR_i$(或$WRSR_i$) 值为因变量，计算直线回归方程，即：

\[RSR=a+b\times Probit \]

回归方程检验：对该回归方程，需要进行检验。回归方程的检验往往围绕以下几点展开：

  • 残差独立性检验：Durbin-Watson 检验
  • 方差齐性检验（异方差性）：Breusch-Pagan 检验和 White 检验
  • 回归系数 $b$ 的有效性检验： $t$ 检验法和置信区间检验法
  • 拟合优度检验：(自由度校正)决定系数、Pearson相关系数、Spearman秩相关系数、交叉验证法等 
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.938
Model:                            OLS   Adj. R-squared:                  0.928
Method:                 Least Squares   F-statistic:                     90.63
Date:                Thu, 16 May 2019   Prob (F-statistic):           7.66e-05
Time:                        19:40:42   Log-Likelihood:                 11.629
No. Observations:                   8   AIC:                            -19.26
Df Residuals:                       6   BIC:                            -19.10
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const         -0.6085      0.125     -4.862      0.003      -0.915      -0.302
Probit         0.2217      0.023      9.520      0.000       0.165       0.279
==============================================================================
Omnibus:                        1.731   Durbin-Watson:                   1.240
Prob(Omnibus):                  0.421   Jarque-Bera (JB):                0.682
Skew:                          -0.695   Prob(JB):                        0.711
Kurtosis:                       2.657   Cond. No.                         30.1
==============================================================================
import statsmodels.api as sm
# 计算回归方差并进行回归分析
r0 = pd.np.polyfit(Distribution['Probit'], Distribution.index, deg=1)
print(sm.OLS(Distribution.index, sm.add_constant(Distribution['Probit'])).fit().summary())
if r0[1] > 0:
        print("\n回归直线方程为：y = {a} Probit + {b}".format(a=r0[0],b=r0[1]))
else:
        print("\n回归直线方程为：y = {a} Probit - {b}".format(a=r0[0],b=abs(r0[1])))

对RSR矫正值进行分档排序

按照回归方程推算所对应的RSR估计值对评价对象进行分档排序，分档数由研究者根据实际情况决定。 

   X1：产前检查率  R1：产前检查率  X2：孕妇死亡率  ...    Probit  RSR Regression  Level
A     99.54      10.0    -60.27  ...  5.385320        0.585320      2
B     96.52       7.0    -59.67  ...  4.325510        0.350371      2
C     99.36       9.0    -43.91  ...  6.281552        0.784005      1
D     92.83       3.0    -58.99  ...  4.325510        0.350371      2
E     91.71       2.0    -35.40  ...  5.385320        0.585320      2
F     95.35       5.0    -44.71  ...  5.841621        0.686477      2
G     96.09       6.0    -49.81  ...  4.746653        0.443734      2
H     99.27       8.0    -31.69  ...  6.959964        0.934402      1
I     94.76       4.0    -22.91  ...  5.000000        0.499898      2
J     84.80       1.0    -81.49  ...  3.718448        0.215792      3
    # 代入回归方程并分档排序
    Result['Probit'] = Result['RSR'].apply(lambda item: Distribution.at[item, 'Probit'])
    Result['RSR Regression'] = pd.np.polyval(r0, Result['Probit'])
    threshold = pd.np.polyval(r0, [2, 4, 6, 8]) if threshold is None else pd.np.polyval(r0, threshold)
    Result['Level'] = pd.cut(Result['RSR Regression'], threshold, labels=range(len(threshold) - 1, 0, -1))

目录

https://zhuanlan.zhihu.com/p/37738503 TOPSIS法(优劣解距离法)介绍及 python3 实现

https://zhuanlan.zhihu.com/p/37780755 如何在Word中写出一篇优美的论文？

https://zhuanlan.zhihu.com/p/38209882 RSR(秩和比综合评价法)介绍及python3实现