Matlab：数据包络分析 (DEA) 入门教程

发布时间：2021-09-11 阅读 17303

Stata连享会主页 || 视频 || 推文 || 知乎

温馨提示： 定期清理浏览器缓存，可以获得最佳浏览体验。

New！ lianxh 命令发布了：
随时搜索推文、Stata 资源。安装命令如下：
. ssc install lianxh
详情参见帮助文件 (有惊喜)：
. help lianxh

✌ 课程详情： https://gitee.com/lianxh/Course

⛳ 课程主页： https://gitee.com/lianxh/Course

⛳ Stata 系列推文：

全部 | Stata入门 | Stata教程 | Stata资源 | Stata命令
计量专题 | 论文写作 | 数据分享 | 专题课程
结果输出 | Stata绘图 | 数据处理 | Stata程序
回归分析 | 面板数据 | 交乘项-调节 | IV-GMM
内生性-因果推断 | 倍分法DID | 断点回归RDD | PSM-Matching | 合成控制法
Probit-Logit | 时间序列 | 空间计量 | 分位数回归 | 生存分析 | SFA-DEA
文本分析-爬虫 | Python-R-Matlab | 机器学习
Markdown | 工具软件 | 其它

作者：曹昊煜 (兰州大学)
邮箱：caohy19@lzu.edu.cn

编者按：本文部分公式来自下文，特此致谢！
Source：成刚. 数据包络分析方法与MaxDEA软件[M]. 知识产权出版社, 2014. -PDF-

0. 写在前面

这篇推文主要为大家介绍一些数据包络分析 (Data Envelopment Analysis, DEA) 基本模型的原理和编程。我本科的时候就初步接触了 DEA 模型，并使用该方法撰写了第一篇真正意义上的学术论文，对于一个大山深处的二年级本科生而言，这无疑是令人振奋的。随着对计量经济理论的深入了解，我逐渐放弃了对数据包络分析的应用，因为以我的能力难以理解复杂的线性规划式和该方法的经济含义。直到近期连享会开设了相关的课程，才促使我重新审视这一非参数方法并尝试了一些基本模型的编程。仿佛回到了那个炎热的午后，我满怀期待地盯着 Matlab 的窗口，一如当年我紧张地看着 DEAP 的界面。

实际上，DEA 方法的编程实现并不复杂，其原因一方面在于目前已经出现了诸如 DEAP、Maxdea、DEA-slover 等菜单式软件，另一方面在于 Matlab 等软件都包含成熟的求解线性规划的工具箱，而不用我们自己编写单纯形法的底层代码。但对于初学者而言，菜单式软件恰如一个“黑箱”，让人难以琢磨，Matlab 在线性规划中的应用又难以与 DEA 的规划式结合起来。撰写本文的目的正在于介绍一些基本 DEA 模型的原理，并详解其 Matlab 编程思想。熟悉这些基本模型，可以极大地缓解学习高等模型的恐惧心理。

1. 线性规划及其 Matlab 实现

DEA 方法的基本模型的形式以线性规划为主，或者可以通过一些变换将非线性规划转换为线性形式。因此在这一部分中，我们首先为大家介绍线性规划及其在 Matlab 中的求解方法。

1.1 线性规划与对偶

线性函数是在一组线性约束条线的限制下，最优化线性目标函数的问题。线性规划可以求目标函数的最大化或者最小化值，在 Matlab 中，一般求解的是最小值 (如果目标是最大化问题，只需要求相反数即可)，其形式是：

\begin{aligned} m i n_{x} & c^{T} x \\ s . t . & A x \leq b \\ A e q \cdot x = b e q \\ l b \leq x \leq u b \end{aligned}

这一形式的目标函数是一个内积，包含了不等式约束、等式约束和取值范围三种约束条件，是一种最一般的表达方式。求解线性规划的过程中有时需要使用对偶技术降低求解难度，如果线性规划原问题的形式为：

\begin{aligned} m i n_{x} & c^{T} x \\ s . t . & A x \geq b \\ x \geq 0 \end{aligned}

则其对偶形式为：

\begin{aligned} m a x_{x} & b^{T} y \\ s . t . & A^{T} y \leq c \\ y \geq 0 \end{aligned}

推导对偶技术需要注意以下要点：

