5分钟实战：用Python的SciPy库处理含非期望产出的超效率SBM模型
当你的研究数据中同时包含常规产出和污染排放等非期望产出时，传统DEA模型往往难以给出准确效率评估。本文将带你用Python的SciPy库快速实现一个考虑非期望产出的超效率SBM模型，从数据准备到结果解读一气呵成。
1. 环境准备与数据预处理
在开始建模前，我们需要确保环境配置正确并准备好标准化的数据格式。这个步骤往往被许多教程忽略，但却是保证后续分析可靠性的关键。
首先安装必要的库（如果尚未安装）：
pip install numpy scipy pandas典型的投入产出数据应包含三类指标：
- 投入指标：如劳动力、资本、能源消耗等
- 期望产出：如GDP、产品产量等
- 非期望产出：如CO2排放、废水排放等
建议使用pandas的DataFrame来组织数据，以下是一个标准化的数据结构示例：
import pandas as pd data = { 'DMU': ['A', 'B', 'C', 'D'], 'Labor': [10, 15, 12, 18], # 投入1：劳动力 'Capital': [200, 300, 250, 400], # 投入2：资本 'GDP': [500, 600, 550, 700], # 期望产出 'CO2': [80, 120, 100, 150] # 非期望产出 } df = pd.DataFrame(data).set_index('DMU')关键预处理步骤：
- 数据标准化（特别是当不同指标量纲差异大时）
- 非期望产出指标取负值（因其与效率负相关）
- 检查数据完整性，处理缺失值
注意：非期望产出指标在建模时需要特殊处理，通常采用方向距离函数或将其视为"负产出"。
2. 超效率SBM模型核心实现
超效率SBM模型的核心是求解一系列线性规划问题。我们利用SciPy的linprog函数来实现这一过程，相比传统DEA模型，它能处理松弛变量并允许效率值超过1。
以下是完整的模型实现代码：
import numpy as np from scipy.optimize import linprog def super_sbm(inputs, good_outputs, bad_outputs, crs=True): """ 考虑非期望产出的超效率SBM模型 :param inputs: 投入指标矩阵 (n_dmu × n_inputs) :param good_outputs: 期望产出矩阵 (n_dmu × n_good_outputs) :param bad_outputs: 非期望产出矩阵 (n_dmu × n_bad_outputs) :param crs: 规模报酬不变假设 (True/False) :return: 各DMU的效率值 """ n_dmu = inputs.shape[0] efficiencies = np.zeros(n_dmu) # 将非期望产出视为负产出 outputs = np.hstack([good_outputs, -bad_outputs]) for k in range(n_dmu): # 目标函数系数 c = np.zeros(n_dmu + 1) c[-1] = 1 # 效率值的倒数 # 约束条件 A_ub = np.vstack([ np.hstack([outputs, -outputs[k].reshape(1, -1)]), # 产出约束 np.hstack([-inputs, inputs[k].reshape(1, -1)]), # 投入约束 ]) b_ub = np.zeros(A_ub.shape[0]) if crs: # 规模报酬不变约束 A_ub = np.vstack([A_ub, np.hstack([np.ones(n_dmu), 0])]) b_ub = np.hstack([b_ub, 1]) # 变量边界 bounds = [(0, None)] * n_dmu + [(None, None)] # 求解线性规划 res = linprog(c, A_ub=A_ub, b_ub=b_ub, bounds=bounds, method='highs') if res.success: efficiencies[k] = 1 / res.x[-1] else: efficiencies[k] = np.nan return efficiencies关键参数说明：
crs：决定是否采用规模报酬不变假设method='highs'：使用高性能线性规划求解器bounds：确保权重非负，但允许效率值超过1
3. 实战案例：区域碳排放效率评估
让我们用一个实际案例演示如何使用上述函数。假设我们有四个地区的投入产出数据：
# 投入指标：劳动力(万人), 资本(亿元), 能源(万吨标准煤) inputs = np.array([ [10, 200, 50], [15, 300, 70], [12, 250, 60], [18, 400, 90] ]) # 期望产出：GDP(亿元) good_outputs = np.array([ [500], [600], [550], [700] ]) # 非期望产出：CO2排放(万吨) bad_outputs = np.array([ [80], [120], [100], [150] ]) # 计算效率 eff_crs = super_sbm(inputs, good_outputs, bad_outputs, crs=True) eff_vrs = super_sbm(inputs, good_outputs, bad_outputs, crs=False) print("CRS假设下的效率值:", eff_crs) print("VRS假设下的效率值:", eff_vrs)典型输出结果可能如下：
CRS假设下的效率值: [1.25 0.85 1.10 0.75] VRS假设下的效率值: [1.30 0.90 1.15 0.80]结果解读要点：
- 效率值>1表示该DMU处于生产前沿面上方（超高效）
- 效率值<1表示存在改进空间
- CRS与VRS结果的差异反映规模效率影响
4. 高级应用与常见问题处理
在实际应用中，我们经常会遇到各种特殊情况。以下是几个常见问题的解决方案：
4.1 处理零值或负值数据
当数据中存在零或负值时，传统DEA模型可能失效。解决方法包括：
- 对数据进行平移变换
- 使用特殊的距离函数
- 采用非径向的SBM模型
# 数据平移示例（当存在负值时） shift = np.abs(np.min(inputs)) + 1 shifted_inputs = inputs + shift4.2 模型选择指南
不同场景下的模型选择策略：
| 场景特征 | 推荐模型 | 优势 |
|---|---|---|
| 有非期望产出 | 非期望产出SBM | 准确处理负向指标 |
| 需要区分高效DMU | 超效率SBM | 允许效率值>1 |
| 规模效应显著 | VRS模型 | 考虑规模报酬可变 |
| 面板数据 | 窗口分析或Malmquist | 分析效率动态变化 |
4.3 结果可视化技巧
使用matplotlib可以直观展示效率分布：
import matplotlib.pyplot as plt plt.figure(figsize=(10, 5)) plt.bar(range(len(eff_crs)), eff_crs, alpha=0.7, label='CRS') plt.bar(range(len(eff_vrs)), eff_vrs, alpha=0.7, label='VRS') plt.axhline(1, color='red', linestyle='--') plt.xlabel('DMU') plt.ylabel('Efficiency Score') plt.title('Efficiency Comparison (CRS vs VRS)') plt.legend() plt.show()4.4 性能优化建议
当处理大量DMU时，可采取以下优化措施：
- 使用稀疏矩阵存储约束条件
- 并行化计算各DMU的效率
- 采用更高效的求解器如Gurobi或CPLEX
from multiprocessing import Pool def parallel_sbm(args): """并行计算单个DMU的效率""" k, inputs, outputs = args # ... (省略具体实现) return efficiency # 使用多进程并行计算 with Pool() as p: results = p.map(parallel_sbm, [(k, inputs, outputs) for k in range(n_dmu)])在实际项目中，我发现数据质量往往比模型选择更重要。一个常见误区是过度关注模型复杂度，而忽视了基础数据的清洗和验证工作。特别是在处理非期望产出指标时，确保指标方向和数据范围的一致性至关重要。
