窄带包络解调python实现_对数据包络分析法DEA的再理解，以及python 实现

上一篇：什么是数据包络分析法 DEA 中，对数据包络分析发进行了初步的了解和认识，对其具体的原理并未深究。本文继续对DEA进行学习，并用python实现。

基础视频了解

视频来自公众号：九转星河

综述

绩效评估是评估组织或个人如何以较少的资源获得较多的产出结果的多属性评估，也称之为成本效益分析。数据包络分析是 A.Charnes, W.W.Copper 和 E.Rhodes 在 1978 年提出的评价多指标输入输出，衡量系统有效性的方法。将属性划分为投入项、产出项(成本型、效益型指标)，不预先设定权重，只关心总产出与总投入，以其比率作为相对效率。数据包络分析有多种模式，主要为：CCR 模式，BBC 模式、交叉模式、A&P 模式。

构建决策单元

在 DEA 中，每个评估的对象称为“决策单元”(记为 DMU)，设共有n，每个决策单元都有 m_1种投入和 m_2 种产出，x_{ij} (i=1,2,...,n; j= 1,2,....,m_2) 表示第 i 个决策单元的第 j项产出，u= (u_1,u_2,....,u_{m1}),v= (v_1,v_2,....,v_{m2}) 分别表示投入产出的权值向量。

建立 DEA效率评估模型(CCR 模式)

为啥叫这个奇怪的名字呢？因为这个模型符号是曲子研究者姓氏的第一个字母。比如 CCR模型是由：A.Charnes , W.W.Cooper 及 E.Rhodes 三人建立的。

定义：决策单元 k的效益评价指数为：

从投资资源的角度来看，在当前产出的水准下，比较投入资源的使用情况，以此作为效益评价的依据，这所模式称之为”投入导向模式“。输出结果 e 则表示的是相对有效性的评价值。我们此处要用的也正是这个值。

现在，我们的求值变成了求 e_k 的最大值，但这种模式不是传统的线性规划模型，难以求得最优解，因此将其线性化后取对偶模型。如下：

min OE_k

这个模型可以理解为，将决策单元 k 的投入、产出表示为其它决策单元的线性组合。

如果有某个(某些)决策单元的产出量达到决策单元 k 的水平(达到，不一定要超过，第二个约束的含义)，而投入量要尽可能的小(第一个约束的含义)，那么就会出现 OE_k< 1 的情况，此时说明，该决策单元还存在资源浪费。

繁殖，如果该决策单元的效率已经是最高了，任何的决策单元都不能使用比他还小的投入获得同样(甚至更多)的产出，这时候 OE_k= 1

为加速模型求解及分析产能效应，为上述方程引入了松弛变量 s_i^-,s_j^+。其中 s_i^- 成为差额变数，表示该决策单元为达到 ”DEA有效“ 应减少的投入量， s_j^+ 成为超额变数，代表为达到”DEA有效“ 应增加的产出量。引入非阿基米德数 \varepsilon，通常设为 10^{-6}，则CCR模型线性规划模型化为：

尽管 \varepsilon 很小，但我们无法获取 s_i^-,s_j^+ 的大小数量级，直接进行计算可能会产生误差，因此求解该问题我们使用分层序列法。第一阶段求解 OE_k 的最小值，第二阶段在 OE_k 已知的情况下求 \displaystyle\sum_{i=1}^{m_1}s_i^- + \displaystyle\sum_{j=1}^{m_2}s_j^+ 的最大值，即：

有效性分析

示例

总结：因为对每个决策单元来说，投入、产出的数值都是确定的，那么求解过程，实际上就是在已有的数据和条件限定下，求解对应的每个标签的权重系数。(所以，计算的结果可能并不是某个确定值，只要是在限定条件内的都有可能是结果)

相当于，求出一个DEA有效的决策单元，然后以此决策单元为基准，对其它DMU进行相对有效评价。

DEA有效的经济学含义

python程序

import gurobipy

import pandas as pd

# 分页显示数据, 设置为 False 不允许分页

pd.set_option('display.expand_frame_repr', False)

# 最多显示的列数, 设置为 None 显示全部列

pd.set_option('display.max_columns', None)

# 最多显示的行数, 设置为 None 显示全部行

pd.set_option('display.max_rows', None)

class DEA(object):

def __init__(self, DMUs_Name, X, Y, AP=False):

self.m1, self.m1_name, self.m2, self.m2_name, self.AP = X.shape[1], X.columns.tolist(), Y.shape[1], Y.columns.tolist(), AP

self.DMUs, self.X, self.Y = gurobipy.multidict({DMU: [X.loc[DMU].tolist(), Y.loc[DMU].tolist()] for DMU in DMUs_Name})

