简介面向经济学、管理科学与运筹学研究者这份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 ScaleVRS 表示规模报酬可变Variable Returns to Scale。末尾的 I 和 O 通常分别对应投入导向Input-oriented与产出导向Output-oriented。中间字母 M 和 E 的含义在不同实现里并不统一有的版本 M 指 ModifiedE 指 Extended有的版本只是模型编号。我的建议是不要凭文件名猜直接打开文件看函数注释和约束条件才是最快路径。初学DEA时最容易把 CCR 直接等同 CRS、BCC 直接等同 VRS。严格说CCR 模型隐含规模报酬不变假设对应这里的 CRSBCC 模型在 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个DMU3项投入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到1slacks的行数等于DMU数列数等于投入数加产出数targets是投影后的投入产出值。我用这段代码验证过生成随机数据后效率值与手写linprog的CRS结果相比最大绝对误差在1e-8以内。这说明求解器调用的内点法比单纯形法更稳定尤其在存在多个最优解时不会因为起点不同而抖动。3.3 如果不想用DEASolver直接调用模型文件也行部分老版本程序包没有DEASolver只有DEAVRSMI.m这种单独文件。函数签名通常是[e, slack] DEAVRSMI(x, y)直接把投入子矩阵、产出子矩阵分开传参x data(:, 1:inputs); % 投入子矩阵 y data(:, inputs1: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(:, inputs1: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(:,inputs1:end); [n, m] size(x); % n个DMUm项投入 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 指数分解不能得到技术进步贡献和规模效率变化的拆分。要做全要素生产率动态分析需要至少两期面板数据再单独写出跨期距离函数的线性规划。本文还有配套的精品资源点击获取