用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)关键预处理步骤数据标准化特别是当不同指标量纲差异大时非期望产出指标取负值因其与效率负相关检查数据完整性处理缺失值注意非期望产出指标在建模时需要特殊处理通常采用方向距离函数或将其视为负产出。2. 超效率SBM模型核心实现超效率SBM模型的核心是求解一系列线性规划问题。我们利用SciPy的linprog函数来实现这一过程相比传统DEA模型它能处理松弛变量并允许效率值超过1。以下是完整的模型实现代码import numpy as np from scipy.optimize import linprog def super_sbm(inputs, good_outputs, bad_outputs, crsTrue): 考虑非期望产出的超效率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_ubA_ub, b_ubb_ub, boundsbounds, methodhighs) if res.success: efficiencies[k] 1 / res.x[-1] else: efficiencies[k] np.nan return efficiencies关键参数说明crs决定是否采用规模报酬不变假设methodhighs使用高性能线性规划求解器bounds确保权重非负但允许效率值超过13. 实战案例区域碳排放效率评估让我们用一个实际案例演示如何使用上述函数。假设我们有四个地区的投入产出数据# 投入指标劳动力(万人), 资本(亿元), 能源(万吨标准煤) 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, crsTrue) eff_vrs super_sbm(inputs, good_outputs, bad_outputs, crsFalse) 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, alpha0.7, labelCRS) plt.bar(range(len(eff_vrs)), eff_vrs, alpha0.7, labelVRS) plt.axhline(1, colorred, linestyle--) plt.xlabel(DMU) plt.ylabel(Efficiency Score) plt.title(Efficiency Comparison (CRS vs VRS)) plt.legend() plt.show()4.4 性能优化建议当处理大量DMU时可采取以下优化措施使用稀疏矩阵存储约束条件并行化计算各DMU的效率采用更高效的求解器如Gurobi或CPLEXfrom 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)])在实际项目中我发现数据质量往往比模型选择更重要。一个常见误区是过度关注模型复杂度而忽视了基础数据的清洗和验证工作。特别是在处理非期望产出指标时确保指标方向和数据范围的一致性至关重要。