第八章 评价模型

层次分析法 AHP

在上一篇文章的决策模型中，我们简单讲解了一下不确定性决策、风险决策和多属性决策。在多属性中我们使用了一道购买飞机的例题作为讲解。现在我们还继续用这道例题进行讲解，只不过我们的侧重点有所改变。我们的侧重点改成了如何求解不同决定因素的权值？求解这类问题就要用到我们的层次分析法。

层次分析法是一种决策方法，它可以帮助你在面对复杂的问题时，将问题分解为不同的层次和因素，并根据它们的相对重要性进行比较和评价。层次分析法的模型结构如下：

graph TD
	subgraph 目标层
	H
	end
	subgraph 准则层
	C1
	C2
	C3
	end
	subgraph 方案层
	A1
	A2
	A3
	end
	目标层 --- 准则层
	准则层 --- 方案层

• 最高层：只有一个元素，为预定目标或理想结果；
• 准则层：目标所涉及的中心环节，可以有多个层次。
• 方案层：实现目标的各种措施和解决方案。

构建好模型后，使用判断矩阵分门别类的标注各个准则，每个准则的标号见下表：

建立好判断矩阵后进行一致性检验，最后使用权重做出评价。

例题 购买飞机（续）

某航空公司在国际市场上购买飞机，按6个决策指标对不同型号的飞机进行综合评价，这6个指标是：最大速度、最大范围、最大负载、价格、可靠性、灵敏度。现在4种型号的飞机可供选择，具体指标值见表。

现要求分析不同影响因素的权重。

首先我们确定模型的结构图：

graph TD
  H[选择最优的飞机]
	C1[最大速度]
	C2[最大范围]
	C3[费用]
	C4[最大负载]
	C5[可靠性]
	C6[灵敏度]
	
	subgraph 目标层
	H
	end
	subgraph 准则层
	C1
	C2
	C3
	C4
	C5
	C6
	end
	subgraph 方案层
	A1[1号飞机]
	A2[2号飞机]
	A3[3号飞机]
	A4[4号飞机]
	end
	目标层 --- 准则层
	准则层 --- 方案层

对于准则层，我们需要根据主观来对每个准则进行评价。构建评价矩阵，第n行第m列代表n比m指标的重要程度。构建好的矩阵如下：

A=\begin{bmatrix} 1&3&5&1/5&1/3&3\\ 1/3&1&3&1/7&1/5&1/3\\ 1/5&1/3&1&1/9&1/7&1/5\\ 5&7&9&1&3&5\\ 3&5&7&1/3&1&3\\ 1/3&3&5&1/5&1/3&1\\ \end{bmatrix}

通过分析这个判断矩阵可以发现，行列翻转的值对应互为倒数。接下来我们需要对这个矩阵进行一致性检验。

首先定义一致性指标：

CI=\frac{\lambda - n}{n-1}

其中：\lambda为矩阵的特征值，n为唯一非零特征根。

这时候我们再引入随机一致性指标RI，这个指标是给定的，如下表：

然后我们确定一致性比率CR：

CR=\frac{CI}{RI}

一般来说，当我们的一致性比率CR<0.1的时候，在认为A不一致程度在允许的范围内，反之则没有通过一致性检验。我们下面来计算上面构建的矩阵的一致性比率。

# 计算矩阵的一致性比率
import numpy as np

def caculate_consistency_ratio(m: np.array):
    """
    计算矩阵的一致性比率，方阵的n小于等于9。
    :param m: 方阵
    :return: 一致性比率
    """
    # 计算最大特征值
    eigenvalue = np.linalg.eig(m)[0].max()
    # 计算一致性指标
    consistency_index = (eigenvalue - m.shape[0]) / (m.shape[0] - 1)
    # 计算随机一致性指标
    random_consistency_index = np.array(
        [0, 0, 0.58, 0.90, 1.12, 1.24, 1.32, 1.41, 1.45])
    # 计算一致性比率
    consistency_ratio = consistency_index / \
        random_consistency_index[m.shape[0] - 1]
    # 返回一致性比率，为模值
    return abs(consistency_ratio)