原问题和对偶问题的优化方向是相反的；
原问题目标函数的决策变量个数等于对偶问题的约束个数；
原问题的约束个数等于对偶问题的决策变量个数；
原问题不等式约束系数矩阵的转置等于对偶问题中不等式约束的系数矩阵；
原问题的目标函数参数等于对偶问题约束式右侧的系数；
原问题不等式约束右侧的系数等于对偶问题目标函数的参数。

初步接触对偶变化时，很难一次性掌握其变化规律，需要多次推导对偶问题，才能保证熟练使用。关于对偶的具体例子可以参考「线性规划对偶问题」。

1.2 线性规划在 Matlab 中的实现

对一个最小化的线性规划问题：

\begin{aligned} m i n_{x} & c^{T} x \\ s . t . & A x \leq b \\ A e q \cdot x = b e q \\ l b \leq x \leq u b \end{aligned}

在 Matlab 中可以使用一条简单的代码求解：

[x,fval]=linprog(c,A,b,Aeq,beq,LB,UB,X0,OPTIONS)

linprog 是单纯形法求解线性规划的命令，c 为目标函数参数向量，A 为约束条件左侧的系数矩阵，b 为约束条件右侧的参数向量，Aeq 和 beq 分别是等式约束的系数矩阵和参数向量，LB 和 UB 同时指定决策变量的上下界，暂时不考虑初值和其他高级选项。该命令的返回值有两个，分别是决策变量的最优值和目标函数的最大值。

再次给出一个具体的示例，假设现在需要求解如下问题：

\begin{aligned} m a x & z = 2 x_{1} + 3 x_{2} - 5 x_{3} \\ s . t . & x_{1} + x_{2} + x_{3} = 7 \\ 2 x_{1} - 5 x_{2} + x_{3} \geq 10 \\ x_{1} + 3 x_{2} + x_{3} \leq 12 \\ x_{1}, x_{2}, x_{3} \geq 0 \end{aligned}

根据 Matlab 中的代码要求，我们首先需要完成两件事：一是将最大化目标函数变为最小化，二是要将除上下界外的不等号全部变化为小于等于号，变换后的问题是：

\begin{aligned} m i n & z = - 2 x_{1} - 3 x_{2} + 5 x_{3} \\ s . t . & x_{1} + x_{2} + x_{3} = 7 \\ - 2 x_{1} + 5 x_{2} - x_{3} \leq - 10 \\ x_{1} + 3 x_{2} + x_{3} \leq 12 \\ x_{1}, x_{2}, x_{3} \geq 0 \end{aligned}

由此可以写出求解该问题的代码：

c = [-2 -3 5];           % 指定目标函数参数
A = [-2 5 -1
     1  3 1];            % 指定不等式约束系数
b = [-10 12]';           % 指定不等式约束右侧参数

Aeq = [1 1 1];           % 指定等式约束系数
beq = 7;                 % 指定等式约束右侧参数

UB = [];                 % 决策变量无上界
LB = [0 0 0]';           % 指定决策变量下界

[x,fval] = linprog(c,A,b,Aeq,beq,LB,UB)    % 单纯形法求解

决策变量的求解结果为 $[6.43, 0.57, 0]^{T}$ ，目标函数的最小化结果为 -14.57。通过这一个简单的例子我们可以熟悉线性规划在 Matlab 中的实现方法，其关键在于将目标函数和约束中的线性表达转换为矩阵表达。

2. 技术效率与 DEA 基本编程思想

2.1 什么是技术效率

技术效率指的是一个生产单元的生产过程达到本行业技术水平的程度。一般来说，技术效率可以使用产出和投入的比例衡量，但这种衡量方式一般仅适用于单投入单产出的情形。如果有 m 种投入和 q 种产出，则可以使用加权方式确定综合的投入产出情况：

\begin{aligned} v = v_{1} x_{1} + v_{2} x_{2} + \dots + v_{m} x_{m} \\ u = u_{1} y_{1} + u_{2} y_{2} + \dots + u_{q} y_{q} \end{aligned}

现在假设存在 k 个决策单元 (DMU)，则每个决策单元的产出投入比可以表示为：

