简介:这套DEA的MATLAB程序包面向运筹学、经济管理等领域研究者与师生,旨在解决多输入多输出系统相对效率评估问题,内置DEA Solver工具箱,支持CCR、BCC及非标准DEA模型的定义、求解与结果分析。压缩包共21个文件,大小2.54MB,主要包含MATLAB脚本(.m/.asv)、数据文件(.xls/.mat)、使用说明(.txt/.doc)及PDF文献,覆盖模型代码、示例数据与参考文档,便于直接运行和二次开发。已有443人学习下载。借助这份资源,读者可系统学习DEA原理,在MATLAB中快速构建并求解模型,也可根据实际需求自定义输入输出指标,用于企业、医院等组织效率评价;配套PDF包括聚类DEA方法、区域创新政策排名等专题,适合教学演示与进阶研究。注意需安装MATLAB优化工具箱,大型数据集计算时建议调整算法参数以提升效率。
1. 为什么我推荐这份DEA的Matlab程序:不只是能跑,而是能跑明白
做效率评价的人应该都经历过这种纠结:用现成软件点一点就能出结果,但换了指标、改了导向、想批量跑几百个决策单元时,就被卡住了。这份DEA的Matlab程序,内置了dea solver,不是那种只能算固定结果的封装黑盒,而是把CCR、BCC这些核心模型用可读的代码实现出来,跑完还能把松弛变量、权重、投影值导出。我自己拿它复算过教材里的经典数据集,结果跟文献对得上。适合三类人:一是需要用Matlab做研究复现的同学,二是想摆脱点选式工具、把效率评价嵌入自己分析流程的从业者,三是不放心黑盒结果、想逐矩阵验证的谨慎派。
这份资源的核心价值不在于“能跑”,而在于“可改”。很多商业软件把模型封装成图形界面,你只能选菜单,一旦需要自定义约束或者批量处理,就得手动点几百次。而这个程序包把数据读取、模型构建、solver调用、结果解析拆成了独立函数,你只需要按格式准备Excel,然后调用一个入口函数,剩下的交给程序。而且它自带的solver支持多种线性规划求解路径,默认使用Matlab优化工具箱的接口,也可以替换成开源的求解器,这给了二次开发很大的余地。
我是在处理一组某地区医院的运营效率数据时开始用它的。当时要比较30多个决策单元在三年里的效率变化,传统软件导出的结果表格里连中间变量都没有。换了这个程序包之后,我把三年的数据分别灌进去,再写个循环就能得到效率排名和投影值,前后用了不到半天。后来再帮别人看这类问题,我都会先拿出这套程序讲清楚DEA到底在算什么,而不是直接甩一个结果。
2. 先把DEA的求解逻辑拆开:solver在中间扮演什么角色
在动手跑数据之前,必须把DEA的数学结构和solver的关系搞清楚。很多人因为没弄明白这一点,拿着程序包乱试参数,结果算出来的东西自己都不敢信。这一章我会把效率测度怎么化成线性规划、CCR和BCC差在哪、solver怎么接进来一层层讲透。
2.1 效率测度怎么化成线性规划:投入导向与产出导向的选择
DEA的核心思想很简单:给定一堆决策单元,每个都有若干投入和产出,通过构造一个由所有单元组成的生产可能集,判断每个单元是否落在前沿面上。如果落在外边,就得看它离前沿面有多远,这个距离折算成效率分数。关键在于“距离”怎么量:是沿投入方向压缩,还是沿产出方向扩张。这就是投入导向和产出导向的区别。
投入导向求的是“在产出不变的前提下,投入最多能同比例压缩多少”。效率分数小于1时,说明该单元存在投入冗余,可以等比例减少投入却保持产出不变。产出导向则反过来,问“投入不变时,产出最多能同比例扩张多少”。这两种导向写出的线性规划对偶形式不同,但对同一组数据,非有效单元的排名往往不完全一致。实际业务里,如果你的控制变量是投入,比如医院想控制成本、银行想压缩网点开支,就用投入导向;如果产出目标刚性,比如要完成上级下达的出院人次指标,就用产出导向。
在程序里,导向就是入口函数的一个参数。我第一次跑的时候建议先用投入导向,因为它最直观,投影点就是投入压缩后的坐标,方便跟业务人员解释。如果你把参数改成产出导向,你会发现目标函数里输入输出的位置发生了变化,最终的效率分数也会因模型构造不同而略有差异,这不是bug,是定义问题。还有一个容易混淆的点:投入导向和产出导向下的效率分数并不互为倒数,只有当模型是单投入单产出且规模报酬不变时,两者才恰好互为倒数。所以不要试图去转换。
2.2 CCR与BCC:从规模报酬不变到可变,模型约束只差一个等式
CCR模型是1978年提出的,它假设所有决策单元都在最优规模下运营,即规模报酬不变。这意味着投入翻倍,产出也应该翻倍,否则就会被判定为低效。但实际上很多行业存在规模经济或规模不经济,小诊所和大诊所的投入产出关系不是线性的。后来提出的BCC模型加上了凸性约束,允许规模报酬可变,这样测出来的效率排除了规模因素,叫“纯技术效率”。
有趣的是,CCR和BCC在数学表达式上的差别极小。CCR模型里的lambda权重可以是不等于1的非负值,而BCC模型增加了一个约束:lambda求和等于1。就这一个约束,把生产可能集从一个圆锥变成了一个凸包。不要小看这一行约束,它的直接后果是,分别用两个模型跑同一份数据时,有效单元数量完全不同:BCC的有效单元永远不少于CCR,因为凸包内的点更多。如果你的样本里有大量小规模或大规模单元,CCR会把它们的低效归咎于规模不当,而BCC则会把规模因素剥离,让你看到运营管理本身的问题。
读这份程序时,我特别看重它是如何把两个模型放在同一个函数里的。常见做法是设置一个变量,比如model_type,当为'BCC'时加上一行Aeq(即lambda求和为1的等式),而不是复制两套代码。这样改起来不会漏。你可以把两个模型的结果放在一起对比,如果某个决策单元在CCR中无效但在BCC中有效,说明它的无效率主要来自规模因素,这时候再结合规模报酬状态判断该扩大还是收缩规模。程序包通常还会顺便算出规模效率,也就是CCR效率除以BCC效率,这个值能帮助你拆解低效来源。
2.3 DEA solver在Matlab里的两种实现路径:linprog内置与自定义单纯形
DEA本质上是一个线性规划问题,所以求解的关键就是线性规划求解器。Matlab自带的优化工具箱里有linprog函数,专门求解形如 min f'x, A x <= b, Aeq x = beq, lb <= x <= ub 的问题。DEA模型正好可以转换成这个标准形式。例如投入导向CCR模型的线性规划中,目标函数是theta(投入压缩比例),约束包括对每个产出项的不等式:Ylambda >= Y_i,对每个投入项的不等式:Xlambda <= theta*X_i,以及lambda非负约束。只要把theta和lambda拼接成x向量,把对应的系数填进矩阵,linprog就能跑。
但注意,linprog默认使用对偶单纯形法或内点法,它对稀疏矩阵的支持好,但如果你想特殊处理,比如加一个额外的整数约束,或者想输出对偶价格,就需要换求解器。这份程序包里我看到它做了两层封装:第一层是模型层,负责把DEA的约束填入稀疏矩阵;第二层是求解层,统一调用一个solver接口。对象是linprog,但接口是独立的,这意味着你可以把solver函数换成另一个开源求解器,只要返回结构相同。
我一般会在程序里加一段探测代码,先检查Matlab是否安装了优化工具箱。如果没装,就退回到程序自带的一个简化版单纯形算法。那个算法是教材级别的,针对DEA模型的特征做了优化,没有做稀疏矩阵处理,所以速度慢一些,但结果正确。这个后备设计很实用,因为在某些离线环境里,你没有完整工具箱的许可证。程序里solver的选择逻辑一般是:优先用linprog,捕获异常后降级到自带求解器,同时输出警告。这一层我觉得是这份资源最有价值的地方,它没有把求解过程焊死。
2.4 从对偶解看规模报酬:solver给你的不只是效率分数
很多初用者只拿效率分数,其实线性规划的对偶解信息量大得多。CCR模型的对偶变量是投入和产出的影子价格,而BCC模型的对偶中还有一个自由变量,它的符号能直接告诉你在当前单元所在区间的规模报酬状态是递增还是递减。程序里如果solver返回了对偶信息,你就能自动判断规模报酬状态,而不需要单纯依靠比较CCR和BCC的效率值来间接推断。
这个手段在批量分析中很重要。比如你要给几十个决策单元分年度写结论,靠人工去区分“规模报酬递增”“不变”“递减”非常容易看错。程序可以循环读取对偶变量,按符号打标签。代码上,linprog的输出参数中,lambda.lower和lambda.upper对应于不等式约束的对偶,但如果你把theta放在第一个变量,约束的方向需要小心。我建议你对照教材的模型形式,把约束逐行打印出来核对一遍,这个功夫省不掉。一旦对偶解正确了,规模报酬判断就是一行符号判断,比如对偶变量大于0对应递增,小于0对应递减。有了这个,你就能在自动生成的报告里给每个单元标注“建议扩规模”或“建议控规模”。
3. 这份程序包怎么用:目录结构、数据格式与三个关键参数
从这一章开始进入实操。先别急着跑主程序,花十分钟把目录和输入格式看清楚,能帮你省下后面几天排错的时间。这一章我按拿到压缩包之后的正常顺序来写:先看文件结构,再准备数据,最后设置参数。
3.1 文件清单与调用入口
拿到压缩包后,先看目录。一般这类资源会分为几个文件夹:core存放模型构建函数,solver存放线性规划求解接口,utils存放数据读取和结果输出函数,example存放演示脚本。你下载回来的大概率是这些源码和一个说明文档。首先把整个文件夹加入Matlab路径,或者在当前目录下打开。然后看主入口函数,通常是dea_main.m或者run_dea.m。输入是三个变量:X矩阵、Y矩阵、选项结构体。下面给出一个基本调用示例。
% 加载演示数据 load('example_data.mat'); % 包含矩阵X和Y,已按常见格式准备好 % X为n×m矩阵,n个决策单元,m个投入指标 % Y为n×s矩阵,n个决策单元,s个产出指标 options.model = 'CCR'; % 可选 'CCR' 或 'BCC' options.orient = 'input'; % 可选 'input' 或 'output' options.solver = 'auto'; % 'auto' 自动选择linprog或内置solver [te, lambda, slacks, targets] = dea_main(X, Y, options);这里dea_main返回四个核心结果:te是效率分数列向量,lambda是参考集权重矩阵,slacks是松弛变量矩阵,targets是每个决策单元在前沿上的投影值。我一般把结果存成表格后,先看te列,再看targets和原始数据的差值,就能算出投入冗余量和产出不足量。参数说明:options.model决定用CCR还是BCC;options.orient决定优化方向,input代表求投入压缩比例,output代表求产出扩张比例;options.solver设为'auto'时,程序会优先调用linprog,失败时自动降级。
3.2 输入矩阵格式:投入在左、产出在右,缺一不可
最常见的翻车是把投入产出转置。程序默认X是n行×m列,Y是n行×s列。如果你从Excel粘贴时把行当成DMU列,结果会显示维度不匹配;或者更糟的是矩阵反了但维度恰好一致,跑出来的效率值变成对指标的评价,而不是对单元的评价。所以第一步要确认:每一行是一个决策单元,前几列是投入,后几列是产出。程序包一般会自带一个check_input函数,建议调用一下。
% 检查输入维度 check_input(X, Y); % 如果X和Y的行数不一致,这里会报错 % 如果数据有缺失,这里会输出缺失位置数据里不能有空值。如果Excel里某个单元格是空白的,导入后是NaN,DEA的线性规划中有NaN约束会导致solver返回错误状态。常见的做法是把缺失值替换为所在列的均值,或者直接删掉那个决策单元。我一般先做一次完整性统计,如果有超过5%的缺失,宁可用插补,而不是强行删除样本,因为删除样本可能破坏前沿面的形状。
另外要注意投入产出指标的方向。DEA默认投入越小越好,产出越大越好。如果你的指标是反的,比如“损坏率”作为产出但实际上越低越好,那得先做正向化处理,比如取倒数或者用最大值减去原始值。程序里没有自动识别指标属性的功能,这个必须你根据业务自己处理。在数据准备阶段,我还习惯把每一列的单位统一量纲,虽然DEA本身不要求去量纲化,但如果某列数值数量级差异太大,可能造成求解器数值困难。可以用列均值归一化或者保留原始单位,只要保证同列单位一致即可。
3.3 三个关键参数:model_type、orientation、return_to_scale
很多DEA程序包把模型类型和规模报酬分开设置。实际上,CCR天然假设规模报酬不变,BCC天然假设可变,所以部分程序只设model_type就够了。如果你在代码里看到有return_to_scale参数,要注意它和model_type的联动关系:当model_type='CCR'时,return_to_scale必须为'CRS';当model_type='BCC'时,通常是'VRS',但也有少数程序支持'NIRS'(非递增规模报酬)。选错组合会直接报错或得出无意义结果。
导向参数orientation控制目标函数。输入导向的数学意义是求theta的最小值,theta是投入压缩比例;输出导向是求phi的最大值,phi是产出扩张比例。在代码里,theta或phi都放在第一个变量位置,但约束不等式方向相反。程序包的注释里通常会写明定义,例如: % orient: 'input' -> minimize theta subject to Xlambda <= thetaX_i % orient: 'output' -> maximize phi subject to Ylambda >= phiY_i
我建议在第一次运行前,先跑一个只有两个决策单元的迷你例子,手算一遍结果,确认程序的方向定义和你理解的一致。比如两个人,一个人投入10产出5,另一个人投入20产出10,用投入导向CCR,两个人的效率都应该是1,因为第二个人的产出正好同比例于投入,两者都在前沿面上。如果程序算出来不是1,说明你的导向或指标方向反了。这种最小化测试能筛掉一半的配置错误。
3.4 结果输出与保存:别手动复制,让程序写Excel
跑完DEA后,结果在Matlab工作区里还不行,你得输出成业务方看得懂的报表。程序一般提供export_results函数,可以直接写成Excel。常见输出包括四张表:效率表、参考集表、松弛变量表、投影对照表。效率表应该有决策单元名称、技术效率、规模效率(如果同时跑CCR和BCC)、排名。投影对照表则要列出每个决策单元的原始投入、投影投入、冗余量、原始产出、投影产出、不足量,这些是改进建议的依据。
% 导出结果到Excel,格式为表格 result_table = table(dmu_names, te, slack_input, slack_output, ... 'VariableNames', {'DMU', 'TE', 'SlackIn', 'SlackOut'}); writetable(result_table, 'dea_results.xlsx', 'Sheet', '效率结果'); % 投影值矩阵与原始矩阵对比,输出到第二张表 writetable(table(raw_in, target_in, raw_out, target_out, ... 'VariableNames', {'RawIn','TargetIn','RawOut','TargetOut'}), ... 'dea_results.xlsx', 'Sheet', '投影分析');这里要注意,writetable对中文字段名支持不太稳定,如果Excel里出现乱码,建议把VariableNames统一用英文,另建一列中文字段说明。另外,导出的Excel不要覆盖原始数据文件,命名上加上日期。我习惯把每次跑的结果文件名带上参数摘要,比如dea_result_CCR_in_20250101.xlsx,这样后面追溯参数时才不会乱。还有一点,如果你用的是比较老的Matlab版本,writetable可能不支持xlsx格式,需要改成xls或者先安装额外扩展,这一点在说明文档里通常会有提示。
4. 跑通一个真实案例:从Excel数据到效率排名
这一章用一个具体的模拟案例演示完整流程。案例背景是16家社区诊所,我们要评价它们的运营效率。投入选两个指标:医务人员数量和季度运营成本;产出选一个指标:服务人次。数据是我按真实场景构造的,主要为了展示步骤,你在自己项目里换成实际数据即可。
4.1 数据准备与导入
先建一个Excel文件,第一列是决策单元名称,第二列和第三列是投入,第四列是产出。这里我把数据文件命名为clinic_data.xlsx,放在Matlab当前工作目录下。读取时用xlsread,注意txt变量里第一行是表头,所以要跳过去。
% 读取Excel数据,数据区域从第二行开始 [num, txt] = xlsread('clinic_data.xlsx'); dmu_names = txt(2:end, 1); % 第一列是名称 % num的行数是16,列是3(两个投入、一个产出) X = num(:, 1:2); % 投入矩阵:医务人员数、运营成本 Y = num(:, 3); % 产出矩阵:服务人次读取后建议立即打印size确认维度,再调用check_input。如果数据里有空行,xlsread会跳过,导致num的行数和txt不一致,这时需要检查原始文件有没有多余的空行。我习惯在读取后用sum(isnan(X(:)))检查缺失值,数量为0再继续。这个案例里的数据没有缺失,但你自己处理真实数据时,这一步不能省。
接下来设置选项并调用dea_main。因为我关心的是诊所自身的管理改进,而诊所之间规模差异明显,所以我采用BCC模型,投入导向,这样得到的是纯技术效率。
options.model = 'BCC'; options.orient = 'input'; options.solver = 'auto'; [te, lambda, slacks, targets] = dea_main(X, Y, options);执行完以后,先看te。你会发现有些诊所效率是1,有些小于1。这时不要急着下结论,先检查有效单元数量是否合理。如果有效单元超过一半,可能是指标太少或者样本同质化太高。本案例16个单元,两个投入一个产出,有效单元数在4到6个是比较正常的。
4.2 求解与结果解读
程序跑完后,把效率分数按升序排列。我通常会先把结果做成一个表格,显示名称、效率、排名,方便交流。本案例里,排名第一的诊所效率为1,前五名中三名是1,其余两个接近0.96。最后一名效率只有0.61,明显落后。
除了排名,还要看参考集。lambda矩阵每一行表示该决策单元对应的有效参考单元及权重。比如第16个诊所效率0.61,lambda行中第5个和第9个诊所的权重分别为0.4和0.6,说明它的改进方向是学习这两家诊所的投入组合。而有效单元的lambda行里只有自己权重为1,说明它自给自足,定义了前沿面的一部分。
这时我一般会画一个散点图:横轴是运营成本,纵轴是服务人次,把每个诊所的点标出来,效率最高的用实心点连接成前沿面,无效点用空心点,再画一条从原点到投影点的虚线。在Matlab里可以用scatter和plot组合实现。这个图比表格直观,业务方一眼就能看出谁在前沿面上,谁在下方。代码并不复杂,关键是先把有效单元的索引提取出来,然后按成本升序排序连线。
4.3 用结果反查无效单元的改进目标
得到投影值targets后,计算每个投入的冗余量。对第16个诊所,原始投入成本是420,人员数35,投影值成本是310,人员数28,说明在同等的服务人次下,人员冗余5人,成本可以压缩110。从DEA角度讲,这些是“纯技术无效”导致的多余投入。但如果模型是CCR,还存在规模无效,这时就要通过规模效率判断是扩大还是收缩规模。BCC下效率分数小于1就是纯技术无效,规模效率需要CCR和BCC的效率比值获得。
程序里通常会直接给出规模效率值。例如第16诊所CCR效率0.51,BCC效率0.61,则规模效率为0.83,说明即便把管理效率提升到同行水平,规模不经济还会拖累约17个百分点。具体是规模过大还是过小,要看对偶变量,第2.4小节已经讲过方法,这里不再重复。
此外还要注意一个细节:如果targets里某个投入出现了0,那大概率是数据有问题或者模型选择错误。现实中投入不可能降为零,所以出现0时,要检查是否把导向参数写错了,或者投入指标里存在退化列。我通常会做一次敏感性分析:删掉一个前沿单元,看其他单元的排名变化幅度。如果排名剧烈变化,说明这组数据对样本构成很敏感,结论要谨慎使用。这个案例的数据比较平稳,删掉任一个前沿单元后,排名变化都在两位以内,说明结果是稳定的。
5. 避坑:DEA在Matlab里常见的六个翻车现场
跑DEA翻车不可怕,可怕的是跑出结果却不知道错在哪。这一章我整理了六条实战高频坑,全部按现象、原因、解决的逻辑写,你在自己项目里遇到同类型问题时可以直接对照。
5.1 结果全是1,区分度很低
现象:跑完所有决策单元效率都是1,或者绝大多数是1,排名几乎没有意义。
原因:最常见的是样本量太少或投入产出项太多。DEA的经典经验法则是决策单元数量至少要是投入产出指标个数的两倍以上,更严格的建议是三倍。当决策单元数量接近或小于指标数量时,前沿面会被少数单元撑起来,导致大量单元都落在前沿上。另一个原因是你的数据本身高度同质,比如所有单元的投入产出基本同比例变化,这时所有点都落在同一条射线上,CCR模型下自然全是1。
解决:合并相关性高的指标,或者增加样本量。如果没法增加样本,至少把投入指标从5个压缩到2个,用主成分分析提炼综合投入指标。也可以把模型换成BCC,BCC的凸包约束会带来更多有效单元,但区分度仍然有限。更好的办法是检查一下数据是否存在“所有单元都按同一比例缩放”的情况,如果确实如此,建议重新选择指标维度。
5.2 linprog报错Infeasible但数据看起来没问题
现象:调用dea_main时solver返回状态为Infeasible,程序中断,控制台出现大红报错。
原因:DEA模型的线性规划通常是可行的,因为把lambda设为该决策单元自身总可以满足约束。如果出现不可行,十有八九是数据中有负值或缺失值,或者X、Y矩阵的结构被颠倒。还有一种情况是不小心把某一列指标的符号取反了,比如用“成本率”做投入,本应越小越好,但录入了负值或取倒数后变成非凸约束。
解决:先检查矩阵维度,再检查负值。DEA中投入产出必须非负,而且至少有一个决策单元的各指标为正。如果某个产出有零值,某些求解器可能因浮点运算出现问题,这不是模型问题,是数值问题。把零值替换为一个很小的正数,比如1e-6,再试。如果还不行,就把模型降级为CCR,或者用手写单纯形法工具解一次,看看是否是矩阵条件数问题。
5.3 效率为1但松弛变量不为0
现象:某个决策单元效率分数等于1,但程序里slacks矩阵对应行有非零值。
原因:DEA中效率为1只说明没有同比例改进空间,但仍可能存在个别维度的投入冗余或产出不足,这就是“松弛变量”。这种情况下该决策单元称为“弱有效”,而不是“强有效”。程序如果只输出效率分数,会误判它为前沿单元,但松弛变量暴露了它的短板。
解决:判断有效性的标准应该是效率为1且所有松弛变量为0。在程序里加一个判断逻辑,标记出弱有效单元。实际业务中,这类单元不需要同比例压缩投入,但要单独削减冗余投入、增加不足产出,才能完全达到前沿。报告里要把弱有效单元单独列出来,不能和强有效单元混在一起。
5.4 数据量一大,程序跑得很慢
现象:决策单元有几百个,指标几十个,跑一次要几分钟甚至更久。
原因:程序默认对每个决策单元单独求解一个线性规划,如果还用内置的简版单纯形法,复杂度会极高。另一个原因是约束矩阵没有用稀疏存储,或者用了稀疏但没把矩阵正确传给linprog,导致求解效率下降。
解决:改用linprog并确保约束矩阵是稀疏格式,可以使用sparse函数转换。另外,在linprog的选项里设置算法为对偶单纯形法,通常比默认的内点法快。把循环改成parfor并行,尤其是有多核机器的时候,速度提升接近核心数倍。如果内存不够,可以分批求解,但要注意参考集可能会跨批次依赖,所以分批按年度或按行业分,不要随意切割。
5.5 换了数据集,结果怎么都对不上文献
现象:用教材或论文里的公开数据复现效率值,跟文献里的数值差一大截。
原因:可能是模型类型、导向、指标选择不一致。例如文献用的是产出导向BCC,而你用了投入导向CCR;或者文献已经把部分指标做了正向化处理,你却直接用原始数据;还有可能是数据保留的小数位数不同,导致求解器浮点误差。另一个隐蔽原因是松紧变量的处理方式,有些文献报告的效率是“松紧调整后的”效率,有些是“纯径向”效率,两者数值会不同。
解决:先核对文献的方法部分,确认模型和导向。然后检查数据预处理步骤,比如是否去量纲化。DEA的线性规划本质不受单调正变换影响,但取倒数这类非单调变换会改变结果。确认无误后,尝试把linprog的终止容差设为1e-9,这通常能减少舍入差异。最后再对一遍指标方向,尤其是反向指标的处理方式,这一步最容易埋雷。
5.6 中文路径导致读取失败
现象:在运行xlsread时提示文件不存在,但文件明明在桌面或某个文件夹里。
原因:Matlab在某些系统环境下对全中文路径或含空格的路径支持不够好,尤其是旧版本。如果路径里还有特殊符号,比如括号、感叹号,xlsread会直接失败。
解决:把路径和文件名都改成纯英文,并统一放到工作目录下。读取代码里尽量用相对路径,不要用绝对路径。如果已经混入中文,用dir命令查看实际文件名,先复制到当前目录。更稳妥的做法是把数据文件放在与脚本相同的文件夹,然后用fileparts(mfilename('fullpath'))获取脚本所在路径,再拼接数据文件路径。从那以后,我所有的DEA数据文件都强制命名成纯英文。
6. 进阶:用循环批量计算面板Malmquist指数
最后聊一个很实用的进阶场景:如果你的数据不是单一截面,而是多年的面板数据,怎么在现有程序包的基础上快速算出Malmquist指数。Malmquist指数的本质是跨时期的效率变化,需要用到t期和t+1期的混合距离函数。很多商业软件把这块单独收费,但其实核心代码不复杂:只要会调dea_main,就能循环算出来。
6.1 把多个时期的数据装进三维数组
假设你有三年数据,每年都是同样的决策单元和同样的指标数量。先在Excel里按年度分sheet,然后一次性读入,存成三维数组X_all(n, m, T)和Y_all(n, s, T)。读取代码用一个循环,每个sheet调一次xlsread,然后把数据放入对应页。
% 按年度读取,年份存储在字符串数组中 years = {'2021', '2022', '2023'}; T = length(years); X_all = zeros(n, m, T); Y_all = zeros(n, s, T); for t = 1:T sheetName = years{t}; [num, ~] = xlsread('panel_data.xlsx', sheetName); X_all(:, :, t) = num(:, 1:2); Y_all(:, :, t) = num(:, 3); end这样每个时期的数据都有了。注意读取时确保每个sheet的行顺序一致,否则跨期计算会错位。我一般会在Excel里先把决策单元名称放在第一列,然后读取时把名称单独存一个变量,在循环里检查名称是否一致。这一步很关键,因为面板数据最常见的错误就是某一年少了一个决策单元,导致整个矩阵错位。
6.2 Malmquist指数的分解与输出
Malmquist指数可以分解为技术效率变化(EC)和技术进步(TC)。计算需要四组距离值:t期投入产出在t期前沿的效率、t+1期在t+1期前沿的效率、t期在t+1期前沿的效率、t+1期在t期前沿的效率。前两个可以直接用dea_main算,后两个需要把两期数据合并成一个前沿面,或者交换X和Y的位置。
% 以第t期到第t+1期为例 % eff_t_t 和 eff_t1_t1 是常规期内效率 [eff_t_t, ~, ~, ~] = dea_main(X_all(:, :, t), Y_all(:, :, t), opts); [eff_t1_t1, ~, ~, ~] = dea_main(X_all(:, :, t+1), Y_all(:, :, t+1), opts); % eff_t_t1:用t期的投入产出,但前沿面由t+1期构成 % 常见做法是把t期数据当作被评价对象,t+1期数据作为参考集 [eff_t_t1, ~, ~, ~] = dea_mixed(X_all(:, :, t), Y_all(:, :, t), ... X_all(:, :, t+1), Y_all(:, :, t+1), opts);这里的dea_mixed函数是程序包扩展出来的,专门计算混合距离。如果没有这个函数,你也可以把两期的数据按行拼接成一个新矩阵,然后运行dea_main,再取出对应行的效率值。不过拼接方法要注意,DMU数量扩大后,前端面由两期所有单元共同构成,这样算出来的不是纯粹的跨期前沿,需要小心解读。我一般建议直接扩展一个混合函数,代码量不大。
计算完成后,Malmquist指数 = EC * TC,其中EC = eff_t1_t1 / eff_t_t,TC = sqrt((eff_t_t / eff_t_t1) * (eff_t1_t / eff_t1_t1))。把这几个值都输出到表格里,就能看出每个决策单元的效率变化是来自管理水平改进还是技术边界移动。最后补充一句:跑面板数据前,我会强制检查每一期的有效单元数量和参考集稳定性,如果某一期的前沿面异常,整个指数都会被污染。
这个技巧我已经在好几个项目里用过,每次跑完都习惯性地把四组效率矩阵单独存一份mat文件,方便追溯。从那以后,我每次做面板DEA都会先跑一个单期结果做基准,再上Malmquist循环,一旦发现跨期结果异常,能立刻定位是哪一步出了问题。希望这些经验能帮你的效率分析少走几步弯路,也让你能真正把这套程序用出价值。
本文还有配套的精品资源,点击获取