if __name__ == '__main__':
    m = np.array([[1, 3, 5, 1/5, 1/3, 3],
                  [1/3, 1, 3, 1/7, 1/5, 1/3],
                  [1/5, 1/3, 1, 1/9, 1/7, 1/5],
                  [5, 7, 9, 1, 3, 5],
                  [3, 5, 7, 1/3, 1, 3],
                  [1/3, 3, 5, 1/5, 1/3, 1]])
    print(caculate_consistency_ratio(m))

计算得到的一致性比率为0.0641，通过一致性检验。

下面我们来看一下一致性检验不通过会发生什么？还是下面的表格:

我们交换了费用和最大负载的数值，这个表格意味着：费用相对于最大速度重要，最大速度相对于最大负载重要，然而费用却没有最大负载重要，这不符合逻辑，我们来重新计算得到新判断矩阵的一致性比率：

计算得CR=0.9894，超过了0.1，故没有通过一致性检验。

之后我们对构造的判断矩阵进行归一化处理，使用算术平均法求权重：

w_{i}=\frac{1}{n}\sum_{j=1}^{n}\frac{a_{ij}}{\sum\limits_{k=1}^{n}a_{kj}}(i=1,2,3,\ldots,n)

# 计算权重
import numpy as np

def weighting(m: np.array):
    """求权重，算术平均法

    Parameters
    ----------
    m : np.array
        一致性矩阵

    Returns
    -------
    np.array
        权重向量
    """
    # 列向量归一化
    m = m / np.sum(m, axis=0)
    # 行向量相加
    row_sum = np.sum(m, axis=1)
    # 求平均
    return row_sum / len(row_sum)

if __name__ == '__main__':
    m = np.array([[1, 3, 5, 1/5, 1/3, 3],
                  [1/3, 1, 3, 1/7, 1/5, 1/3],
                  [1/5, 1/3, 1, 1/9, 1/7, 1/5],
                  [5, 7, 9, 1, 3, 5],
                  [3, 5, 7, 1/3, 1, 3],
                  [1/3, 3, 5, 1/5, 1/3, 1]])
    print(weighting(m))

求得权重矩阵：

w=[0.1383,0.0540,0.0285 ,0.4450,0.2338,0.1004]

对应各个指标的重要性。我们还可以用特征值法，即将特征向量进行归一化：

# 计算权重
import numpy as np

def weighting(m: np.array):
    """求权重，特征值法。
    param: m: np.array, 一个方阵
    return: np.array, 权重
    """
    # 求特征值和特征向量
    eig_val, eig_vec = np.linalg.eig(m)
    # 找到最大特征值的索引
    max_index = np.argmax(eig_val)
    # 取出最大特征值对应的特征向量
    max_eig_vec = eig_vec[:, max_index]
    # 归一化
    max_eig_vec = max_eig_vec / max_eig_vec.sum()
    # 返回模值
    return abs(max_eig_vec)

if __name__ == '__main__':
    m = np.array([[1, 3, 5, 1/5, 1/3, 3],
                  [1/3, 1, 3, 1/7, 1/5, 1/3],
                  [1/5, 1/3, 1, 1/9, 1/7, 1/5],
                  [5, 7, 9, 1, 3, 5],
                  [3, 5, 7, 1/3, 1, 3],
                  [1/3, 3, 5, 1/5, 1/3, 1]])
    print(weighting(m))

w=[0.1365,0.0502,0.0273 ,0.4553 ,0.2376 ,0.0931]

和上一种计算方法算出的结果相差不大。

TOPSIS法

我们在概率模型中用一道例题解释了TOPSIS法的基本步骤，在这里我们会用另一道例题练习TOPSIS法。

TOPSIS法（Technique for Order Preference by Similarity to Ideal Solution）可翻译为逼近理想解排序法，国内常简称为优劣解距离法

TOPSIS法是一种常用的综合评价方法，其能充分利用原始数据的信息，其结果能精确地反映各评价方案之间的差距，此外层次分析法对指标或方案过多的情况处理较为麻烦。

