用Python的SciPy库5分钟搞定超效率SBM模型（含非期望产出处理）

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')

关键预处理步骤：

1. 数据标准化（特别是当不同指标量纲差异大时）
2. 非期望产出指标取负值（因其与效率负相关）
3. 检查数据完整性，处理缺失值

注意：非期望产出指标在建模时需要特殊处理，通常采用方向距离函数或将其视为"负产出"。

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. 效率值>1表示该DMU处于生产前沿面上方（超高效）
2. 效率值<1表示存在改进空间
3. CRS与VRS结果的差异反映规模效率影响

4. 高级应用与常见问题处理

在实际应用中，我们经常会遇到各种特殊情况。以下是几个常见问题的解决方案：

4.1 处理零值或负值数据

当数据中存在零或负值时，传统DEA模型可能失效。解决方法包括：

• 对数据进行平移变换
• 使用特殊的距离函数
• 采用非径向的SBM模型

# 数据平移示例（当存在负值时） shift = np.abs(np.min(inputs)) + 1 shifted_inputs = inputs + shift

4.2 模型选择指南

不同场景下的模型选择策略：

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)])

在实际项目中，我发现数据质量往往比模型选择更重要。一个常见误区是过度关注模型复杂度，而忽视了基础数据的清洗和验证工作。特别是在处理非期望产出指标时，确保指标方向和数据范围的一致性至关重要。