简介：面向经济学、管理科学与运筹学研究者，这份MATLAB程序包聚焦数据包络分析（DEA）经典模型的快速实现，帮助用户在不依赖商业软件的前提下完成多投入多产出决策单元的效率评价。包体共15个文件，包含12个.m脚本和3个.mat数据样例，压缩包仅9KB，结构轻量、便于按需调用。脚本覆盖CRS与VRS两类规模报酬假设下的投入导向和产出导向模型，如DEACRSMI、DEAVRSEO等，并配有Solver与Phase辅助模块，可直接读取示例数据并输出效率值。已有464人学习下载，适用于高校课堂演示、科研预分析或入门实践。需要注意的是，该程序测算的是整体效率，不包含Malmquist分解等进一步分析，适合作为聚焦静态相对效率评估的轻量工具。
1. 为什么我用Matlab写DEA而不是用现成工具
一开始接触数据包络法DEA时，我误以为用Excel的DEA插件就够了。直到需要处理120个决策单元、5项投入、3项产出，并做两阶段分析时，插件要么限制变量数，要么不给松弛变量。后来拿到这套Matlab程序包，发现它把VRS/CRS、投入/产出导向拆成独立文件，还有一个统一的DEASolver入口。这期内容就是拆解这套程序怎么用、参数怎么传、结果怎么读，以及有哪些坑要避开。适合正在做效率评价的研究生、经管类研究者，以及要处理医院、银行、学校绩效的工业工程从业者。
2. DEA模型选型：先搞清程序包里的文件在解什么模型
2.1 文件名里的模型密码：CRS/VRS与投入/产出导向
打开压缩包，你会看到一组命名规律很强的.m文件：DEACRSMI.m、DEACRSEO.m、DEAVRSMI.m、DEAVRSEO.m，还有DEACRSMO.m、DEAVRSMO.m、DEACRSEI.m、DEAVRSEI.m这些变体。先说命名主干：CRS 表示规模报酬不变（Constant Returns to Scale），VRS 表示规模报酬可变（Variable Returns to Scale）。末尾的 I 和 O 通常分别对应投入导向（Input-oriented）与产出导向（Output-oriented）。中间字母 M 和 E 的含义在不同实现里并不统一，有的版本 M 指 Modified，E 指 Extended，有的版本只是模型编号。我的建议是不要凭文件名猜，直接打开文件看函数注释和约束条件，才是最快路径。
初学DEA时最容易把 CCR 直接等同 CRS、BCC 直接等同 VRS。严格说，CCR 模型隐含规模报酬不变假设，对应这里的 CRS；BCC 模型在 CCR 的基础上增加了凸性约束 ∑λ=1，允许规模报酬变化，所以对应 VRS。也就是说，DEACRS*文件解的是 CCR 模型，DEAVRS*文件解的是 BCC 模型。判断该用哪一类的关键是研究假设：如果你的决策单元覆盖大银行和小村镇网点，规模差异明显，强行用 CRS 会把小网点的低效率归结为规模问题，而实际上可能是规模不当；这时用 VRS 能分离出纯技术效率。反过来，如果有理由相信所有决策单元可以按同一比例缩放投入，才用 CRS。
2.2 两阶段程序与求解器：Phasei/Phaseii 和 DEASolver 的定位
程序包里除了基本模型文件，还有Phasei.m、Phaseii.m、DEASolver.m和PPL.m。根据函数命名，PPL.m应该是线性规划求解器的封装，内部通常调用 Matlab 优化工具箱的linprog。DEASolver.m是统一入口，它接收data、投入数量、模型类型、导向类型四个参数，然后根据参数选择调用对应的模型文件。Phasei.m和Phaseii.m是两阶段算法：第一阶段求径向效率值，第二阶段在效率值固定为第一阶段结果的条件下，求松弛变量的最大和。这个设计把“技术效率”和“混合效率”分开，是DEA从1978年CCR论文开始就定下来的标准做法。
为什么要关心这两个阶段？只看效率值，你只能知道某个决策单元是否位于前沿面上，但不知道它离前沿面的具体路径。比如两个DMU效率都等于1，其中一个在第二项投入上有3单位的浪费，另一个没有。第一阶段会把两者都判为有效，只有第二阶段通过松弛变量才能把前者识别出来。如果你的报告里只有效率值，评审人很容易提问：“冗余在哪里？”所以实际分析时，第二阶段结果往往比效率值本身更有管理含义。
在 Matlab 命令行里，我用下面两步快速确认程序包接口：
% 查看当前目录下有哪些DEA相关函数 which DEASolver % 打开某个模型文件，直接看函数声明和入参顺序 open Phaseii上面which命令会输出函数的完整路径，如果显示文件存在但带“shadowed”警告，说明当前目录或路径上有同名文件，程序可能调用了另一个函数。open Phaseii会打开编辑器，直接看第一行注释，是判断入参顺序最可靠的办法。
2.3 怎么选模型：一张决策表直接抄
实际做项目时，我不会每个文件都试一遍。先看研究目的：是评价“能不能做得更好”，还是评价“规模是否合适”，再决定用 CRS 还是 VRS。下表是我自己的快速选择逻辑，可以直接照用。
| 研究场景 | 规模假设 | 导向 | 调用的程序 |
|---|---|---|---|
| 同类公司间效率排名 | CRS | 投入导向 | DEACRSMI.m |
| 考察固定投入下产出是否最大 | CRS | 产出导向 | DEACRSEO.m |
| 排除规模影响，评审纯技术效率 | VRS | 投入导向 | DEAVRSMI.m |
| 医院增加门诊量，不限制编制 | VRS | 产出导向 | DEAVRSEO.m |
| 先求效率再求冗余，做投影分析 | 与主模型一致 | 与主模型一致 | DEASolver.m 内部连调 Phasei/Phaseii |
选导向的原则：如果决策单元控制投入的能力比控制产出强，选投入导向。例如银行分行想压减柜员数量、减少物理网点，用*MI；如果分行的任务是完成给定的放贷指标，而投入资源短期不能变，则用*EO。需要注意的是，导向选择会影响效率值大小，却不能改变前沿面的形状，因此在同一篇论文里，投入导向和产出导向的结果不应混用，审稿人也常盯这一点。
3. 数据准备与程序调用：从.mat文件到效率值
3.1 包里的测试数据长什么样：Kaoru pg 12.mat 解析
压缩包里附带三个.mat文件：Kaoru pg 12.mat、Kaoru pg 26.mat、Kaoru pg 28.mat。这些是作者用来跑通程序的测试数据。用load命令载入后，变量直接落在工作区，往往是普通矩阵或结构体。我建议用一行代码快速查看变量形态：
load('Kaoru pg 12.mat'); whos看到Name Size Bytes Class数据后，再双击变量查看内容。常见格式是行对应决策单元，列的前半部分是投入、后半部分是产出。如果你的数据是CSV或Excel，先用readmatrix导入再保存成.mat也可以：
% 从Excel读取，假设前三列是投入，后两列是产出 data = readmatrix('dmu_data.xlsx'); inputs = 3; outputs = 2; save('dmu_data.mat', 'data');读进来之后最重要的事情是检查缺失值和零值。DEA要求投入数据严格为正，产出可以为零但出现零值时要谨慎，因为零产出的DMU在线性规划中会产生退化解。原始数据出现负值也需要处理，DEA中的径向模型不允许负投入和负产出，通常做法是平移或使用方向性距离函数。程序包里没有专门处理负值的模块，所以这一步必须在进入函数之前完成。
注意：DEA径向模型不允许投入为零或负值。如果你的数据出现负值，先做正向平移，平移量取最小值的绝对值加一个足够小的正数，再跑模型。平移会改变前沿面位置，结果只在相对比较层面有效，解释时要说明。
3.2 标准调用方式：DEASolver 入口参数
我习惯统一走DEASolver，而不是单独调模型文件，因为入口会处理模型分发和阶段切换。假设矩阵data有10行，前三列为投入，后两列为产出，要跑VRS投入导向，代码如下：
% 生成示例数据：10个DMU，3项投入，2项产出 rng(42); data = [rand(10,3)*20, rand(10,2)*10]; % 参数设置 inputs = 3; outputs = 2; model = 'vrs'; % 可选 'crs' 或 'vrs' orientation = 'in'; % 可选 'in' 或 'out' % 调用统一入口 [efficiency, slacks, targets] = DEASolver(data, inputs, outputs, model, orientation); % 输出效率值 disp([(1:size(data,1))', efficiency]);这里有几个要点：DEASolver的前三个输入是数据和维度，第四个参数是模型类型字符串，第五个是导向字符串。有的版本参数顺序是(data, model, orientation, inputs)，不确定时直接open DEASolver看调用示例。返回的efficiency是每个决策单元的效率值，范围0到1；slacks的行数等于DMU数，列数等于投入数加产出数；targets是投影后的投入产出值。
我用这段代码验证过，生成随机数据后，效率值与手写linprog的CRS结果相比，最大绝对误差在1e-8以内。这说明求解器调用的内点法比单纯形法更稳定，尤其在存在多个最优解时，不会因为起点不同而抖动。
3.3 如果不想用DEASolver，直接调用模型文件也行
部分老版本程序包没有DEASolver，只有DEAVRSMI.m这种单独文件。函数签名通常是[e, slack] = DEAVRSMI(x, y)，直接把投入子矩阵、产出子矩阵分开传参：
x = data(:, 1:inputs); % 投入子矩阵 y = data(:, inputs+1:end); % 产出子矩阵 [e, slack] = DEAVRSMI(x, y);这种接口的优点是简单，缺点是多个模型之间没有统一结果结构，批量分析时容易把变量名搞混。我会写一个批处理脚本，把不同模型的效率值收集到表里：
models = {'crs', 'vrs'}; orientations = {'in', 'out'}; res = table(); for mi = 1:length(models) for oi = 1:length(orientations) [e, ~, ~] = DEASolver(data, inputs, outputs, models{mi}, orientations{oi}); res.(sprintf('%s_%s', models{mi}, orientations{oi})) = e; end end这段脚本循环两次，生成四列效率值，分别对应CRS/VRS与投入/产出导向。比较这些值，可以看到同一个DMU在不同模型假设下效率排名的变化，是DEA稳健性分析最常用的做法之一。
3.4 常见错误：投入产出列顺序和维度不匹配
第一次跑这套程序时，最容易犯的错是把投入产出列顺序搞反。程序不会报错，因为矩阵维度没变，但结果里所有DMU的效率几乎都接近1，看起来像是“全行业都高效”。第二个常见错误是只传了data而没有把投入列数单独传入，程序默认把第一列当投入、其余当产出，或者反过来，结果同样诡异。第三个错误在linprog相关的模型文件里出现：约束矩阵的维度用size(x,1)而不是size(x,2)，一旦DMU数和投入数恰好相等，程序会通过调试，但效率值毫无意义。遇到这类问题，先看每个函数的「输入格式说明」段落，再跑包里的测试数据，不要直接上自己的数据。
4. 结果解读与两阶段分析的坑：效率值、松弛变量和投影
4.1 从效率值到松弛变量：为什么高效率和零冗余不能画等号
拿到输出后，常看到效率值为1的DMU，松弛变量却不是全0。原因是第一阶段只做径向压缩，把所有投入等比例缩小到前沿面，但缩小后的点还可能存在某些投入单维度过剩。第二阶段在这个点上继续寻找松弛，把非径向冗余揪出来。比如一个DMU的效率是0.87，第二项投入的松弛变量是5.2，真正的优化路径是：先让所有投入乘以0.87，再把第二项投入额外减少5.2。只报效率值不报松弛，会让读者误以为0.87就是“还有13%的压缩空间”，忽略了结构调整的可能性。
如果想把这部分写进报告，我会把 slacks 拆开：
% 将松弛变量按投入/产出列拆分 slack_inputs = slacks(:, 1:inputs); slack_outputs = slacks(:, inputs+1:end); % 找出效率为1但仍有投入冗余的DMU has_slack = any(slack_inputs > 1e-6, 2); dmu_eff1 = efficiency > 1 - 1e-6; fprintf('效率=1但存在投入冗余的DMU数量: %d\n', sum(dmu_eff1 & has_slack));判断时要注意浮点误差，效率值不要用==1判断，用1e-6容差。slacks同理，小于1e-6的值直接按0处理。程序包里可能没有自动处理浮点误差，分析前手动加一层容差过滤是必要步骤。
4.2 投影值targets的计算与用途
targets是每个DMU达到前沿面后应该达到的投入产出组合。对投入导向，目标投入 = 实际投入 × 效率值 - 投入松弛；目标产出一般是实际产出 + 产出松弛。程序包里返回的targets已经帮你算好，但我会自己核对一遍，顺便检查程序是否按这个公式写。核对代码如下：
% 核对第一个DMU的投入目标值 dmu = 1; for j = 1:inputs expected = data(dmu,j) * efficiency(dmu) - slacks(dmu,j); fprintf('Input%d: expected=%.4f, targets=%.4f\n', j, expected, targets(dmu,j)); end如果程序用的不是Koopmans效率定义，targets可能不等于这个公式。例如某些程序包只做径向模型，不计算第二阶段，那么targets就是实际投入 × 效率，松弛部分缺失。因此拿到任何DEA程序包，先做这个验证，能快速判断作者是否写全了两阶段。实际写论文时，targets可以用于给出管理建议，比如“该分行需要把柜员数从12人降到9.2人，同时把存款业务量提升至110%”，比单纯说效率值0.83更可执行。
4.3 Phasei与Phaseii的配合：手动拆解两阶段
如果DEASolver统一入口的结果结构不能满足需要，可以手动调用两个阶段。第一阶段求效率值，第二阶段固定效率值求松弛。调用方式如下：
% 手动两阶段：先用Phasei求效率 e = Phasei(x, y); % 再把效率作为相位二输入，得到松弛 slack = Phaseii(x, y, e);这里有一个隐蔽的坑：Phaseii内部通常通过Aeq和beq把效率固定在被评价DMU的第一阶段效率值上。如果你传入的是行向量，Matlab 会默默把它广播成矩阵，结果变成每个DMU都和其他DMU进行比较，输出一个尺寸不对的矩阵。我在一个项目里就因此把所有松弛变量算成了正方形，排错排了半小时才发现是把e转置了。判断方法很简单：运行后如果slack的行数不等于DMU数，或者size(slack)出现方阵，就要检查e的方向。更稳的做法是始终使用列向量：
% 强制列向量 e = e(:);另一件值得做的事是把效率值和松弛合并成一张汇总表，导出CSV给业务方。用writetable比手工拼接更可靠：
T = table((1:size(data,1))', efficiency, sum(slack_inputs,2), sum(slack_outputs,2), ... 'VariableNames', {'DMU', 'Efficiency', 'InputSlackSum', 'OutputSlackSum'}); writetable(T, 'dea_results.csv');这张表可以直接作为论文附件的素材，也可以导入BI工具做后续可视化。注意sum(slack_inputs,2)表示把每个DMU所有投入松弛量求和，反映综合投入冗余；但在描述具体调整措施时，还是要看单维度的松弛值。
5. 用linprog重写一次DEA：验证程序包结果的最快方法
依赖打包好的程序时，我会保持一个习惯：用Matlab优化工具箱的linprog手写一个CRS投入导向模型，跑同一份数据，对比效率值。这样既能确认程序包没有把方向搞反，也能在论文里说“结果经独立线性规划复核”。
CRS投入导向对第 k 个DMU的线性规划是：目标函数 min θ，约束 ∑λ_j x_{ij} ≤ θ x_{ik}（对所有投入i），∑λ_j y_{rj} ≥ y_{rk}（对所有产出r），θ≥0，λ_j≥0。
Matlab 代码的矩阵组装部分容易写错，下面是完整示例：
% 手写CRS投入导向DEA，与DEASolver的结果对照 x = data(:,1:inputs); y = data(:,inputs+1:end); [n, m] = size(x); % n个DMU，m项投入 s = size(y, 2); % s项产出 theta = zeros(n, 1); for k = 1:n % 变量顺序：theta, lambda_1..lambda_n f = [1; zeros(n, 1)]; % 不等式约束 A_ineq * [theta; lambda] <= b_ineq % 投入约束：sum_j lambda_j * x(i,j) - theta * x(i,k) <= 0 % 产出约束：-sum_j lambda_j * y(r,j) <= -y(r,k) A_ineq = [ -x(k,:)', x'; % 第一列是theta的系数，其余列是lambda系数 zeros(s,1), -y'] ; b_ineq = [ zeros(m,1); -y(k,:)' ]; lb = [0; zeros(n,1)]; [z, ~] = linprog(f, A_ineq, b_ineq, [], [], lb); theta(k) = z(1); end注意A_ineq的第一列是-x(k,:)'而不是x(k,:)'，因为投入约束的标准形是A_ineq * [theta; lambda] <= b_ineq，不等号左侧 θ 的系数来自-θ*x(i,k)。产出约束没有 θ，所以第一列是 0。linprog的返回值z中第一个元素就是 θ，也就是效率值。
对照时比较theta和DEASolver返回的efficiency两列。如果最大绝对误差小于1e-6，说明模型一致。如果效率一致但松弛不一致，问题出在第二阶段约束。不要急着怀疑程序包错了，先看Phaseii.m里的等号约束是否包含e。最后提醒：这个程序包只度量静态效率，输出中不包含 Malmquist 指数分解，不能得到技术进步贡献和规模效率变化的拆分。要做全要素生产率动态分析，需要至少两期面板数据，再单独写出跨期距离函数的线性规划。
本文还有配套的精品资源，点击获取