层次分析法的缺点是我们需要主观的对每个指标进行评价，使用TOPSIS法做决策，我们可以不用计算每个指标的权重从而进行决策。

在使用TOPSIS法之前，我们首先需要对数据进行预处理。预处理指标分为四类，分别是：

每个指标的不同，对应的数据处理方式不同，我们需要将所有的指标转换成极大型指标：

1. 极小型指标的正向化处理

   对于极小值类型：

   \hat{x}_i=\max(x)-x_i

2. 对于中间性指标的正向化处理

   最佳值为x_{best}，取M=\max\left\{|x_{i}-x_{best}|\right\}，再使用公式：

   \hat{x}_i=1-\frac{x_i-x_{best}}M

   来进行转化。

3. 对于区间型指标，如果最佳区间是[a,b]，取M=\max\left\{a-\min\left\{x_{i}\right\},\max\left\{x_{i}\right\}-b\right\}，然后使用以下公式：

   \hat{x}_i=\begin{cases}1-\dfrac{a-x_i}{M},&x_i<a\\1,&a<x_i<b\\1-\dfrac{x_i-b}{M},&x_i>b\end{cases}

至此，所有数据都变成了极大型指标，接下来对数据进行标准化处理，之前我们介绍了很多种方法来标准化，这里我们使用向量归一化法。

y_{ij}=\frac{{x_{ij}}}{\sqrt{\sum\limits_{i=1}^{m}x^2{}_{ij}}}\left(1\leq i\leq m,1\leq j\leq n\right)

接下来就是TOPSIS的核心内容，我们需要计算理想最优解和理想最劣解，然后计算不同指标距离最优解和最裂解的距离（相对贴近度）来进行评分：

S_i^*=\sqrt{\sum_{j=1}^n\left(x_{ij}-x_j^*\right)^2}\quad S_i^-=\sqrt{\sum_{j=1}^n\left(x_{ij}-x_j^-\right)^2}

相对贴近度：C_i^*=\frac{S_i^-}{S_i^-+S_i^*}

例题 河流水质评价

水质测量表.xlsx

现要求使用TOPSIS法对每个河流的水质情况进行排序。

# 河流水质评价
import numpy as np
import pandas as pd

def normalization(data: np.array, type: list):
    """标准化，根据不同的类型，使用不同的标准化方法

    Parameters
    ----------
    data : np.array
        处理的数据
    type :  list
        标准化类型，对应每一列的类型。min：极小型，max：极大型，range：区间型，mid：中间型
    """
    for i in range(len(type)):
        if type[i] == 'min':
            # 极小型，最大值减去各值
            data[:, i] = data[:, i].max() - data[:, i]
        elif type[i] == 'max':
            # 极大型不变
            pass
        elif type[i][0] == 'range':
            # 区间型，获取a,b值，然后进行标准化
            a = type[i][1]
            b = type[i][2]
            M = max(a - data[:, i].min(), data[:, i].max() - b)
            for j in range(len(data[:, i])):
                if data[j, i] < a:
                    data[j, i] = (data[j, i] - a) / M
                elif data[j, i] > b:
                    data[j, i] = 1 - (data[j, i] - b) / M
                else:
                    data[j, i] = 1
        elif type[i][0] == 'mid':
            # 中间型，获取中间值，然后进行标准化
            mid = type[i][1]
            M = max(abs(data[:, i]-mid))
            data[:, i] = 1 - (data[:, i] - mid) / M
        else:
            raise ValueError('type must be min, max, range, mid')
    # 向量归一化
    return abs(data) / np.sqrt(np.sum(abs(data) ** 2, axis=0))

def create_array(i):
    return np.array([i['含氧量(ppm)'], i['PH值'], i['细菌总数(个/mL)'], i['植物性营养物量(ppm)']])

