用 Python 做 DEA 效率评价的同行，应该都有过这种体会：CCR、BCC 这类径向模型处理常规的投入产出数据还算顺手，一旦数据里出现二氧化碳排放、废水、不良贷款这类非期望产出，径向模型就特别别扭。这几年能源经济、绿色金融、产业效率领域用得越来越多的，是基于 NDDF（非径向方向性距离函数）的 Malmquist-Luenberger 指数及其分解，也就是标题里那串关键词的组合。这篇文章我准备把 NDDF 的数学规划模型、ML 指数的分解逻辑、完整的 Python 实现代码，以及我在实际跑数据时踩过的坑一次性写清楚。你要是正在做绿色全要素生产率测算、碳排放绩效评价、银行效率动态分析这类课题，这篇文章可以直接拿来当操作手册用。
1. 为什么要从径向模型转向NDDF
1.1 径向DEA在非期望产出面前的尴尬
先讲一个基础判断：DEA 的核心价值在于不需要事先假设生产函数，而是通过线性规划直接包络出生产前沿。CCR 和 BCC 这两个最经典的模型都是径向模型，意思是所有投入或产出都按同一个比例进行缩放。这种“一个比例打天下”的处理方式，在对只有好投入好产出的传统场景中还能接受，一旦加入非期望产出，问题就暴露得特别明显。
拿一个火电厂举例：假设我们要评价它的环境效率，投入有资本、劳动、煤炭消耗，期望产出是发电量，非期望产出是二氧化碳排放。径向模型在测算低效程度时，会要求资本、劳动、煤炭、发电量、排放全部按照同一个 β 去压缩或扩张。但现实中的减排路径很少是“所有要素同比例齐步走”的，更常见的是煤炭多削减一点、劳动少调整一点、发电量尽量不降甚至还要增长。径向模型对这种异质性完全无能为力。
学界早期处理非期望产出的办法有两种：一种是把排放物直接当成投入来处理，另一种是先做一个单调变换再当作期望产出处理。这两类做法都有明显缺陷。把排放当投入，隐含了“排放越少越好”的线性关系，却破坏了实物型生产过程的因果含义；做数据变换，又会破坏 DEA 模型对数据单调性的要求，结果容易被个别异常值带偏。所以后来方向性距离函数（DDF，Directional 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_k
z_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到ML
DEA 的截面效率测算只能回答“在同一个时间点上，谁更高效”这个问题。但政策研究里更关心的往往是动态变化：今年的生产效率相比去年是提升了还是退步了？提升多少来自管理改善，多少来自技术进步？这类问题需要引入生产率指数。
经典 Malmquist 指数就是干这个的。它通过比较同一决策单元在两个时期的投入产出变化，把全要素生产率（TFP）变动分解为效率变化（EC）和技术变化（TC）。但传统 Malmquist 指数是基于径向距离函数的，同样面临非期望产出难以纳入的问题。于是 Chung、Färe 和 Grosskopf 在 1997 年把方向性距离函数引入 Malmquist 指数的框架，提出了 Malmquist-Luenberger 指数，也就是我们常说的 ML 指数。
更准确地说，ML 指数是对 Malmquist 指数的一种“绿色改造”。它用方向性距离函数替代普通距离函数，使得生产率测算能够同时容纳期望产出的增加和非期望产出的减少。再往后，研究者把非径向的 NDDF 嵌入 ML 指数的公式，就得到了基于 NDDF 的 ML 指数及其分解——这正是标题里那套方法的完整技术路线。
2.2 ML指数的构造公式
要计算第 t 期到第 t+1 期的 ML 指数，需要四个距离函数值，分别对应不同时期技术前沿和不同时期观测值的组合：
- D_t(x_t, y_t, b_t)：以 t 期技术为前沿，测 t 期观测值的方向距离
- D_t(x_{t+1}, y_{t+1}, b_{t+1})：以 t 期技术为前沿，测 t+1 期观测值的方向距离
- D_{t+1}(x_t, y_t, b_t)：以 t+1 期技术为前沿，测 t 期观测值的方向距离
- D_{t+1}(x_{t+1}, y_{t+1}, b_{t+1})：以 t+1 期技术为前沿，测 t+1 期观测值的方向距离
为了写公式方便，先定义 E_t(s) = 1 + D_t(x_s, y_s, b_s)，表示以 t 期技术测度 s 期观测值时，1 加上 NDDF 距离值。那么第 t 期到第 t+1 期的 ML 指数可以写成：
ML = sqrt[ (E_t(t) / E_t(t+1)) × (E_{t+1}(t) / E_{t+1}(t+1)) ]
ML 指数进一步分解为效率变化 EC 和技术变化 TC：
EC = E_t(t) / E_{t+1}(t+1)
TC = sqrt[ (E_{t+1}(t+1) / E_t(t+1)) × (E_{t+1}(t) / E_t(t)) ]
需要注意的是，这里可能出现 E 值大于 1、等于 1、小于 1 三种情况。E 等于 1 意味着该观测值正好落在前沿上；E 大于 1 说明该观测值在前沿内部（低效）；E 小于 1 说明该观测值超越了当前前沿。允许 E 小于 1 这一点非常重要，后面讲代码的时候我会再强调。
2.3 分解结果的经济解读
拿到 EC、TC、ML 三个值以后，怎么判断趋势？规律很简单：
- ML 大于 1：绿色全要素生产率从第 t 期到第 t+1 期是提升的
- 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, 2019 |
| dmu | 决策单元名称或编号 | 北京、河北、上海 |
| K | 资本存量 | 亿元 |
| L | 劳动力 | 万人 |
| E | 能源消耗 | 万吨标准煤 |
| Y | 期望产出 | 亿元 |
| B | 非期望产出 | 万吨 CO2 |
变量名用英文短名是因为 Python 代码里操作方便，展示结果时再映射成中文即可。数据量纲不影响线性规划求解结果，因为 DEA 模型本质上是按比例计算的，但量纲差异过大会影响求解器数值稳定性，所以如果各列数量级差别特别大，比如 GDP 到千亿、劳动力只有几百，可以考虑先做无量纲化处理，不过要说明的是这不会改变 DEA 距离值。
我写代码时习惯把所有列名做成常量，避免字符串笔误。特别是 ML 指数要处理四个距离函数，稍不留神引用错列名，结果就全错了。
3.2 用scipy.optimize.linprog求解NDDF
Python 生态里求解线性规划的常用工具是 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_x=None, w_y=1/3, w_b=1/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_ub=np.array(A_ub), b_ub=np.array(b_ub), bounds=bounds, method='highs') 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年到第t+1年的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_{t+1}(x_{t+1}, y_{t+1}, b_{t+1}) d_t_cross1 = solve_nddf(ref_t, obj_tn) # D_t(x_{t+1}, y_{t+1}, b_{t+1}) d_tn_cross2 = solve_nddf(ref_tn, obj_t) # D_{t+1}(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}-{t+1}', '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_index=True) print(final_result.head())以两个相邻年份 2018-2019 的输出为例，结果表格大致长这样：
| dmu | period | ML | EC | TC | D_tt | D_tn_tn | D_t_cross1 | D_tn_cross2 |
|---|---|---|---|---|---|---|---|---|
| A地区 | 2018-2019 | 1.052 | 1.018 | 1.033 | 0.024 | 0.006 | 0.131 | 0.118 |
| B地区 | 2018-2019 | 0.982 | 0.995 | 0.987 | 0.003 | 0.008 | 0.035 | 0.051 |
| C地区 | 2018-2019 | 1.137 | 1.074 | 1.058 | 0.047 | -0.027 | 0.159 | 0.097 |
从这份结果可以快速读取信息：A 地区的 ML 大于 1，其中 EC 和 TC 都大于 1，说明效率和前沿扩张同时贡献了增长；B 地区的 ML 小于 1，EC 和 TC 都小于 1，整体呈现轻微退步；C 地区 ML 改善明显，EC 是主要推手，说明 C 地区在向前沿追赶，但技术进步幅度相对小一些。另外注意 C 地区 D_tn_tn 出现了负值，这说明它在 t+1 期就已经跑到了当年参考前沿的外部，属于超高效单元，这种情况在允许 β 取负值的设定下是正常的。
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_{t+1}, y_{t+1}, b_{t+1}) 这种交叉距离函数时，频繁出现求解失败。
原因在前面提过：当归属 DMU 在第 t+1 期已经优于第 t 期前沿时，它相对第 t 期前沿的距离应该表现为负值。这个负值不是错误，而是“超效率”的体现。如果强行把 β 限制成非负，模型找不到可行解，返回值就只能是 NaN。
所以我的建议是：模型公式里写 β ≥ 0 是为了推导方便，实际编码时把 β 的 bounds 设为 (None, None)，让求解器自由处理。这一点对 ML 指数的计算尤其重要，因为四个距离函数里有交叉项，交叉项出现负值的概率相当高。
4.2 参考技术集合不能搞混
ML 指数看起来只是四个距离函数的组合，但代码里最容易出错的位置就是“谁是参考技术，谁是被评价对象”。D_t(x_t, y_t, b_t) 和 D_t(x_{t+1}, y_{t+1}, b_{t+1}) 的区别仅仅在于参考年份相同、被评价对象不同；D_t(x_t, y_t, b_t) 和 D_{t+1}(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 为负数 | 跨期距离函数组合出现异常比值 | 检查四个距离值是否满足基本逻辑关系，比如交叉项是否拿到 NaN |
4.5 关于负距离值的一点心得
我刚开始跑跨期距离函数时，看到 D 出现负值总觉得是不是写错了，后来反复对照文献才发现，负值恰恰说明该单元相比某个时期的前沿“更高效”。这种超效率现象在环境效率评价里很常见，比如某个地区在 t+1 期大幅淘汰了落后产能、提高了清洁能源比例，它相对于 t 期前沿自然就是“领先”的。
但在论文里报告结果时，如果出现负值，一定不要顺手就删掉或替换成 0。正确的做法是保留原值，并在方法部分写清楚：NDDF 允许距离值取负，负值表示观测点位于前沿外侧。否则审稿人可能找出计算逻辑漏洞，那才叫得不偿失。
根据自己的实操经验，最后再分享两个小技巧。第一个是写代码时把 solve_nddf 单独封装，不要把所有逻辑都堆在主循环里。这看起来只是一个代码风格问题，但当你需要换一组权重重新跑结果、或者扩展成多人份面板时，单独封装的函数改动成本极低。第二个是运算完成后，把 D_tt、D_tn_tn 这几个中间距离值一起输出到结果表里，然后用 Excel 或 pandas 做一次手工抽查，随机挑两三个 DMU，按公式逐步手算 ML 和 EC，确认和代码输出对得上。这套代码的核心逻辑看起来不复杂，但四个距离函数之间的交叉引用非常容易在迭代中改出偏差，抽检是最笨也最可靠的方法。