h_{k} = \frac{u_{1} y_{1 k} + u_{2} y_{2 k} + \dots + u_{q} y_{q k}}{v_{1} x_{1 k} + v_{2} x_{2 k} + \dots + v_{m} x_{m k}} = \frac{\sum_{r = 1}^{q} u_{r} y_{r k}}{\sum_{i = 1}^{m} v_{i} y_{i k}}

2.2 投入与产出导向

在径向 DEA 中，无效率往往是通过投入和产出的等比例变化定义的，因此既可以在给定投入的情况下最大化产出，也可以在给定产出的情况下最小化投入。投入导向的 DEA 模型指的是在给定产出的情况下，选择最小化的投入，而产出导向的 DEA 模型则正好相反。

对于不同的规模报酬，不同导向的效率分析结果可能存在一定差异。如果模型是规模报酬不变的模型时 (CRS)，两种导向的效率分析结果是一样的。而在可变规模报酬 (VRS) 模型中，二者是不同的。

在实践中，投入和产出导向的选择没有明确的要求，但是需要符合模型的基本含义。比如，在一项生产活动中已经给定了要素投入的数量，那么我们就不必再选择最小的产出，而是应该在给定投入的情况下最大化产出。

2.3 模型编程的基本步骤

DEA 模型的编程本质上是调用软件自带的工具箱求解线性规划问题，因此，给定相应的公式后，其编程过程遵循一定的步骤：

第一步：确定参数列向量，并将目标函数写成内积形式，如果目标函数是最大化问题，则对参数系数取相反数；
第二步：将所有约束条件重写为符合最小化问题的形式，即不等式均为小于等于；
第三步：使用矩阵表示约束条件右侧的系数矩阵，需要注意正负号和分块矩阵之间的行列对应；
第四步：将约束条件右侧系数写为列向量；
第五步：根据特定语法编写程序。

经过这些特定的步骤，我们可以编写出任意 DEA 基本模型的代码，开始的练习可以在草稿纸上将线性规划式展开，以便将线性规划的求和形式转变为矩阵形式。在本文中，我们将主要对两个 CCR 模型的矩阵形式进行详细推导和介绍。

3. CCR 模型：CRS 径向 DEA

3.1 投入导向的 CRS 模型

首先介绍规模报酬不变的 CRS 模型，使用上一小节中加权产出投入比的概念，如果进一步将技术效率的范围限定在 $[0, 1]$ ，我们就可以的到一个最基本的 DEA 模型：

\begin{aligned} m a x & \frac{\sum_{r = 1}^{q} u_{r} y_{r k}}{\sum_{i = 1}^{m} v_{i} y_{i k}} \\ s . t . & \frac{\sum_{r = 1}^{q} u_{r} y_{r k}}{\sum_{i = 1}^{m} v_{i} y_{i k}} \leq 1 \\ v \geq 0; u \geq 0 \\ i = 1, 2, \dots, m; r = 1, 2, \dots, q; j = 1, 2, \dots, n \end{aligned}

该模型的特征是规模报酬不变的。但由于该形式是一个非线性规划，因此需要将其转化线性规划形式，变换后的结果为：

\begin{aligned} m a x & \sum_{r = 1}^{q} μ_{r} y_{r k} \\ s . t . & \sum_{r = 1}^{q} μ_{r} y_{r j} - \sum_{i = 1}^{m} ν_{i} x_{i j} \leq 0 \\ \sum_{i = 1}^{m} ν_{i} x_{i k} = 1 \\ μ \geq 0; ν \geq 0 \\ i = 1, 2, \dots, m; r = 1, 2, \dots, q; j = 1, 2, \dots, n \end{aligned}

相对于这一形式，我们更关心对偶模型，因为对偶模型的决策变量中包含效率值。上述问题的对偶形式为：

\begin{aligned} m i n & θ \\ s . t . & \sum_{j = 1}^{n} λ_{j} x_{i j} \leq θ x_{i k} \\ \sum_{j = 1}^{n} λ_{j} y_{r j} \geq y_{r k} \\ λ \geq 0 \\ i = 1, 2, \dots, m; r = 1, 2, \dots, q; j = 1, 2, \dots, n \end{aligned}