def ideal_solution(v: np.array):
    """理想解法
    """
    # 求理想解，为每一列的最大值
    ideal = np.max(v, axis=0)
    # 求负理想解，为每一列的最小值
    negative_ideal = np.min(v, axis=0)
    # 求每个向量到理想解的距离
    distance_ideal = np.sqrt(np.sum((v - ideal) ** 2, axis=1))
    # 求每个向量到负理想解的距离
    distance_negative_ideal = np.sqrt(
        np.sum((v - negative_ideal) ** 2, axis=1))
    # 求每个向量到理想解的距离与到负理想解的距离的比值
    ratio = distance_negative_ideal / \
        (distance_ideal + distance_negative_ideal)
    return ratio

# 读取数据
data_pd = pd.read_excel('水质测量表.xlsx', sheet_name='Sheet1')
data = data_pd.to_dict(orient='records')
print(data)

# 构建矩阵
m = np.array(list(map(create_array, data)))

# 标准化
m = normalization(m, ['max', ['mid', 7], 'min', ['range', 10, 20]])

# 打印计算结果，保留两位小数，不使用科学计数法
np.set_printoptions(precision=2, suppress=True)

# 理想解法
res = ideal_solution(m)

for i in range(len(data)):
    data[i]['评分'] = res[i]

data_sorted = sorted(data, key=lambda x: x['评分'], reverse=True)
for i in range(len(data_sorted)):
    print(data_sorted[i]['河流'], data_sorted[i]['评分'])

排序结果为：

T 0.6529779260140586
I 0.6357386964756354
O 0.6197806950761303
P 0.5617656502001313
K 0.5484878572629038
M 0.5353536350682877
L 0.5331457294360594
A 0.5033203435158075
E 0.49261415942593856
J 0.49239464964527435
H 0.4847933538204373
N 0.4808595081063982
G 0.4671278748215505
R 0.4670343349581396
C 0.45199296375275866
F 0.4430550337717477
D 0.4361082313866471
Q 0.4112678087637759
B 0.3925029704900129
S 0.3339513984752299

熵权法

我们在信息论中学过信息熵的概念：信息熵反映了信息不确定度的大小，不确定性越大，信息熵就越大。所谓熵权法就是根据判断指标的不确定度来赋权。例如一个指标的数据十分地分散，那么它的权重就应该很高，而另一个指标的数据分散程度很低，例如所有候选人的身高都是180，那么这个判断指标就变得没有意义，所以权重就应该很低。如何来计算各个指标的权重呢，这就是熵权法的目标。

熵权法的赋权步骤如下：

1. 标准化

   在TOPSIS法中我们讨论了如何对数据进行标准化，这里就不再阐述。

2. 计算概率，为求熵做准备：

   \begin{aligned}p_{ij}&=\frac{Y_{jj}}{\sum\limits _{i=1}^nY_{ij}},i=1,\cdots,n,j=1,\cdots,m\\\end{aligned}

3. 求各指标的信息熵：

   E_j=-\frac{\sum\limits^n_{i=1}p_{ij}\ln p_{ij}}{\ln n}

4. 定义信息效用值d_j=1-E_j，信息效用值越大，代表权重越大。

5. 将信息效用值归一化，然后得到熵权：

   w_j=\frac{d_j}{\sum\limits^n_{j=1}d_j}

下面我们结合一道例题来详细解释：

例题 医院科室考核

题目来源：熵权法评价估计详细原理讲解 - 知乎 (zhihu.com)，作者：子木

某医院为了提高自身的护理水平，对拥有的11个科室进行了考核，考核标准包括9项整体护理，并对护理水平较好的科室进行奖励。下表是对各个科室指标考核后的评分结果。

医院评测.xlsx

但是由于各项护理的难易程度不同，因此需要对9项护理进行赋权，以便能够更加合理的对各个科室的护理水平进行评价。

# 医院科室考核，使用熵权法
import numpy as np
import pandas as pd

