本人建模小白。。。详细模型思路可参考《清风数学建模交流》坚持学习,加油!!!
优劣解距离法(TOPSIS)
(处理多方案打分的问题)(优点:数据多与少都可以分析,在系统数据资料较少和条件不满足统计要求的情况下,更具有实用性。
缺点:要求样本数据集具有时间序列特性,只是对评判对象的优劣做出鉴别,并不反映绝对水平,具有“相对评价”的全部缺点)
算法步骤:
第一步: 将原始矩阵正向化
第二步: 正向化矩阵标准化
第三步: 计算得分并归一化
实战案例(清风视频)
题目:评价下表中20条河流的水质情况。
注:含氧量越高越好(极大型指标),PH值越接近7越好(中间型指标),细菌总数越少越好(极小型指标),植物性营养物量介于10~20之间最佳,超过20或低于10均不好(区间型指标)
河流 | 含氧量 | PH 值 | 细菌总数 | 植物性营养物量 |
---|---|---|---|---|
A | 4.69 | 6.59 | 51 | 11.94 |
B | 2.03 | 7.86 | 19 | 6.46 |
C | 9.11 | 6.31 | 46 | 8.91 |
D | 8.61 | 7.05 | 46 | 26.43 |
E | 7.13 | 6.50 | 50 | 23.57 |
F | 2.39 | 6.77 | 38 | 24.62 |
G | 7.69 | 6.79 | 38 | 6.01 |
H | 9.30 | 6.81 | 27 | 31.57 |
I | 5.45 | 7.62 | 5 | 18.46 |
J | 6.19 | 7.27 | 17 | 7.51 |
K | 7.93 | 7.53 | 9 | 6.52 |
L | 4.40 | 7.28 | 17 | 25.30 |
M | 7.46 | 8.24 | 23 | 14.42 |
N | 2.01 | 5.55 | 47 | 26.31 |
O | 2.04 | 6.40 | 23 | 17.91 |
P | 7.73 | 6.14 | 52 | 15.72 |
Q | 6.35 | 7.58 | 25 | 29.46 |
R | 8.29 | 8.41 | 39 | 12.02 |
S | 3.54 | 7.27 | 54 | 3.16 |
T | 7.44 | 6.26 | 8 | 28.41 |
python计算结果图
第一步: 将原始矩阵正向化
最常见的四种指标:
指标名称 | 指标特点 | 例子 |
---|---|---|
极大型(效益型)指标 | 越大(多)越好 | 成绩、GDP增速、企业利润 |
极小型(成本型)指标 | 越小(少)越好 | 费用、坏品率、污染程度 |
中间型指标 | 越接近某个值越好 | 水质量评估时的PH值 |
区间型指标 | 落在某个区间最好 | 体温、水中植物性营养物量 |
所谓的将原始矩阵正向化,就是要将所有的指标类型统一转化为极大型指标。(转化的函数形式可以不唯一)
1>极小型指标转化为极大型指标
2>中间型指标转化为极大型指标
3>区间型指标转化为极大型指标
# 第一步: 将原始矩阵正向化
def dataDirection(data):
# 极小型指标转化为极大型指标
_max_ = data['细菌总数'].max()
data['细菌总数'] = _max_ - data['细菌总数']
# 中间型指标转化为极大型指标
x_best = 7
data['PH值'] = abs(data['PH值'] - x_best)
M = data['PH值'].max()
data['PH值'] = 1 - data['PH值'] / M
# 区间型指标转化为极大型指标
a, b = 10, 20
max1 = data['植物性营养物量'].max()
min1 = data['植物性营养物量'].min()
M = max(a - min1, max1 - b)
for i in range(len(data)):
if data.iloc[i, 3] < a:
data.iloc[i, 3] = 1 - (a - data.iloc[i, 3]) / M
elif a <= data.iloc[i, 3] <= b:
data.iloc[i, 3] = 1
else:
data.iloc[i, 3] = 1 - (data.iloc[i, 3] - b) / M
return data
第二步: 正向化矩阵标准化
# 正向化矩阵标准化
def standardization(data):
# 遍历每一列
for i in range(data.shape[1]):
# 求每一列元素的平方和并开根号
a = np.sqrt(sum(data.iloc[:, i] ** 2))
data.iloc[:, i] = data.iloc[:, i] / a
return data
第三步: 计算得分并归一化
# 计算得分并归一化
def normalization(data):
scores = []
length_max = []
length_min = []
_max_ = data.max()
_min_ = data.min()
for i in range(data.shape[0]):
tp1 = tp2 = 0
for j in range(data.shape[1]):
tp1 += ((_max_[j] - data.iloc[i, j]) ** 2)*w[j]
tp2 += ((_min_[j] - data.iloc[i, j]) ** 2)*w[j]
tp1 = np.sqrt(tp1)
tp2 = np.sqrt(tp2)
length_max.append(tp1)
length_min.append(tp2)
scores.append(length_min[i] / (length_max[i] + length_min[i]))
return scores / sum(scores)
指标带有权重的情况
上面计算的得分默认各指标权重都是1,但实际情况下权重不一定都是1,此时便可利用层次分析法计算各个指标的权重,得到加权后的得分。
指标带有权重的计算方法
根据如下计算公式可知只需要在第三步函数中加上权重值即可(权值通过AHP计算得到)
完整代码
# coding=gbk
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
# 第一步: 将原始矩阵正向化
def dataDirection(data):
# 极小型指标转化为极大型指标
_max_ = data['细菌总数'].max()
data['细菌总数'] = _max_ - data['细菌总数']
# 中间型指标转化为极大型指标
x_best = 7
data['PH值'] = abs(data['PH值'] - x_best)
M = data['PH值'].max()
data['PH值'] = 1 - data['PH值'] / M
# 区间型指标转化为极大型指标
a, b = 10, 20
max1 = data['植物性营养物量'].max()
min1 = data['植物性营养物量'].min()
M = max(a - min1, max1 - b)
for i in range(len(data)):
if data.iloc[i, 3] < a:
data.iloc[i, 3] = 1 - (a - data.iloc[i, 3]) / M
elif a <= data.iloc[i, 3] <= b:
data.iloc[i, 3] = 1
else:
data.iloc[i, 3] = 1 - (data.iloc[i, 3] - b) / M
return data
# 正向化矩阵标准化
def standardization(data):
# 遍历每一列
for i in range(data.shape[1]):
# 求每一列元素的平方和并开根号
a = np.sqrt(sum(data.iloc[:, i] ** 2))
for j in range(data.shape[0]):
data.iloc[j, i] = data.iloc[j, i] / a
return data
# 计算得分并归一化
def normalization(data):
scores = []
length_max = []
length_min = []
_max_ = data.max()
_min_ = data.min()
for i in range(data.shape[0]):
tp1 = tp2 = 0
for j in range(data.shape[1]):
tp1 += (_max_[j] - data.iloc[i, j]) ** 2 # 此处可以根据其他如熵权法确定权重再*w[j]即可
tp2 += (_min_[j] - data.iloc[i, j]) ** 2
tp1 = np.sqrt(tp1)
tp2 = np.sqrt(tp2)
length_max.append(tp1)
length_min.append(tp2)
scores.append(length_min[i] / (length_max[i] + length_min[i]))
return scores / sum(scores)
# 绘图
def drawPic(df):
x = df.index
y = df.得分
plt.barh(x, y, height=0.5, align='center')
plt.show()
if __name__ == '__main__':
file_path = '../datas/data1.xlsx'
# 从excel文件中读取数据
data = pd.read_excel(file_path, index_col='河流', sheet_name='Sheet1')
# 第一步: 将原始矩阵正向化
data1 = dataDirection(data)
# 第二步: 正向化矩阵标准化
data1 = standardization(data1)
# 第三步: 计算得分并归一化
scores = np.array(normalization(data1))
df_empty = pd.DataFrame(scores, columns=['得分'], index=data.index)
df = pd.concat([data1, df_empty], axis=1)
# 按得分降序排序
df.sort_values(by=['得分'], ascending=False, inplace=True)
df.to_excel('../datas/data2.xlsx', float_format='%.4f')
# 绘制图像
drawPic(df)
层次分析法(AHP)
层次分析法主要就是构建出判断矩阵求解出权重,但感觉主观性有些强…代码用到的数据取自数学建模清风视频。本代码可根据自己设定的判断矩阵生成最终的权值数据,再用Excel进一步计算结果即可(目标层采用特征值法生成权重,可再添加另两个方法,最好三种方法都做一下)
核心要点
假设方案层数量有N个,准则层包含的指标有M个,则需要构建的判断矩阵的个数为(1 + M)
计算方法
原始数据
景色 | 苏杭 | 北戴河 | 桂林 |
---|---|---|---|
苏杭 | 1 | 2 | 5 |
北戴河 | 1/2 | 1 | 2 |
桂林 | 1/5 | 1/2 | 1 |
1>算术平均法求权重
第一步: 将判断矩阵按照列归一化;
第二步: 将归一化的列相加(按行求和);
第三步: 将相加后得到的向量中的每个元素除以 n 即可得到权重向量.
方法2:几何平均法求权重
第一步: 将A的元素按照行相乘得到一个新的列向量
第二步: 将新的向量的每个分量开n次方
第三步: 对该列向量进行归一化即可得到权重向量
方法3:特征值法求权重
直接将特征向量归一化即可求得特征向量。
求的结果为:最大特征值为 3.0055,
一致性比例 CR = 0.0053 对应的 特征向量:[-0.8902,-0.4132,-0.1918]
归一化后得到权重向量 : [0.5954,0.2764,0.1283]
最终得到如下图所示的权重矩阵:
完整代码
# coding=gbk
import numpy as np
import pandas as pd
import xlrd
class AHP:
# 相关信息的传入和准备
def __init__(self, arr):
# 记录矩阵相关信息
self.arr = arr
# 记录矩阵大小
self.n = arr.shape[0]
# 初始化RI值,用于一致性检验
self.RI_list = [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]
# 矩阵的特征值和特征向量
self.eig_val, self.eig_vector = np.linalg.eig(self.arr)
# 矩阵的最大特征值
self.max_eig_val = np.max(self.eig_val)
# 矩阵最大特征值对应的特征向量
self.max_eig_vector = self.eig_vector[:, np.argmax(self.eig_val)].real
# 矩阵的一致性指标CI
self.CI_val = (self.max_eig_val - self.n) / (self.n - 1)
# 矩阵的一致性比例CR
self.CR_val = self.CI_val / (self.RI_list[self.n - 1])
# 一致性判断
def test_consist(self):
# 打印矩阵的一致性指标CI和一致性比例CR
print("判断矩阵的CI值为:" + str(self.CI_val))
print("判断矩阵的CR值为:" + str(self.CR_val))
# 进行一致性检验判断
if self.n == 2: # 当只有两个子因素的情况
print("仅包含两个子因素,不存在一致性问题")
else:
if self.CR_val < 0.1: # CR值小于0.1,可以通过一致性检验
print("判断矩阵的CR值为" + str(self.CR_val) + ",通过一致性检验")
return True
else: # CR值大于0.1, 一致性检验不通过
print("判断矩阵的CR值为" + str(self.CR_val) + "未通过一致性检验")
return False
# 算术平均法求权重
def cal_weight_by_arithmetic_method(self):
# 求矩阵的每列的和
col_sum = np.sum(self.arr, axis=0)
# 将判断矩阵按照列归一化
array_normed = self.arr / col_sum
# 计算权重向量
array_weight = np.sum(array_normed, axis=1) / self.n
# 返回权重向量的值
return array_weight
# 几何平均法求权重
def cal_weight__by_geometric_method(self):
# 求矩阵的每列的积
col_product = np.product(self.arr, axis=0)
# 将得到的积向量的每个分量进行开n次方
array_power = np.power(col_product, 1 / self.n)
# 将列向量归一化
array_weight = array_power / np.sum(array_power)
# 返回权重向量的值
return array_weight[::-1]
# 特征值法求权重
def cal_weight__by_eigenvalue_method(self):
# 将矩阵最大特征值对应的特征向量进行归一化处理就得到了权重
array_weight = self.max_eig_vector / np.sum(self.max_eig_vector)
# 返回权重向量的值
return array_weight
if __name__ == "__main__":
index = np.array(['a1', 'a2', 'a3', 'a4', 'a5', 'a6', 'a7', 'a8']).T
b1 = pd.read_excel('../../data/判断矩阵.xlsx', index_col=0)
# 算术平均法求权重
weight1 = AHP(b1).cal_weight_by_arithmetic_method()
# 几何平均法求权重
weight2 = AHP(b1).cal_weight__by_geometric_method()
# 特征值法求权重
weight3 = AHP(b1).cal_weight__by_eigenvalue_method()
# 一致性判断
print(AHP(b1).test_consist())
temp = np.hstack((weight1, weight2))
temp = np.hstack((temp, weight3)).reshape(3, 8).T
# 打开文件(若读取xlsx出错降低库的版本就行...)
work_book = xlrd.open_workbook('../../data/判断矩阵.xlsx')
# 获取工作簿所有sheet表对象名称
sheets_name = work_book.sheet_names()
temp1 = []
temp2 = []
temp3 = []
for i in range(1, 9):
b = pd.read_excel('../../data/判断矩阵.xlsx', index_col=0, sheet_name=sheets_name[i]).values
# 一致性判断
print(i)
print(AHP(b).test_consist())
# 算术平均法求权重
weight1 = AHP(b).cal_weight_by_arithmetic_method()
temp1.append(weight1)
# 几何平均法求权重
weight2 = AHP(b).cal_weight__by_geometric_method()
temp2.append(weight2)
# 特征值法求权重
weight3 = AHP(b).cal_weight__by_eigenvalue_method()
temp3.append(weight3)
temp = np.hstack((temp, temp3))
temp = np.hstack((temp, temp2))
temp = np.hstack((temp, temp1))
df_empty = pd.DataFrame(temp,
columns=['算术平均法(1)', '几何平均法(2)', '特征值法(3)', '合资品牌(3)', '自主品牌(3)', '新势力品牌(3)', '合资品牌(2)',
'自主品牌(2)', '新势力品牌(2)', '合资品牌(1)', '自主品牌(1)', '新势力品牌(1)'],
index=['a1', 'a2', 'a3', 'a4', 'a5', 'a6', 'a7', 'a8'])
df_empty.index.name = '各指标'
df_empty.to_excel('../../ans0/data0.xlsx', encoding='utf-8', float_format='%.4f')
df = pd.read_excel('../../ans0/data0.xlsx', index_col=0)
scores = []
sc = ['算术平均法(1)', '几何平均法(2)', '特征值法(3)']
column_names = df.columns
for i in range(3, 12):
if 3 <= i <= 5:
temp1 = sum(df[sc[2]] * df[column_names[i]])
elif 6 <= i <= 8:
temp1 = sum(df[sc[1]] * df[column_names[i]])
else:
temp1 = sum(df[sc[0]] * df[column_names[i]])
scores.append(temp1)
print(scores)
灰色关联分析法(GRA) 参考博客
(处理多方案打分的问题)(优点:原理简单,能同时进行多个对象评价,计算快捷,结果分辨率高、评价客观,具有较好的合理性和适用性,实用价值较高
缺点:只能国赛,只能反映各评价对象内部的相对接近度,并不能反映与理想的最优方案的相对接近程度(我的理解是只能判断已有方案那个最好,不能得出最优方案)
基本思想
首先根据序列曲线几何图形的相似程度来判断其联系是否紧密。曲线越接近,相应的序列之间关联度越大,反之越小。灰色关联度分析对于一个系统发展变化态势提供了量化的度量,非常适合动态历程分析。其基本思想是通过确定参考数据列和若干个比较数据列的几何形状相似程度来判断其联系是否紧密,它反映了曲线间的关联程度。
通常可以运用此方法来分析各个因素对于结果的影响程度,也可以运用此方法解决随时间变化的综合评价类问题,其核心是按照一定规则确立随时间变化的母序列,把各个评估对象随时间的变化作为子序列,求各个子序列与母序列的相关程度,依照相关性大小得出结论
算法步骤:
第一步:确定分析数列;
确定反映系统行为特征的参考数列和影响系统行为的比较数列。反映系统行为特征的数据序列,称为参考数列。影响系统行为的因素组成的数据序列,称为比较数列。
第二步: 变量的无量纲化;
预处理:先求出每个指标的均值,再用该指标中的每个元素都除以该均值。(约定俗成)
第三步:计算关联系数;
第四步:计算关联度;
第五步:关联度排序
实战案例
下表为某地区国内生产总值的统计数据(以百万元计),问该地区从2000年到2005年之间哪一种产业对GDP总量影响最大
年份 | 国内生产总值 | 第一产业 | 第二产业 | 第三产业 |
---|---|---|---|---|
2000 | 1988 | 386 | 839 | 763 |
2001 | 2061 | 408 | 846 | 808 |
2002 | 2335 | 422 | 960 | 953 |
2003 | 2750 | 482 | 1258 | 1010 |
2004 | 3356 | 511 | 1577 | 1268 |
2005 | 3806 | 561 | 1893 | 1352 |
步骤1:确立母序列
在此需要分别将三种产业与国内生产总值比较计算其关联程度,故母序列为国内生产总值。若是解决综合评价问题时则母序列可能需要自己生成,通常选定每个指标或时间段中所有子序列中的最佳值组成的新序列为母序列。
步骤2:无量纲化处理
在此采用均值化法,即将各个序列每年的统计值与整条序列的均值作比值,可以得到如下结果:
年份 | 国内生产总值 | 第一产业 | 第二产业 | 第三产业 |
---|---|---|---|---|
2000 | 0.731959 | 0.836101 | 0.682761 | 0.743906 |
2001 | 0.758837 | 0.883755 | 0.688458 | 0.787780 |
2002 | 0.859720 | 0.914079 | 0.781229 | 0.929152 |
2003 | 1.012518 | 1.044043 | 1.023735 | 0.984725 |
2004 | 1.235641 | 1.106859 | 1.283331 | 1.236269 |
2005 | 1.401325 | 1.215162 | 1.540486 | 1.318167 |
步骤3:计算每个子序列中各项参数与母序列对应参数的关联系数( i 表示子序列个数,这里就是三个产业;k 表示每个序列包含的样本数,这里就是年份数目(可以理解成指标数目))
可得到如下结果:
年份 | S01(t) | S02(t) | S03(t) |
---|---|---|---|
2000 | 0.475145 | 0.658636 | 0.892228 |
2001 | 0.429863 | 0.573289 | 0.767955 |
2002 | 0.635577 | 0.546182 | 0.576630 |
2003 | 0.752048 | 0.898480 | 0.775266 |
2004 | 0.422378 | 0.665686 | 1.000000 |
2005 | 0.335584 | 0.403502 | 0.531718 |
步骤4:计算关联度
第一产业 0.508432
第二产业 0.624296
第三产业 0.757300
通过比较三个子序列与母序列的关联度可以得出结论:该地区在2000年到2005年期间的国内生产总值受到第三产业的影响最大。
完整代码
# coding=gbk
import pandas as pd
data = pd.read_excel('../datas/原始数据.xlsx', index_col='年份')
# 数据无量纲化
mean_col = data.mean()
for i in range(data.shape[1]):
data.iloc[:, i] = data.iloc[:, i] / mean_col[i]
# 提取参考队列和比较队列
ck = data.iloc[:, 0]
cp = data.iloc[:, 1:4]
for i in range(cp.shape[1]):
data.iloc[:, i + 1] = abs(ck[:] - cp.iloc[:, i])
cp = cp = data.iloc[:, 1:4]
# 求最值
min_ = cp.min(axis=0).min()
max_ = cp.max(axis=0).max()
p = 0.5
ksi = (min_ + p * max_) / (cp + p * max_)
r = ksi.mean()
# 5、关联度降序排序
ans = r.sort_values(ascending=False)
print(ksi)
print(r)
print(ans)