在对偶规划中， $λ$ 表示 DMU 的线性组合系数，k 表示待评价的 DMU，参数 $θ^{*}$ 即为效率值，其范围在 0 到 1 之间。该模型的含义相当于使用加权方法构造一个不存在的 DMU，其投入不大于待评价的 DMU，产出不小于待评价的 DMU，即 $x = \sum_{j = 1}^{n} λ_{j} x_{i j}$ ， $y = \sum_{j = 1}^{n} λ_{j} y_{r j}$ 。

为了更清晰地说明 CRS 模型的编程方式，我们首先将包含求和符号的公式展开为线性求和形式：

\begin{aligned} m i n & 0 \cdot λ_{1} + 0 \cdot λ_{2} + \dots + 0 \cdot λ_{n} + 1 \cdot θ \\ s . t . & λ_{1} x_{11} + λ_{2} x_{12} + \dots + λ_{n} x_{1 n} \leq θ x_{1 k} \\ λ_{1} x_{21} + λ_{2} x_{22} + \dots + λ_{n} x_{2 n} \leq θ x_{2 k} \\ \dots \\ λ_{1} x_{m 1} + λ_{2} x_{m 2} + \dots + λ_{n} x_{m n} \leq θ x_{m k} \\ λ_{1} y_{11} + λ_{2} y_{12} + \dots + λ_{n} y_{1 n} \geq y_{1 k} \\ λ_{1} y_{21} + λ_{2} y_{22} + \dots + λ_{n} y_{2 n} \geq y_{2 k} \\ \dots \\ λ_{1} y_{q 1} + λ_{2} y_{q 2} + \dots + λ_{n} y_{q n} \geq y_{q k} \\ λ \geq 0 \\ i = 1, 2, \dots, m; r = 1, 2, \dots, q; j = 1, 2, \dots, n \end{aligned}

该模型的目标函数具有一定的迷惑性，很容易让人认为最优化的决策参数只有一个，但实际上决策变量向量为 $[λ_{1}, λ_{2}, \dots, λ_{n}, θ]$ ，一共有 n+1 个参数。我们使用向量和矩阵形式表达这一问题，并将不等式约束条件全部转化为小于等于。

\begin{aligned} m i n & 0 \cdot λ_{1} + 0 \cdot λ_{2} + \dots + 0 \cdot λ_{n} + 1 \cdot θ \\ s . t . & [\begin{array}{c} x_{11} & x_{12} & \dots & x_{1 n} & - x_{1 k} \\ x_{21} & x_{22} & \dots & x_{2 n} & - x_{2 k} \\ ⋮ & ⋮ & \dots & ⋮ & ⋮ \\ x_{m 1} & x_{m 2} & \dots & x_{m n} & - x_{m k} \end{array}] [\begin{array}{c} λ_{1} \\ λ_{2} \\ ⋮ \\ λ_{n} \\ θ \end{array}] \leq [\begin{array}{c} 0 \\ 0 \\ ⋮ \\ 0 \end{array}] \\ [\begin{array}{c} - y_{11} & - y_{12} & \dots & - y_{1 n} & 0 \\ - y_{21} & - y_{22} & \dots & - y_{2 n} & 0 \\ ⋮ & ⋮ & \dots & ⋮ & ⋮ \\ - y_{q 1} & - y_{q 2} & \dots & - y_{q n} & 0 \end{array}] [\begin{array}{c} λ_{1} \\ λ_{2} \\ ⋮ \\ λ_{n} \\ θ \end{array}] \leq [\begin{array}{c} - y_{1 k} \\ - y_{2 k} \\ ⋮ \\ - y_{q k} \end{array}] \\ λ \geq 0 \\ i = 1, 2, \dots, m; r = 1, 2, \dots, q; j = 1, 2, \dots, n \end{aligned}

将该形式与 Matlab 要求的函数形式对应，即可编写投入导向的 CRS 模型代码。使用 31 个省份的医疗资源数据作为示例 (示例数据下载)，分析各省份的医疗资源的投入产出效率。

% 导入数据
data = xleread('dea_data.xlsx',1);

% 投入导向模型
X= data(:,1:2)';  % 投入指标数据，每一列代表每个决策单元的投入数据
Y= data(:,3:4)';  % 产出指标数据，每一列代表每个决策单元的产出数据
n=size(X,2);      % 决策单元数
m=size(X,1);      % 投入指标数
q=size(Y,1);      % 产出指标数