def entropy_weight_method(data: np.array):
    """熵权法计算权重

    Parameters
    ----------
    data : np.array
        输入的数据矩阵

    Returns
    -------
    np.array
        权重向量
    """
    # 1. 数据标准化
    # 这里都是极大值指标，所以都是正向指标，只需要归一化即可
    data = np.apply_along_axis(normalize, 0, data)
    # print(data)
    # 2. 计算熵值
    data = data / data.sum(axis=0)
    # print(data)
    # 3. 计算熵值
    entropy = -np.nansum(data * np.log(data), axis=0)/np.log(data.shape[0])
    # print(entropy)
    # 4. 计算权重
    weight = (1 - entropy) / np.sum(1 - entropy)
    return weight

def normalize(v):
    """极差变换法进行归一化

    Parameters
    ----------
    v : np.array
        输入的向量

    return
    ------
    np.array
        归一化后的向量
    """
    # 求极差
    r = np.max(v, axis=0) - np.min(v, axis=0)
    # 极差变换
    v = (v - np.min(v, axis=0)) / r
    return v

def create_matrix(data: dict):
    return np.array([data['X1'], data['X2'], data['X3'], data['X4'], data['X5'], data['X6'], data['X7'], data['X8'], data['X9']])

if __name__ == "__main__":
    data_pd = pd.read_excel("医院评测.xlsx", sheet_name="Sheet1")
    data = data_pd.to_dict(orient="records")
    data_matrix = np.array(list(map(create_matrix, data)))
    # print(data_matrix)
    # print(entropy_weight_method(data_matrix))
    # 计算评分
    weight = entropy_weight_method(data_matrix)
    score = np.dot(data_matrix, weight)
    print(score)
    for i in range(len(data)):
        data[i]['分数'] = score[i]
    data_sorted = sorted(data, key=lambda x: x['分数'], reverse=True)
    for i in range(len(data_sorted)):
        print(f"第{i+1}名：{data_sorted[i]['科室']}，分数：{data_sorted[i]['分数']}")

结果：

第1名：F，分数：98.0100557213313
第2名：J，分数：97.80841297792979
第3名：K，分数：97.02126899780376
第4名：I，分数：95.96929203365586
第5名：E，分数：95.84064938103302
第6名：A，分数：95.70696209543694
第7名：H，分数：95.17203465759951
第8名：C，分数：93.17273780796536
第9名：B，分数：93.14062354088951
第10名：D，分数：92.77037549333703
第11名：G，分数：90.20508544958052

熵权法的缺点：

• 忽略了指标本身重要程度，有时确定的指标权数会与预期的结果相差甚远，同时熵值法不能减少评价指标的维数，也就是熵权法符合数学规律具有严格的数学意义，但往往会忽视决策者主观的意图；
• 如果指标值的变动很小或者很突然地变大变小，熵权法用起来有局限。

模糊综合评价法

在我们分析数据的时候，经常会碰见一般、好、优秀等这些模糊的评价，我们的处理方法就是按等级进行赋值。这叫做去模糊化。

但是反过来，我们提供一大堆翔实有效的数据就好了吗，就比如说我想给一款游戏打分，我是要具体的落实在1-10这个分数上，通过一些数据处理，让我们的结论别的不再那么清楚（落实到分数上），这样别人阅读的时候能很快的抓到重点。

如何用科学的方法做出模糊的评价？这就是我们今天要讨论的模糊综合评价法。

模糊综合评价法，便是应用模糊系统的原理，从多个因素对被评判事物的隶属度等级状况进行综合评判的方法。模糊综合评价法的一般步骤如下：

首先需要我们定义隶属函数，例如在年龄的这个评价分布中，分有年轻、中年、年老，但具体30岁算年轻还是年老众说纷纭。我们只需定义三个函数来处理不同年龄对三个称谓的隶属度：

\mu_A\left(x\right)=\begin{cases}1,0<x<20\\\dfrac{40-x}{20},20\leq x\leq40\\0,40<x<150\end{cases}

\mu_B\left(x\right)=\begin{cases}0,0<x<20\\\dfrac{x-20}{20},20\leq x\leq40\\1,40<x<50\\\dfrac{70-x}{20},50\leq x\leq70\\0,70<x<150\end{cases}

\mu_C\left(x\right)=\begin{cases}0,0<x<50\\\dfrac{x-50}{20},50\leq x\leq70\\1,70<x<150\end{cases}