print(f'DEA(AP={AP}) MODEL RUNING...')

def __CCR(self):

for k in self.DMUs:

MODEL = gurobipy.Model()

OE, lambdas, s_negitive, s_positive = MODEL.addVar(), MODEL.addVars(self.DMUs), MODEL.addVars(self.m1), MODEL.addVars(self.m2)

MODEL.update()

MODEL.setObjectiveN(OE, index=0, priority=1)

MODEL.setObjectiveN(-(sum(s_negitive) + sum(s_positive)), index=1, priority=0)

MODEL.addConstrs(gurobipy.quicksum(lambdas[i] * self.X[i][j] for i in self.DMUs if i != k or not self.AP) + s_negitive[j] == OE * self.X[k][j] for j in range(self.m1))

MODEL.addConstrs(gurobipy.quicksum(lambdas[i] * self.Y[i][j] for i in self.DMUs if i != k or not self.AP) - s_positive[j] == self.Y[k][j] for j in range(self.m2))

MODEL.setParam('OutputFlag', 0)

MODEL.optimize()

self.Result.at[k, ('效益分析', '综合技术效益(CCR)')] = MODEL.objVal

self.Result.at[k, ('规模报酬分析', '有效性')] = '非 DEA 有效' if MODEL.objVal < 1 else 'DEA 弱有效' if s_negitive.sum().getValue() + s_positive.sum().getValue() else 'DEA 强有效'

self.Result.at[k, ('规模报酬分析', '类型')] = '规模报酬固定' if lambdas.sum().getValue() == 1 else '规模报酬递增' if lambdas.sum().getValue() < 1 else '规模报酬递减'

for m in range(self.m1):

self.Result.at[k, ('差额变数分析', f'{self.m1_name[m]}')] = s_negitive[m].X

self.Result.at[k, ('投入冗余率', f'{self.m1_name[m]}')] = 'N/A' if self.X[k][m] == 0 else s_negitive[m].X / self.X[k][m]

for m in range(self.m2):

self.Result.at[k, ('差额变数分析', f'{self.m2_name[m]}')] = s_positive[m].X

self.Result.at[k, ('产出不足率', f'{self.m2_name[m]}')] = 'N/A' if self.Y[k][m] == 0 else s_positive[m].X / self.Y[k][m]

return self.Result

def __BCC(self):

for k in self.DMUs:

MODEL = gurobipy.Model()

TE, lambdas = MODEL.addVar(), MODEL.addVars(self.DMUs)

MODEL.update()

MODEL.setObjective(TE, sense=gurobipy.GRB.MINIMIZE)

MODEL.addConstrs(gurobipy.quicksum(lambdas[i] * self.X[i][j] for i in self.DMUs if i != k or not self.AP) <= TE * self.X[k][j] for j in range(self.m1))

MODEL.addConstrs(gurobipy.quicksum(lambdas[i] * self.Y[i][j] for i in self.DMUs if i != k or not self.AP) >= self.Y[k][j] for j in range(self.m2))

MODEL.addConstr(gurobipy.quicksum(lambdas[i] for i in self.DMUs if i != k or not self.AP) == 1)

MODEL.setParam('OutputFlag', 0)

MODEL.optimize()

self.Result.at[k, ('效益分析', '技术效益(BCC)')] = MODEL.objVal if MODEL.status == gurobipy.GRB.Status.OPTIMAL else 'N/A'

return self.Result

def dea(self):

columns_Page = ['效益分析'] * 3 + ['规模报酬分析'] * 2 + ['差额变数分析'] * (self.m1 + self.m2) + ['投入冗余率'] * self.m1 + ['产出不足率'] * self.m2

columns_Group = ['技术效益(BCC)', '规模效益(CCR/BCC)', '综合技术效益(CCR)','有效性', '类型'] + (self.m1_name + self.m2_name) * 2

self.Result = pd.DataFrame(index=self.DMUs, columns=[columns_Page, columns_Group])

self.__CCR()

self.__BCC()

self.Result.loc[:, ('效益分析', '规模效益(CCR/BCC)')] = self.Result.loc[:, ('效益分析', '综合技术效益(CCR)')] / self.Result.loc[:,('效益分析', '技术效益(BCC)')]

return self.Result

def analysis(self, file_name=None):

Result = self.dea()

file_name = 'DEA 数据包络分析报告.xlsx' if file_name is None else f'\\{file_name}.xlsx'

Result.to_excel(file_name, 'DEA 数据包络分析报告')

以上内容主要转自：https://zhuanlan.zhihu.com/p/60853027

另有参考文档：数据包络分析 (点击下载)