w = [];
for i = 1:n
    f = [zeros(1,n) 1];   % 定义目标函数
    Aeq = [];             % 没有等式约束
    beq = [];
    LB = zeros(n+1,1);    % 指定下界
    UB = [];
       
    A = [X -X(:,i);-Y zeros(q,1)];         % 设定不等式约束
    b = [zeros(m,1);-Y(:,i)];
    
    w(:,i) = linprog(f,A,b,Aeq,beq,LB,UB); % 模型求解
end
    
theta_CCR_input = w(n+1,:)';               % 结果输出

数据中共包含三个变量，前两个变量为投入数据，分别是 “床位数” 和 “卫技人员数”，后三个变量为产出数据，有两个期望产出 “诊疗人次” 和 “入院人数”，一个非期望产出 “医疗废弃物”。在 CRS 模型中，投入和产出导向的效率分析结果相同，所以计算结果与产出导向模型一起分析。

3.2 产出导向的 CRS 模型

产出导向的模型与投入导向的模型设定正好相反，其含义是给定产出的条件下最小化投入，线性规划式为：

\begin{aligned} m i n & \sum_{i = 1}^{n} ν_{i} x_{i k} \\ s . t . & \sum_{r = 1}^{q} μ_{r} y_{r j} - \sum_{i = 1}^{m} ν_{i} x_{i j} \leq 0 \\ \sum_{r = 1}^{q} μ_{r} y_{r k} = 1 \\ μ \geq 0; ν \geq 0 \\ i = 1, 2, \dots, m; r = 1, 2, \dots, q; j = 1, 2, \dots, n \end{aligned}

其对偶模型为：

\begin{aligned} m a x & ϕ \\ s . t . & \sum_{j = 1}^{n} λ_{j} x_{i j} \leq x_{i k} \\ \sum_{j = 1}^{n} λ_{j} y_{r j} \geq ϕ y_{r k} \\ λ \geq 0 \\ i = 1, 2, \dots, m; r = 1, 2, \dots, q; j = 1, 2, \dots, n \end{aligned}

该模型的效率定义为 $1 / ϕ$ ，现在需要将该模型的形式转换为我们所熟悉的形式。首选需要将最大化问题改为最小化，其次需要将不等式约束全部变换为小于等于的形式：

\begin{aligned} m i n & 0 \cdot λ_{1} + 0 \cdot λ_{2} + \dots + 0 \cdot λ_{n} - 1 \cdot ϕ \\ s . t . & λ_{1} x_{11} + λ_{2} x_{12} + \dots + λ_{n} x_{1 n} \leq x_{1 k} \\ λ_{1} x_{21} + λ_{2} x_{22} + \dots + λ_{n} x_{2 n} \leq x_{2 k} \\ \dots \\ λ_{1} x_{m 1} + λ_{2} x_{m 2} + \dots + λ_{n} x_{m n} \leq x_{m k} \\ - λ_{1} y_{11} - λ_{2} y_{12} - \dots - λ_{n} y_{1 n} \leq - ϕ y_{1 k} \\ - λ_{1} y_{21} - λ_{2} y_{22} - \dots - λ_{n} y_{2 n} \leq - ϕ y_{2 k} \\ \dots \\ - λ_{1} y_{q 1} - λ_{2} y_{q 2} - \dots - λ_{n} y_{q n} \leq - ϕ y_{q k} \\ λ \geq 0 \\ i = 1, 2, \dots, m; r = 1, 2, \dots, q; j = 1, 2, \dots, n \end{aligned}

其矩阵形式为：