其中\mu_A、\mu_B、\mu_C分别对应年轻、中年、年老的隶属函数。

上述年轻隶属函数可以用以下几种方式来表示：

1. zadch表示法

   A=\frac{\mu_A(x_1)}{x_1}+\frac{\mu_A(x_2)}{x_2}+\frac{\mu_A(x_3)}{x_3}

   注意这里的分号不是真的做了分式运算。

2. 序偶表示法

   A=\{(x_1,\mu_A(x_1)),((x_2,\mu_A(x_2)),((x_3,\mu_A(x_3))\}

3. 向量表示法

   A=(\mu_A(x_1),\mu_A(x_2),\mu_A(x_3))

确定隶属函数的最优方法是模糊统计法，但在数学建模竞赛中这种方法并不适用。在竞赛中常用指派法，也就是人为的定义隶属函数。

在模糊评价模型中，需要引入三个集合：

• 因素集（评价指标集）：U=\{u_1,u_2,...,u_n\}
• 评语集（评价的结果）：V=\{v_1,v_2,...,v_n\}
• 权重集（指标的权重）：A=\{a_1,a_2,...,a_n\}

我们将结合一个实例进行讲解：

例题 露天煤矿开采方案选取

煤矿方案.xlsx

现有5个边坡设计方案，求解最优的方案。

首先找到5个指标的隶属函数：

1. 可采矿量：越多越好，用偏大型：

   \mu_A(x)=\frac x {8800}

2. 基建投资：越少越好，用偏小型：

   \mu_A(x)=1-\frac x {8000}

3. 采矿成本，根据实际情况,其中a_1=5.5，a_2=8.0：

   \mu_c\left(x\right)=\begin{cases}1,&0\le x\le a_1\\\dfrac{a_2-x}{a_2-a_1},&a_1\le x\le a_2\\0,&a_2<x\end{cases}

4. 不稳定费用的隶属函数，越低越好：

   \mu_D(x)=1-\frac x{200}

5. 净现值的隶属函数：

   \mu_E=\frac{x-50}{1500-50}

然后计算每个数值对应的隶属度：

# 露天煤矿方案选取
import numpy as np
import pandas as pd

class MembershipFunctions:

    def muA(self, x):
        return x/8800

    def muB(self, x):
        return 1 - x/8000

    def muC(self, x):
        if 0 <= x <= 5.5:
            return 1
        elif 5.5 < x <= 8:
            return (8-x)/(8-5.5)
        else:
            return 0

    def muD(self, x):
        return 1 - x/200

    def muE(self, x):
        return (x-50)/(1500-50)

    def process_data(self, v: list):
        res = []
        for i in v:
            if i['项目'] == '可采矿量/万吨':
                res.append(list(map(self.muA, list(i.values())[1:])))
            elif i['项目'] == '基建投资/万元':
                res.append(list(map(self.muB, list(i.values())[1:])))
            elif i['项目'] == '采矿成本/(元/吨)':
                res.append(list(map(self.muC, list(i.values())[1:])))
            elif i['项目'] == '不稳定费用/万元':
                res.append(list(map(self.muD, list(i.values())[1:])))
            elif i['项目'] == '净现值/万元':
                res.append(list(map(self.muE, list(i.values())[1:])))
            else:
                raise Exception('项目名称错误')
        res = np.array(res)
        return res

if __name__ == '__main__':
    data_pd = pd.read_excel('煤矿方案.xlsx', sheet_name='Sheet1')
    data = data_pd.to_dict(orient='records')
    data_matrix = MembershipFunctions().process_data(data)
    print(data_matrix)

得到隶属度矩阵：

R=\begin{bmatrix}{} 0.5341 & 0.7614 & 0.6705 & 1.0000 & 0.8636\\ 0.3750 & 0.3125 & 0.3375 & 0.1500 & 0.2500\\ 1.0000 & 0.7600 & 1.0000 & 0.4000 & 0.4800\\ 0.8500 & 0.7500 & 0.8000 & 0.0000 & 0.2000\\ 1.0000 & 0.4483 & 0.6552 & 0.0000 & 0.0345 \end{bmatrix}

设权重矩阵为：A=(0.25,0.20,0.20,0.10,0.25)

最后得到综合评价矩阵：B=A\cdot R

B=\begin{bmatrix}{} 0.7435 & 0.5919 & 0.6789 & 0.3600 & 0.3905 \end{bmatrix}

		data_matrix = MembershipFunctions().process_data(data)
    w = np.array([0.25, 0.2, 0.2, 0.1, 0.25])
    print(w)
    # w乘以矩阵
    res = np.dot(w, data_matrix)
    print(res)

故选择方案1最佳。

灰色关联分析法

有时候我们想知道到底是哪个指标影响总指标最严重，我们假设以及知道某一个指标可能是与其他的某几个因素相关的，那么我们想知道这个指标与其他哪个因素相对来说更有关系，而哪个因素相对关系弱一点，依次类推，把这些因素排个序，得到一个分析结果，我们就可以知道我们关注的这个指标，与因素中的哪些更相关。这时候我们需要用到灰色关联分析法。

使用灰色关联分析的方法如下：

1. 首先找到母序列，就是我们能反映系统特征的序列，类似于因变量：

   Y=[y_1,y_2,...,y_n]^T

2. 确定子序列，是一个矩阵，类似于自变量。

3. 进行数据预处理，对于母序列来说：

   \widetilde{y_{k}}=\frac{y_{k}}{\overline{y_{i}}},\overline{y_{i}}=\frac{1}{n}\sum_{k=1}^{n}y_{k}

   对于子序列来说：

   \stackrel{\sim}{x_{ki}}=\frac{x_{ki}}{\overline{x_{i}}},\overline{x_{i}}=\frac{1}{n}\sum_{k=1}^{n}x_{ki}\big(i=1,2,\cdots,m\big)

4. 计算灰色关联系数

   先计算两极最小差和最大差：

   a=\min_i\min_k|x_0(k)-x_i(k)|,\quad b=\max_i\max_k|x_0(k)-x_i(k)|

   构造：

   \xi_{i}\left(k\right)=y\left(x_{0}\left(k\right),x_{i}\left(k\right)\right)=\frac{a+\rho b}{\left|x_{0}\left(k\right)-x_{i}\left(k\right)\right|+\rho b}

   其中\rho为分辨系数，取0.5。

5. 计算关联度

   r_i=\frac1n\sum_{k=1}^n\xi_i\left(k\right)=\frac1n\sum_{k=1}^ny\Big(x_0\left(k\right),x_i\left(k\right)\Big)

下面对给出一道例题：

例题 什么影响了GDP？

这是某国2016-2019年份的产值和生产总值的数据表，请分析工业产值和农业产值哪个影响国民生产总值更大？

用Python求出关联系数：

if __name__ == '__main__':
    data = np.array([[55, 65, 75, 100], [24, 38, 40, 50], [10, 22, 18, 20]])
    # 数据预处理，每个数据除以每一行的平均数
    data = data / data.mean(axis=1).reshape((-1, 1))
    # 关联系数
    x0 = data[0]
    xi = data[1:]
    # 计算关联系数
    a = min(abs(x0 - xi).min(axis=1))
    b = max(abs(x0 - xi).max(axis=1))
    print(a, b)

求得：a=0.0116，b=0.3758。

代入获得矩阵：

\delta=\begin{bmatrix}{} 0.6605 & 0.6509 & 0.8924 & 0.8749\\ 0.5508 & 0.3540 & 1.0000 & 0.4976 \end{bmatrix}

求平均值得：

r=[0.76966578,0.60058464]

故x_1占比更高，所以工业产值影响国民生产总值更高。

# 灰色关联分析

import numpy as np

def correlation_coefficient_function(x0, xi, a, b, rho=0.5):
    return (a + rho * b) / (abs(x0 - xi) + rho * b)

if __name__ == '__main__':
    data = np.array([[55, 65, 75, 100], [24, 38, 40, 50], [10, 22, 18, 20]])
    # 数据预处理，每个数据除以每一行的平均数
    data = data / data.mean(axis=1).reshape((-1, 1))
    # 关联系数
    x0 = data[0]
    xi = data[1:]
    # 计算关联系数
    a = min(abs(x0 - xi).min(axis=1))
    b = max(abs(x0 - xi).max(axis=1))
    corr_matrix = np.array(
        [correlation_coefficient_function(x0, xi, a, b) for xi in xi])
    # 计算关联度
    corr = corr_matrix.mean(axis=1)
    print(corr)

主成分分析法 PCA

当我们的变量太多时，我们进行评价就会变得异常的麻烦，那么我们能不能找到一种方法，尽可能的减少需要分析的目标，或者用关系紧密的变量合成新的变量，之后我们就能用较少的综合指标的来分别代表各个变量中的各类信息。这就需要主成分分析法。

PCA(Principal Component Analysis)，即主成分分析方法，是一种使用最广泛的数据降维算法。PCA的主要思想是将n维特征映射到k维上，这k维是全新的正交特征也被称为主成分，是在原有n维特征的基础上重新构造出来的k维特征。

主成分分析法的原理很复杂，需要用到方方面面的知识，但是其应用比较简单。

我们假设有n个样本，p个指标，可以构成大小为n\times p的样本矩阵：

x=\begin{bmatrix}x_{11}&x_{12}&\cdots&x_{1p}\\x_{21}&x_{22}&\cdots&x_{2p}\\\vdots&\vdots&\ddots&\vdots\\x_{n1}&x_{n2}&\cdots&x_{np}\end{bmatrix}=(x_1,x_2,...,x_p)

然后假设我们能想到一组新的变量z_1,z_2,...,z_m(m\le p)，满足：

\begin{cases}z_1={l}_{11}x_1+l_{12}x_2+\cdots+l_{1p}x_p\\z_2=l_{21}x_1+l_{22}x_2+\cdots+l_{2p}x_p\\\vdots\\z_m=l_{m1}x_1+l_{m2}x_2+\cdots+l_{mp}x_p\end{cases}

l_{ij}的确定原则如下：

1. z_i和z_j线性无关；
2. z_1是与x_1,x_2,...,x_p的一切线性组合中方差最大者；
3. z_2与z_1不相关的x_1,x_2,...,x_p的一切线性组合中方差最大者；
4. 新变量指标z_1,z_2,...z_m称为第一主成分、第二主成分……第m主成分。

基本的计算步骤如下：

1. 标准化处理

   计算均值和标准差，得到标准化数据：

   X_{ij}=\frac{x_{ij}-\overline{x_j}}{S_j}

2. 计算样本相关系数矩阵：

   R=\begin{bmatrix}r_{11}&r_{12}&\cdots&r_{1p}\\r_{21}&r_{22}&\cdots&r_{2p}\\\vdots&\vdots&\ddots&\vdots\\r_{n1}&r_{n2}&\cdots&r_{np}\end{bmatrix}

   其中：

   r_{ij}=\frac{1}{n-1}\sum_{k=1}^{n}\Big(X_{ki}-\overline{X_{i}}\Big)\Big(X_{kj}-\overline{X_{j}}\Big)

3. 计算R的特征值和特征向量。

4. 计算中主成分贡献率和累计贡献率：

   \text{贡献率}\alpha_i\frac{\lambda_i}{\sum_{k=1}^p\lambda_k}(i=1,2,\ldots,p)\quad\\\text{累计贡献率}\sum G\frac{\sum_{k-1}^i\lambda_k}{\sum_{k=1}^p\lambda_k}(i=1,2,\ldots,p)

5. 写出主成分

   第i个主成分:

   F_{i}=a_{1i}X_{1}+a_{2i}X_{2}+\ldots+a_{pi}X_{p}(i=1,2,\ldots,m)

6. 对于某个主成分而言，指标前面的系数越大，代表对于该主成分的影响越大。

7. 利用主成分的结果进行后续的分析：聚类分析、回归分析。