用 Python 做 DEA 效率评价的同行应该都有过这种体会CCR、BCC 这类径向模型处理常规的投入产出数据还算顺手一旦数据里出现二氧化碳排放、废水、不良贷款这类非期望产出径向模型就特别别扭。这几年能源经济、绿色金融、产业效率领域用得越来越多的是基于 NDDF非径向方向性距离函数的 Malmquist-Luenberger 指数及其分解也就是标题里那串关键词的组合。这篇文章我准备把 NDDF 的数学规划模型、ML 指数的分解逻辑、完整的 Python 实现代码以及我在实际跑数据时踩过的坑一次性写清楚。你要是正在做绿色全要素生产率测算、碳排放绩效评价、银行效率动态分析这类课题这篇文章可以直接拿来当操作手册用。1. 为什么要从径向模型转向NDDF1.1 径向DEA在非期望产出面前的尴尬先讲一个基础判断DEA 的核心价值在于不需要事先假设生产函数而是通过线性规划直接包络出生产前沿。CCR 和 BCC 这两个最经典的模型都是径向模型意思是所有投入或产出都按同一个比例进行缩放。这种“一个比例打天下”的处理方式在对只有好投入好产出的传统场景中还能接受一旦加入非期望产出问题就暴露得特别明显。拿一个火电厂举例假设我们要评价它的环境效率投入有资本、劳动、煤炭消耗期望产出是发电量非期望产出是二氧化碳排放。径向模型在测算低效程度时会要求资本、劳动、煤炭、发电量、排放全部按照同一个 β 去压缩或扩张。但现实中的减排路径很少是“所有要素同比例齐步走”的更常见的是煤炭多削减一点、劳动少调整一点、发电量尽量不降甚至还要增长。径向模型对这种异质性完全无能为力。学界早期处理非期望产出的办法有两种一种是把排放物直接当成投入来处理另一种是先做一个单调变换再当作期望产出处理。这两类做法都有明显缺陷。把排放当投入隐含了“排放越少越好”的线性关系却破坏了实物型生产过程的因果含义做数据变换又会破坏 DEA 模型对数据单调性的要求结果容易被个别异常值带偏。所以后来方向性距离函数DDFDirectional Distance Function被引入进来把投入和产出的调整方向显式地写进模型里。方向性距离函数确实解决了“非期望产出往哪个方向调”的问题它通过一个预设的方向向量 g 告诉模型投入往下降、期望产出往上升、非期望产出往下降。但 DDF 本质上还是径向的所有变量共享同一个调整比例 β只是在不同的方向上调整。于是 NDDF 应运而生。1.2 NDDF的数学表达NDDF 的全称是 Non-radial Directional Distance Function中文常译作非径向方向性距离函数。它的核心思想说穿了很简单如果径向模型是“一根绳上的蚂蚱所有变量一起动”那 NDDF 就是“各管各的每个变量有自己独立的调整系数”。我们用数学语言把它写清楚。假设每个决策单元有投入向量 x (K, L, E)其中 K 是资本存量L 是劳动力E 是能源消耗期望产出是 Y非期望产出是 B。对第 k 个决策单元构造一个权重向量w (w_K, w_L, w_E, w_Y, w_B)方向向量取g (−x_k, y_k, −b_k)也就是投入到负方向缩减、期望产出到正方向扩张、非期望产出到负方向压缩。NDDF 的值通过求解下面这个线性规划得到max w_K·β_K w_L·β_L w_E·β_E w_Y·β_Y w_B·β_B约束条件Σ z_i·K_i ≤ (1 − β_K)·K_kΣ z_i·L_i ≤ (1 − β_L)·L_kΣ z_i·E_i ≤ (1 − β_E)·E_kΣ z_i·Y_i ≥ (1 β_Y)·Y_kΣ z_i·B_i ≤ (1 − β_B)·B_kz_i ≥ 0β_K, β_L, β_E, β_Y, β_B ≥ 0模型的目标函数是带权重的 β 之和。z_i 是结构变量用于把被评价单元投射到由所有样本单元张成的生产可能性集合上。β 越大说明被评价单元距离生产前沿越远也就是越低效。这里 I 用“三个投入 一个期望产出 一个非期望产出”的设定作为例子实际操作上完全可以扩展成任意数量的投入和产出只需要按同样的规则往目标函数和约束矩阵里追加变量就行。1.3 权重向量怎么定才合理NDDF 和径向模型最大的区别之一是引入了外生权重向量 w。这意味着不同的权重设定会直接影响距离函数值的大小进而影响后续 ML 指数和分解结果。这个点很多初学者容易忽略。目前文献里最常见的设定是“投入组总权重占 1/3、期望产出占 1/3、非期望产出占 1/3”。在这个框架下如果投入有 K、L、E 三个变量那么每个投入的权重就是 1/9这样三者合计正好 1/3。即w (1/9, 1/9, 1/9, 1/3, 1/3)这样设定的好处是三类变量——投入、期望产出、非期望产出——在目标函数里具有平等的“话语权”。如果你研究的问题是碳排放约束下的经济增长可以适当提高非期望产出的权重体现政策上对减排的更高要求反过来如果关注的重点是就业稳定那就把劳动力的权重抬高。需要提醒的是权重本身是一种价值判断不是模型内生出来的。所以在学术论文里如果你用了 NDDF一定要说明权重设定的依据并且最好做一个敏感性分析看看在几组不同权重下结论是否稳定。1.4 径向DEA、DDF、NDDF三者对比模型调整方式是否区分变量能否直接处理非期望产出权重是否外生典型应用CCR/BCC径向所有投入或产出同比例缩放否较困难需变换处理否传统效率评价DDF径向方向性距离所有变量按同一 β 沿方向向量调整否可以否环境效率评价NDDF非径向每个变量有独立 β是可以是绿色全要素生产率、碳排放绩效、ML指数从这张表能看出NDDF 相比前两者最大的进步在于“区分变量”和“引入权重”这两件事。区分变量让模型结果更贴近实际生产行为引入权重让研究者可以把政策偏好或价值判断嵌入到模型里。这也是它在 ML 指数测算场景下成为主流选择的原因。2. Malmquist-Luenberger指数给绿色生产率加上“时间维”2.1 从Malmquist到MLDEA 的截面效率测算只能回答“在同一个时间点上谁更高效”这个问题。但政策研究里更关心的往往是动态变化今年的生产效率相比去年是提升了还是退步了提升多少来自管理改善多少来自技术进步这类问题需要引入生产率指数。经典 Malmquist 指数就是干这个的。它通过比较同一决策单元在两个时期的投入产出变化把全要素生产率TFP变动分解为效率变化EC和技术变化TC。但传统 Malmquist 指数是基于径向距离函数的同样面临非期望产出难以纳入的问题。于是 Chung、Färe 和 Grosskopf 在 1997 年把方向性距离函数引入 Malmquist 指数的框架提出了 Malmquist-Luenberger 指数也就是我们常说的 ML 指数。更准确地说ML 指数是对 Malmquist 指数的一种“绿色改造”。它用方向性距离函数替代普通距离函数使得生产率测算能够同时容纳期望产出的增加和非期望产出的减少。再往后研究者把非径向的 NDDF 嵌入 ML 指数的公式就得到了基于 NDDF 的 ML 指数及其分解——这正是标题里那套方法的完整技术路线。2.2 ML指数的构造公式要计算第 t 期到第 t1 期的 ML 指数需要四个距离函数值分别对应不同时期技术前沿和不同时期观测值的组合D_t(x_t, y_t, b_t)以 t 期技术为前沿测 t 期观测值的方向距离D_t(x_{t1}, y_{t1}, b_{t1})以 t 期技术为前沿测 t1 期观测值的方向距离D_{t1}(x_t, y_t, b_t)以 t1 期技术为前沿测 t 期观测值的方向距离D_{t1}(x_{t1}, y_{t1}, b_{t1})以 t1 期技术为前沿测 t1 期观测值的方向距离为了写公式方便先定义 E_t(s) 1 D_t(x_s, y_s, b_s)表示以 t 期技术测度 s 期观测值时1 加上 NDDF 距离值。那么第 t 期到第 t1 期的 ML 指数可以写成ML sqrt[ (E_t(t) / E_t(t1)) × (E_{t1}(t) / E_{t1}(t1)) ]ML 指数进一步分解为效率变化 EC 和技术变化 TCEC E_t(t) / E_{t1}(t1)TC sqrt[ (E_{t1}(t1) / E_t(t1)) × (E_{t1}(t) / E_t(t)) ]需要注意的是这里可能出现 E 值大于 1、等于 1、小于 1 三种情况。E 等于 1 意味着该观测值正好落在前沿上E 大于 1 说明该观测值在前沿内部低效E 小于 1 说明该观测值超越了当前前沿。允许 E 小于 1 这一点非常重要后面讲代码的时候我会再强调。2.3 分解结果的经济解读拿到 EC、TC、ML 三个值以后怎么判断趋势规律很简单ML 大于 1绿色全要素生产率从第 t 期到第 t1 期是提升的ML 小于 1绿色全要素生产率呈下降趋势EC 大于 1被评价单元在向生产前沿移动也就是“追赶”效应说明效率在改善EC 小于 1该单元在远离生产前沿效率出现退化TC 大于 1生产前沿整体向外扩张技术进步发生了TC 小于 1前沿向内收缩通常被解读为技术退步实际研究中EC 和 TC 经常出现方向不一致。比如某地区 ML 大于 1主要贡献来自 TC 大于 1说明该地区生产率提升靠的是新技术、新工艺的引入如果 EC 很大但 TC 很小甚至小于 1说明该地区是在追赶现有前沿但行业整体技术没有进步长期增长动力不足。这种区分对政策制定特别有价值。2.4 为什么非径向版本更稳很多人会问既然已经有基于 DDF 的 ML 指数为什么还要用 NDDF 版本我实际跑数据后的体会是径向 DDF 的 ML 在分解时容易出现 EC 和 TC 被“共用一个 β”过度捆绑的问题。因为径向 DDF 只允许所有变量同步缩放模型无法区分“减排促成的改进”和“产出扩张带来的改进”于是技术变化和效率变化的区分对权重结构特别敏感。NDDF 版本给每个变量独立的调整系数后整体效率损失可以被更精细地归因非期望产出的高 β 说明减排拖累了效率期望产出的低 β 说明产出扩张空间有限。这种精细度让 ML 的分解结果更稳定也更经得起同行评议。如果你想在论文里使用这个方法审稿人问到“为什么选 NDDF 而不是传统 DDF 或 SBM 类模型”答案其实就是这四个字区分变量。3. Python实现从面板数据到ML指数3.1 数据准备与结构设计要跑通 NDDF 框架下的 ML 指数第一步是把数据整理成面板结构。我建议使用 pandas 的 DataFrame每一行代表一个决策单元在某一年份的观测值。列名建议这样设置列名含义示例year年份2018, 2019dmu决策单元名称或编号北京、河北、上海K资本存量亿元L劳动力万人E能源消耗万吨标准煤Y期望产出亿元B非期望产出万吨 CO2变量名用英文短名是因为 Python 代码里操作方便展示结果时再映射成中文即可。数据量纲不影响线性规划求解结果因为 DEA 模型本质上是按比例计算的但量纲差异过大会影响求解器数值稳定性所以如果各列数量级差别特别大比如 GDP 到千亿、劳动力只有几百可以考虑先做无量纲化处理不过要说明的是这不会改变 DEA 距离值。我写代码时习惯把所有列名做成常量避免字符串笔误。特别是 ML 指数要处理四个距离函数稍不留神引用错列名结果就全错了。3.2 用scipy.optimize.linprog求解NDDFPython 生态里求解线性规划的常用工具是 scipy.optimize.linprog。它默认求解的是最小化问题min c^T x s.t. A_ub x ≤ b_ub A_eq x b_eq bounds ≤ x ≤ bounds而 NDDF 的目标是最大化加权 β 之和所以要做两件事第一把目标向量 c 取负号把最大化变成最小化 第二把期望产出的不等式约束从“≥”变成“≤”方法是两边同时乘 -1。假设我们设定三个投入 K、L、E一个期望产出 Y一个非期望产出 B那么变量向量可以按顺序拼接为[β_K, β_L, β_E, β_Y, β_B, z_1, z_2, ..., z_N]其中前五个是缩放因子 β后面的 z 是结构变量N 是参考技术集合里的决策单元数量。约束矩阵的系数要按这个顺序逐列填充。下面给出一个完整的求解函数这段代码可以直接用到你自己的项目里。import numpy as np import pandas as pd from scipy.optimize import linprog # 列名常量避免拼写错误 X_COLS [K, L, E] Y_COL Y B_COL B def solve_nddf(ref_df, obj_row, w_xNone, w_y1/3, w_b1/3): 求解单个DMU的NDDF距离值。 参数 ---------- ref_df : DataFrame 参考技术集合通常某一年所有DMU的投入产出数据。 obj_row : Series 被评价DMU的观测值。 w_x : list 各投入变量的权重如果不给默认投入组总权重1/3组内均分。 w_y : float 期望产出权重。 w_b : float 非期望产出权重。 返回 ------- d : float NDDF的目标函数值即加权β之和。 n len(ref_df) # 参考集里DMU的数量 n_x len(X_COLS) # 投入变量个数 if w_x is None: w_x [1 / (3 * n_x) for _ in range(n_x)] # 变量顺序: [beta_x1..beta_xn, beta_y, beta_b, z_0..z_{n-1}] n_beta n_x 2 n_vars n_beta n # 目标函数取负因为linprog默认最小化 c [-w for w in w_x] [-w_y, -w_b] [0] * n A_ub [] b_ub [] # 投入约束: sum z_i * X_i beta_x * obj_x obj_x for j, col in enumerate(X_COLS): row [0.0] * n_vars row[j] obj_row[col] # beta_x 的系数 for i in range(n): row[n_beta i] ref_df.iloc[i][col] # z_i 的系数 A_ub.append(row) b_ub.append(obj_row[col]) # 期望产出约束: sum z_i * Y_i - beta_y * obj_y obj_y # 转换为 -sum z_i * Y_i beta_y * obj_y -obj_y row_y [0.0] * n_vars row_y[n_x] obj_row[Y_COL] # beta_y 的系数 for i in range(n): row_y[n_beta i] -ref_df.iloc[i][Y_COL] A_ub.append(row_y) b_ub.append(-obj_row[Y_COL]) # 非期望产出约束: sum z_i * B_i beta_b * obj_b obj_b row_b [0.0] * n_vars row_b[n_x 1] obj_row[B_COL] # beta_b 的系数 for i in range(n): row_b[n_beta i] ref_df.iloc[i][B_COL] A_ub.append(row_b) b_ub.append(obj_row[B_COL]) # 变量边界: beta 允许取负值z 必须非负 bounds [(None, None)] * n_beta [(0, None)] * n result linprog(c, A_ubnp.array(A_ub), b_ubnp.array(b_ub), boundsbounds, methodhighs) if not result.success: return np.nan return -result.fun这段代码有几个关键细节需要重点说明。第一个细节是 β 变量的边界设置。我用了 (None, None)也就是允许 β 取负值。很多初次写 NDDF 代码的人会下意识写成 (0, None)但这样会造成跨期距离函数大量无解。原因在于当被评价对象在参考技术前沿的外侧也就是比前沿更高效时NDDF 的距离值本来就是负的——它表示该 DMU 在方向向量的反方向还有扩张空间。如果你把 β 限制为非负这种“超级效率”情况就直接被模型排除了求解器会报无界或无可行解。第二个细节是期望产出的约束处理。原始约束是“Σ z_i·Y_i ≥ (1 β_Y)·Y_k”但 linprog 只接受小于等于不等式所以我在代码里把左右两边同时乘了 -1变成“−Σ z_i·Y_i β_Y·Y_k ≤ −Y_k”。这个转换必须在约束矩阵的行里正确体现否则结果偏差极大。第三个细节是方向向量的处理。我在代码里没有单独生成一个方向向量对象因为 NDDF 模型的“方向”是隐含在约束右侧的观测值和 β 前面的系数里的。投入和非期望产出用 obj_row 的原始值作为 β 的系数且约束是小于等于这就实现了“沿着负方向缩减”期望产出同样用 obj_row 的原始值但经过乘 -1 转换后变成“沿正方向扩张”的效果。这一套隐式写法比显式拼方向向量更不易出错。3.3 计算四个距离并生成ML指数有了 solve_nddf 这个核心函数接下来只需要按相邻年份循环对每个 DMU 求解四个距离函数。这里我把计算过程写成一个独立的 compute_ml_for_year_gap 函数它接收完整面板数据、当前年份 t返回该年度所有 DMU 的 ML、EC、TC 值。def compute_ml_for_year_gap(data, t): 计算从第t年到第t1年的ML指数及其分解。 results [] # 当年和下一年的全部数据作为参考技术集合 ref_t data[data[year] t].copy() ref_tn data[data[year] t 1].copy() # 只对两个年份都出现的DMU计算 dmus set(ref_t[dmu]) set(ref_tn[dmu]) for dmu in dmus: obj_t ref_t[ref_t[dmu] dmu].iloc[0] obj_tn ref_tn[ref_tn[dmu] dmu].iloc[0] # 四个距离函数 d_tt solve_nddf(ref_t, obj_t) # D_t(x_t, y_t, b_t) d_tn_tn solve_nddf(ref_tn, obj_tn) # D_{t1}(x_{t1}, y_{t1}, b_{t1}) d_t_cross1 solve_nddf(ref_t, obj_tn) # D_t(x_{t1}, y_{t1}, b_{t1}) d_tn_cross2 solve_nddf(ref_tn, obj_t) # D_{t1}(x_t, y_t, b_t) e_tt 1 d_tt e_tn_tn 1 d_tn_tn e_t_cross1 1 d_t_cross1 e_tn_cross2 1 d_tn_cross2 ml np.sqrt((e_tt / e_t_cross1) * (e_tn_cross2 / e_tn_tn)) ec e_tt / e_tn_tn tc ml / ec results.append({ dmu: dmu, period: f{t}-{t1}, ML: ml, EC: ec, TC: tc, D_tt: d_tt, D_tn_tn: d_tn_tn, D_t_cross1: d_t_cross1, D_tn_cross2: d_tn_cross2 }) return pd.DataFrame(results)调用方式很简单如果你有 2015 到 2020 年的数据只需要写一个循环对 2015、2016、2017、2018、2019 各调用一次这个函数再用 pd.concat 拼接所有年份的结果即可。data pd.read_csv(panel_data.csv) all_results [] for t in range(2015, 2020): res compute_ml_for_year_gap(data, t) all_results.append(res) final_result pd.concat(all_results, ignore_indexTrue) print(final_result.head())以两个相邻年份 2018-2019 的输出为例结果表格大致长这样dmuperiodMLECTCD_ttD_tn_tnD_t_cross1D_tn_cross2A地区2018-20191.0521.0181.0330.0240.0060.1310.118B地区2018-20190.9820.9950.9870.0030.0080.0350.051C地区2018-20191.1371.0741.0580.047-0.0270.1590.097从这份结果可以快速读取信息A 地区的 ML 大于 1其中 EC 和 TC 都大于 1说明效率和前沿扩张同时贡献了增长B 地区的 ML 小于 1EC 和 TC 都小于 1整体呈现轻微退步C 地区 ML 改善明显EC 是主要推手说明 C 地区在向前沿追赶但技术进步幅度相对小一些。另外注意 C 地区 D_tn_tn 出现了负值这说明它在 t1 期就已经跑到了当年参考前沿的外部属于超高效单元这种情况在允许 β 取负值的设定下是正常的。3.4 运行环境与性能提示这套代码只需要 pandas、numpy、scipy 三个基础库Python 3.8 以上版本都能跑。如果你用的是 VS Code给 Python 配置好解释器然后在终端里执行 pip install pandas numpy scipy脚本就能直接跑。数据量在几十个 DMU、五六年面板规模下循环加 linprog 的耗时完全可接受如果你有几百个 DMU建议先跑一个小规模测试确认结果稳定再扩大范围否则逐循环求解时间会比较长。后续如果想提速有两个方向一是把同一年份所有 DMU 的 NDDF 求解改成矩阵化批量求解但代码复杂度明显上升二是换成商业求解器如 Gurobi在小规模数据上提升有限不划算。我的建议是先在当前框架下跑通确认方法论没有问题再做性能优化。4. 实操中容易踩的坑4.1 β变量千万别限制成非负这是我在写 NDDF 代码时踩过最深的一个坑。很多资料在介绍 NDDF 模型时β 后面都写着“≥ 0”初学者照着抄到代码里就会把 bounds 写成 (0, None)。结果就是截面测算还行一到跨期求 D_t(x_{t1}, y_{t1}, b_{t1}) 这种交叉距离函数时频繁出现求解失败。原因在前面提过当归属 DMU 在第 t1 期已经优于第 t 期前沿时它相对第 t 期前沿的距离应该表现为负值。这个负值不是错误而是“超效率”的体现。如果强行把 β 限制成非负模型找不到可行解返回值就只能是 NaN。所以我的建议是模型公式里写 β ≥ 0 是为了推导方便实际编码时把 β 的 bounds 设为 (None, None)让求解器自由处理。这一点对 ML 指数的计算尤其重要因为四个距离函数里有交叉项交叉项出现负值的概率相当高。4.2 参考技术集合不能搞混ML 指数看起来只是四个距离函数的组合但代码里最容易出错的位置就是“谁是参考技术谁是被评价对象”。D_t(x_t, y_t, b_t) 和 D_t(x_{t1}, y_{t1}, b_{t1}) 的区别仅仅在于参考年份相同、被评价对象不同D_t(x_t, y_t, b_t) 和 D_{t1}(x_t, y_t, b_t) 则是被评价对象相同、参考年份不同。任何一个参数传错最终 ML 结果就完全错位。我在代码里刻意用 ref_t 和 ref_tn 来命名参考技术集合用 obj_t 和 obj_tn 来命名被评价对象就是为了降低混淆概率。建议你在自己的代码里也坚持这套命名规则不要精简成 data1、data2 这种无意义的变量名。4.3 相邻年份DMU缺失怎么办面板数据经常存在缺失。比如某一年某个地区的数据没公布或者样本企业在政策合并后消失。如果你的数据是强平衡面板即每个 DMU 在所有年份都存在那没问题但如果存在缺失compute_ml_for_year_gap 里 set(ref_t[dmu]) set(ref_tn[dmu]) 这个交集操作会自动跳过缺失的 DMU。这带来两个后果第一缺失 DMU 的 ML 指数根本算不出来论文里要交代清楚样本口径第二参考技术集合的规模随年份变化有的年份参考集有 30 个 DMU有的年份可能只剩 28 个这会让技术前沿的构建基础不一致跨期比较时需要注意。一般情况下只要缺失量不大同行也能接受如果缺失严重建议用内插或匹配方法补齐再往下走。4.4 结果异常排查速查表最后整理一份我自己用的排查速查表当你发现 ML 结果异常时按表里的顺序逐项检查大多数问题都能定位。问题现象可能原因处理方式所有 ML、EC、TC 均为 NaNβ 边界被限制为非负将 beta bounds 改为 (None, None)所有 DMU 的 ML 恒等于 1参考年份和被评价年份在代码里被写成了同一时期检查四个距离函数的参数传递大量 ML 值大于 1.5 或小于 0.5权重向量设置不当某个变量权重过高或过低重新按 1/3-1/3-1/3 原则设置权重某几个 DMU 的 TC 出现极端值该 DMU 跨期出现超高效负距离值放大了分母核对原始数据是否录入错误关注数据质量结果在不同权重下波动极大样本量太小或投入产出指标选择不当增加样本量或做敏感性分析说明稳健性运行时间过长DMU 数量大循环求解次数多先跑小样本验证再考虑并行或商业求解器某些行出现 ML 为负数跨期距离函数组合出现异常比值检查四个距离值是否满足基本逻辑关系比如交叉项是否拿到 NaN4.5 关于负距离值的一点心得我刚开始跑跨期距离函数时看到 D 出现负值总觉得是不是写错了后来反复对照文献才发现负值恰恰说明该单元相比某个时期的前沿“更高效”。这种超效率现象在环境效率评价里很常见比如某个地区在 t1 期大幅淘汰了落后产能、提高了清洁能源比例它相对于 t 期前沿自然就是“领先”的。但在论文里报告结果时如果出现负值一定不要顺手就删掉或替换成 0。正确的做法是保留原值并在方法部分写清楚NDDF 允许距离值取负负值表示观测点位于前沿外侧。否则审稿人可能找出计算逻辑漏洞那才叫得不偿失。根据自己的实操经验最后再分享两个小技巧。第一个是写代码时把 solve_nddf 单独封装不要把所有逻辑都堆在主循环里。这看起来只是一个代码风格问题但当你需要换一组权重重新跑结果、或者扩展成多人份面板时单独封装的函数改动成本极低。第二个是运算完成后把 D_tt、D_tn_tn 这几个中间距离值一起输出到结果表里然后用 Excel 或 pandas 做一次手工抽查随机挑两三个 DMU按公式逐步手算 ML 和 EC确认和代码输出对得上。这套代码的核心逻辑看起来不复杂但四个距离函数之间的交叉引用非常容易在迭代中改出偏差抽检是最笨也最可靠的方法。