\begin{aligned} m i n & 0 \cdot λ_{1} + 0 \cdot λ_{2} + \dots + 0 \cdot λ_{n} - 1 \cdot ϕ \\ s . t . & [\begin{array}{c} x_{11} & x_{12} & \dots & x_{1 n} & 0 \\ x_{21} & x_{22} & \dots & x_{2 n} & 0 \\ ⋮ & ⋮ & \dots & ⋮ & ⋮ \\ x_{m 1} & x_{m 2} & \dots & x_{m n} & 0 \end{array}] [\begin{array}{c} λ_{1} \\ λ_{2} \\ ⋮ \\ λ_{n} \\ ϕ \end{array}] \leq [\begin{array}{c} x_{1 k} \\ x_{2 k} \\ ⋮ \\ x_{m k} \end{array}] \\ [\begin{array}{c} - y_{11} & - y_{12} & \dots & - y_{1 n} & y_{1 k} \\ - y_{21} & - y_{22} & \dots & - y_{2 n} & y_{2 k} \\ ⋮ & ⋮ & \dots & ⋮ & ⋮ \\ - y_{q 1} & - y_{q 2} & \dots & - y_{q n} & y_{q k} \end{array}] [\begin{array}{c} λ_{1} \\ λ_{2} \\ ⋮ \\ λ_{n} \\ ϕ \end{array}] \leq [\begin{array}{c} 0 \\ 0 \\ ⋮ \\ 0 \end{array}] \\ λ \geq 0 \\ i = 1, 2, \dots, m; r = 1, 2, \dots, q; j = 1, 2, \dots, n \end{aligned}

将该形式与 Matlab 要求的函数形式对应，即可编写产出导向的 CRS 模型代码，其形式为：

% 导入数据
data = xlsread('DEA_data.xlsx',1);

% 产出导向模型
X= data(:,1:2)';  % 投入指标数据，每一列代表每个决策单元的投入数据
Y= data(:,3:4)';  % 产出指标数据，每一列代表每个决策单元的产出数据
n=size(X,2);      % 决策单元数
m=size(X,1);      % 投入指标数
q=size(Y,1);      % 产出指标数

w = [];
for i = 1:n
    f = [zeros(1,n) -1];    % 定义目标函数
    Aeq = [];
    beq = [];               % 无等式约束
    LB = zeros(n+1,1);      % 指定下界
    UB = [];
    
    A = [X zeros(m,1);-Y Y(:,i)];           % 设定不等式约束
    b = [X(:,i)' zeros(1,q)]';
    
    w(:,i) = linprog(f,A,b,Aeq,beq,LB,UB);  % 求解模型
end

w = 1./w;                          % 计算效率
theta_CCR_output = w(n+1,:)';      % 结果输出

计算出两种导向的 DEA 模型后，使用以下命令可以导出计算结果，并与省份名称对应：

name = {'安徽','北京','福建','甘肃','广东','广西','贵州','海南','河北','河南',...
        '黑龙江','湖北','湖南','吉林','江苏','江西','辽宁','内蒙古','宁夏','青海',...
        '山东','山西','陕西','上海','四川','天津','西藏','新疆','云南','浙江',...
        '重庆'};
    
table([theta_CCR_input theta_CCR_output],'RowNames',name)


    安徽     0.94424    0.94424
    北京     0.91743    0.91743
    福建           1          1
    甘肃     0.83383    0.83383
    广东           1          1
    广西     0.79526    0.79526
    贵州     0.96309    0.96309
    海南     0.83313    0.83313
    河北     0.90387    0.90387
    河南     0.88262    0.88262
    黑龙江   0.68934    0.68934
    湖北     0.94694    0.94694
    湖南     0.98402    0.98402
    吉林     0.73382    0.73382
    江苏     0.91013    0.91013
    江西           1          1
    辽宁      0.7232     0.7232
    内蒙古   0.70576    0.70576
    宁夏       0.813      0.813
    青海     0.78401    0.78401
    山东     0.88123    0.88123
    山西     0.61303    0.61303
    陕西     0.82264    0.82264
    上海           1          1
    四川     0.93916    0.93916
    天津     0.98505    0.98505
    西藏     0.68185    0.68185
    新疆           1          1
    云南           1          1
    浙江     0.94798    0.94798
    重庆     0.87968    0.87968

由计算结果可以看出，两种模型的计算结果完全相同。其中福建、广东、江西、上海、新疆、云南六个省份的医疗资源投入产出最有效率，而山西省的效率改进空间最大。

4. BCC 模型：VRS 径向 DEA

4.1 投入导向的 BCC 模型

CCR 模型的假定是规模收益可变，所以该模型的效率测度包含了规模因素。可以在 CCR 模型的基础上增加约束：

\sum_{j = 1}^{n} λ_{j} = 1

即可得到包含规模变化的 BCC 模型，这一约束的作用是令部分生产前沿向内收缩，以得到消除规模因素的纯技术效率。

BCC 模型的线性规划式为：