数据包络分析法(DEA)

“效率不在于你做了多少事，而在于你用最少的资源做出了最好的成果。” —— 管理学家彼得·德鲁克

数据包络分析法（Data Envelopment Analysis, DEA）是一种基于线性规划的非参数效率评价方法，由Charnes、Cooper和Rhodes于1978年首次提出。DEA无需预先设定投入产出之间的函数关系，通过数学规划模型直接计算决策单元的相对效率，是评价具有多投入多产出系统效率的强有力工具。

基本原理

决策单元(DMU)的概念

在DEA方法中，被评价的对象称为决策单元（Decision Making Unit, DMU）。每个DMU都消耗一定数量的投入（inputs）来产生一定数量的产出（outputs）。

例如：

• 医院：投入为医生数量、床位数、设备经费；产出为门诊量、治愈人数
• 学校：投入为教师数量、经费投入；产出为毕业生人数、就业率
• 银行网点：投入为员工数、运营成本；产出为存款额、利润

设有 \(n\) 个DMU，每个DMU有 \(m\) 种投入和 \(s\) 种产出。对于第 \(j\) 个DMU（\(j=1,2,\ldots,n\)）：

• 投入向量：\(\mathbf{x}j = (x{1j}, x_{2j}, \ldots, x_{mj})^T\)
• 产出向量：\(\mathbf{y}j = (y{1j}, y_{2j}, \ldots, y_{sj})^T\)

其中 \(x_{ij} > 0\) 表示第 \(j\) 个DMU对第 \(i\) 种投入的消耗量，\(y_{rj} > 0\) 表示第 \(j\) 个DMU对第 \(r\) 种产出的产生量。

效率前沿面

DEA的核心思想是构造效率前沿面（Efficient Frontier），以此为基准评价各DMU的相对效率。前沿面由表现最优的DMU构成，位于前沿面上的DMU效率值为1（DEA有效），位于前沿面内部的效率值小于1（DEA无效）。

直观理解：假设只有一种投入和一种产出，效率前沿就是“单位投入下产出最大“的包络线。多投入多产出时，前沿面是高维空间中的包络超曲面。

相对效率的度量

DEA度量的是相对效率而非绝对效率。即某个DMU的效率是相对于样本中最优DMU而言的。效率值 \(\theta\) 的含义为：

• \(\theta = 1\)：该DMU位于效率前沿面上，称为DEA有效
• \(\theta < 1\)：该DMU未达到效率前沿面，存在效率损失，可通过等比例缩减投入 \((1-\theta)\) 的比例达到前沿面

CCR模型

模型介绍

CCR模型由Charnes、Cooper和Rhodes于1978年提出，假设规模报酬不变（Constant Returns to Scale, CRS）。该模型评价的是技术效率（Technical Efficiency）与规模效率（Scale Efficiency）的综合效率。

分数规划形式

对于被评价的第 \(k\) 个DMU（记为 \(\text{DMU}_k\)），CCR模型的原始形式为分数规划：

\[ \max \quad \frac{\sum_{r=1}^{s} u_r y_{rk}}{\sum_{i=1}^{m} v_i x_{ik}} \]

\[ \text{s.t.} \quad \frac{\sum_{r=1}^{s} u_r y_{rj}}{\sum_{i=1}^{m} v_i x_{ij}} \leq 1, \quad j=1,2,\ldots,n \]

\[ u_r \geq 0, \quad r=1,2,\ldots,s \]

\[ v_i \geq 0, \quad i=1,2,\ldots,m \]

其中 \(u_r\) 为第 \(r\) 种产出的权重，\(v_i\) 为第 \(i\) 种投入的权重。模型的含义是：在所有DMU的效率比值不超过1的约束下，寻找最有利于 \(\text{DMU}_k\) 的权重组合，使其效率最大化。

线性规划形式（对偶形式）

通过Charnes-Cooper变换，将分数规划转化为线性规划。实际计算中通常使用对偶形式（投入导向）：

\[ \min \quad \theta \]

\[ \text{s.t.} \quad \sum_{j=1}^{n} \lambda_j x_{ij} \leq \theta x_{ik}, \quad i=1,2,\ldots,m \]

\[ \sum_{j=1}^{n} \lambda_j y_{rj} \geq y_{rk}, \quad r=1,2,\ldots,s \]

\[ \lambda_j \geq 0, \quad j=1,2,\ldots,n \]

其中：

• \(\theta\) 为 \(\text{DMU}_k\) 的效率值，\(0 < \theta \leq 1\)
• \(\lambda_j\) 为构造虚拟DMU的组合系数

效率判定准则

设CCR模型的最优解为 \(\theta^, \lambda^\)：

1. 若 \(\theta^* = 1\) 且所有松弛变量为零，则 \(\text{DMU}_k\) DEA有效（技术有效且规模有效）
2. 若 \(\theta^* = 1\) 但存在非零松弛变量，则 \(\text{DMU}_k\) 弱DEA有效
3. 若 \(\theta^* < 1\)，则 \(\text{DMU}_k\) DEA无效

BCC模型

模型介绍

BCC模型由Banker、Charnes和Cooper于1984年提出，假设规模报酬可变（Variable Returns to Scale, VRS）。BCC模型将CCR模型的综合效率分解为纯技术效率和规模效率：

\[ \text{综合效率(CCR)} = \text{纯技术效率(BCC)} \times \text{规模效率} \]

线性规划形式

BCC模型在CCR对偶模型的基础上增加了凸性约束 \(\sum_{j=1}^{n}\lambda_j = 1\)：

\[ \min \quad \theta \]

\[ \text{s.t.} \quad \sum_{j=1}^{n} \lambda_j x_{ij} \leq \theta x_{ik}, \quad i=1,2,\ldots,m \]

\[ \sum_{j=1}^{n} \lambda_j y_{rj} \geq y_{rk}, \quad r=1,2,\ldots,s \]

\[ \sum_{j=1}^{n} \lambda_j = 1 \]

\[ \lambda_j \geq 0, \quad j=1,2,\ldots,n \]

规模报酬判断

通过BCC模型还可以判断各DMU的规模报酬状况：

• 若最优解中 \(\sum_{j=1}^{n}\lambda_j^* < 1\)，则 \(\text{DMU}_k\) 处于规模报酬递增阶段
• 若最优解中 \(\sum_{j=1}^{n}\lambda_j^* = 1\)，则 \(\text{DMU}_k\) 处于规模报酬不变阶段
• 若最优解中 \(\sum_{j=1}^{n}\lambda_j^* > 1\)，则 \(\text{DMU}_k\) 处于规模报酬递减阶段

CCR与BCC模型的关系

计算步骤

DEA分析的完整计算流程如下：

第一步：确定DMU集合与投入产出指标

选择具有可比性的DMU集合（建议DMU数量不少于投入产出指标总数的2倍），确定投入指标（越小越好的资源消耗类）和产出指标（越大越好的成果产出类）。

第二步：数据收集与整理

收集各DMU的投入产出数据，检查非负性要求（DEA要求所有数据严格为正），必要时进行数据预处理。

第三步：构建线性规划模型

对每一个 \(\text{DMU}_k\)（\(k=1,2,\ldots,n\)），分别构建CCR或BCC线性规划模型。

第四步：求解线性规划

使用单纯形法或内点法求解，得到最优效率值 \(\theta^\) 和最优组合系数 \(\lambda^\)。

第五步：效率分析与排序

根据 \(\theta^\) 值判断效率状态；计算规模效率 \(\text{SE} = \theta_{CCR}^ / \theta_{BCC}^*\)；分析松弛变量确定投入冗余和产出不足；判断规模报酬阶段。

第六步：结果解释与改进建议

对DEA无效的DMU提出改进方向，确定参考集（benchmarking），计算投影值（目标值）。

实际案例分析

问题描述

某地区有8家社区医院，管理部门希望评估各医院的运营效率。已知各医院的投入产出数据如下：

投入指标：

• \(x_1\)：医护人员数（人）
• \(x_2\)：床位数（张）
• \(x_3\)：年运营经费（万元）

产出指标：

• \(y_1\)：年门诊量（万人次）
• \(y_2\)：年住院人次（人次）
• \(y_3\)：患者满意度评分（分）

各医院数据如下表所示：

建模分析

以医院A为例（\(k=1\)），构建CCR对偶模型：

\[ \min \quad \theta \]

\[ \text{s.t.} \quad 105\lambda_1 + 90\lambda_2 + 120\lambda_3 + 75\lambda_4 + 130\lambda_5 + 85\lambda_6 + 95\lambda_7 + 110\lambda_8 \leq 105\theta \]

\[ 80\lambda_1 + 60\lambda_2 + 100\lambda_3 + 50\lambda_4 + 110\lambda_5 + 55\lambda_6 + 70\lambda_7 + 90\lambda_8 \leq 80\theta \]

\[ 1200\lambda_1 + 900\lambda_2 + 1500\lambda_3 + 750\lambda_4 + 1800\lambda_5 + 800\lambda_6 + 1000\lambda_7 + 1400\lambda_8 \leq 1200\theta \]

\[ 12.5\lambda_1 + 10.8\lambda_2 + 13.0\lambda_3 + 9.5\lambda_4 + 11.0\lambda_5 + 10.2\lambda_6 + 11.5\lambda_7 + 12.0\lambda_8 \geq 12.5 \]

\[ 3200\lambda_1 + 2800\lambda_2 + 3500\lambda_3 + 2400\lambda_4 + 3000\lambda_5 + 2600\lambda_6 + 3100\lambda_7 + 3300\lambda_8 \geq 3200 \]

\[ 92\lambda_1 + 88\lambda_2 + 85\lambda_3 + 90\lambda_4 + 82\lambda_5 + 91\lambda_6 + 87\lambda_7 + 86\lambda_8 \geq 92 \]

\[ \lambda_j \geq 0, \quad j=1,2,\ldots,8 \]

对所有8家医院分别求解上述模型，即可得到各医院的CCR效率值。类似地，增加约束 \(\sum_{j=1}^{8}\lambda_j = 1\) 即为BCC模型。

Python代码实现

以下代码使用 scipy.optimize.linprog 求解DEA模型，完整实现CCR和BCC两种模型的计算。

import numpy as np
from scipy.optimize import linprog

def dea_ccr(inputs, outputs):
    """
    CCR模型（投入导向）：计算各DMU的综合效率
    
    参数:
        inputs: numpy数组, 形状为 (n, m), n个DMU的m种投入
        outputs: numpy数组, 形状为 (n, s), n个DMU的s种产出
    
    返回:
        efficiencies: 各DMU的效率值列表
        lambdas: 各DMU的最优组合系数
    """
    n, m = inputs.shape  # n个DMU, m种投入
    _, s = outputs.shape  # s种产出
    
    efficiencies = []
    lambdas_list = []
    
    for k in range(n):
        # 决策变量: [theta, lambda_1, lambda_2, ..., lambda_n]
        # 目标函数: min theta
        c = np.zeros(n + 1)
        c[0] = 1  # theta的系数为1
        
        # 不等式约束: A_ub @ x <= b_ub
        # 投入约束: sum(lambda_j * x_ij) - theta * x_ik <= 0
        # 即: -x_ik * theta + sum(lambda_j * x_ij) <= 0
        A_ub = np.zeros((m, n + 1))
        b_ub = np.zeros(m)
        
        for i in range(m):
            A_ub[i, 0] = -inputs[k, i]  # -theta * x_ik
            for j in range(n):
                A_ub[i, j + 1] = inputs[j, i]  # lambda_j * x_ij
        
        # 产出约束: sum(lambda_j * y_rj) >= y_rk
        # 转化为: -sum(lambda_j * y_rj) <= -y_rk
        A_ub2 = np.zeros((s, n + 1))
        b_ub2 = np.zeros(s)
        
        for r in range(s):
            A_ub2[r, 0] = 0  # theta不参与产出约束
            for j in range(n):
                A_ub2[r, j + 1] = -outputs[j, r]
            b_ub2[r] = -outputs[k, r]
        
        # 合并不等式约束
        A_ub_full = np.vstack([A_ub, A_ub2])
        b_ub_full = np.concatenate([b_ub, b_ub2])
        
        # 变量边界: theta无上界, lambda_j >= 0
        bounds = [(None, None)] + [(0, None)] * n  # theta可以为任意值(但最优解在0-1)
        
        # 求解线性规划
        result = linprog(c, A_ub=A_ub_full, b_ub=b_ub_full, bounds=bounds, method='highs')
        
        if result.success:
            efficiencies.append(result.x[0])
            lambdas_list.append(result.x[1:])
        else:
            efficiencies.append(None)
            lambdas_list.append(None)
    
    return efficiencies, lambdas_list


def dea_bcc(inputs, outputs):
    """
    BCC模型（投入导向）：计算各DMU的纯技术效率
    
    参数:
        inputs: numpy数组, 形状为 (n, m), n个DMU的m种投入
        outputs: numpy数组, 形状为 (n, s), n个DMU的s种产出
    
    返回:
        efficiencies: 各DMU的效率值列表
        lambdas: 各DMU的最优组合系数
    """
    n, m = inputs.shape
    _, s = outputs.shape
    
    efficiencies = []
    lambdas_list = []
    
    for k in range(n):
        # 决策变量: [theta, lambda_1, lambda_2, ..., lambda_n]
        c = np.zeros(n + 1)
        c[0] = 1
        
        # 投入约束
        A_ub = np.zeros((m, n + 1))
        b_ub = np.zeros(m)
        
        for i in range(m):
            A_ub[i, 0] = -inputs[k, i]
            for j in range(n):
                A_ub[i, j + 1] = inputs[j, i]
        
        # 产出约束
        A_ub2 = np.zeros((s, n + 1))
        b_ub2 = np.zeros(s)
        
        for r in range(s):
            for j in range(n):
                A_ub2[r, j + 1] = -outputs[j, r]
            b_ub2[r] = -outputs[k, r]
        
        A_ub_full = np.vstack([A_ub, A_ub2])
        b_ub_full = np.concatenate([b_ub, b_ub2])
        
        # BCC模型的等式约束: sum(lambda_j) = 1
        A_eq = np.zeros((1, n + 1))
        A_eq[0, 1:] = 1  # lambda的系数为1, theta的系数为0
        b_eq = np.array([1.0])
        
        bounds = [(None, None)] + [(0, None)] * n
        
        result = linprog(c, A_ub=A_ub_full, b_ub=b_ub_full,
                        A_eq=A_eq, b_eq=b_eq, bounds=bounds, method='highs')
        
        if result.success:
            efficiencies.append(result.x[0])
            lambdas_list.append(result.x[1:])
        else:
            efficiencies.append(None)
            lambdas_list.append(None)
    
    return efficiencies, lambdas_list


def analyze_dea_results(inputs, outputs, dmu_names=None):
    """
    综合DEA分析：同时计算CCR和BCC效率，并给出规模效率和规模报酬判断
    
    参数:
        inputs: 投入数据矩阵
        outputs: 产出数据矩阵
        dmu_names: DMU名称列表
    """
    n = inputs.shape[0]
    if dmu_names is None:
        dmu_names = [f"DMU{i+1}" for i in range(n)]
    
    # 计算CCR效率（综合效率）
    ccr_eff, ccr_lambdas = dea_ccr(inputs, outputs)
    
    # 计算BCC效率（纯技术效率）
    bcc_eff, bcc_lambdas = dea_bcc(inputs, outputs)
    
    # 计算规模效率
    scale_eff = []
    for i in range(n):
        if ccr_eff[i] is not None and bcc_eff[i] is not None and bcc_eff[i] > 0:
            scale_eff.append(ccr_eff[i] / bcc_eff[i])
        else:
            scale_eff.append(None)
    
    # 判断规模报酬
    returns_to_scale = []
    for i in range(n):
        if bcc_lambdas[i] is not None:
            lambda_sum = np.sum(bcc_lambdas[i])
            if abs(lambda_sum - 1.0) < 1e-6:
                returns_to_scale.append("不变")
            elif lambda_sum < 1.0:
                returns_to_scale.append("递增")
            else:
                returns_to_scale.append("递减")
        else:
            returns_to_scale.append("未知")
    
    # 打印结果
    print("=" * 80)
    print(f"{'DEA效率分析结果':^80}")
    print("=" * 80)
    print(f"\n{'DMU':<8}{'综合效率(CCR)':<15}{'纯技术效率(BCC)':<17}"
          f"{'规模效率':<12}{'规模报酬':<10}{'是否有效':<10}")
    print("-" * 80)
    
    for i in range(n):
        ccr = f"{ccr_eff[i]:.4f}" if ccr_eff[i] is not None else "N/A"
        bcc = f"{bcc_eff[i]:.4f}" if bcc_eff[i] is not None else "N/A"
        se = f"{scale_eff[i]:.4f}" if scale_eff[i] is not None else "N/A"
        rts = returns_to_scale[i]
        
        # 判断是否DEA有效
        if ccr_eff[i] is not None and abs(ccr_eff[i] - 1.0) < 1e-6:
            status = "有效"
        else:
            status = "无效"
        
        print(f"{dmu_names[i]:<8}{ccr:<15}{bcc:<17}{se:<12}{rts:<10}{status:<10}")
    
    print("-" * 80)
    
    # 对DEA无效的DMU给出改进建议
    print("\n改进建议：")
    print("-" * 80)
    for i in range(n):
        if ccr_eff[i] is not None and ccr_eff[i] < 1.0 - 1e-6:
            print(f"\n{dmu_names[i]}（效率值={ccr_eff[i]:.4f}）：")
            print(f"  投入应缩减比例: {(1-ccr_eff[i])*100:.2f}%")
            print(f"  具体目标值:")
            for j in range(inputs.shape[1]):
                target = ccr_eff[i] * inputs[i, j]
                print(f"    投入{j+1}: 当前={inputs[i,j]:.1f}, "
                      f"目标={target:.1f}, 缩减={inputs[i,j]-target:.1f}")
    
    return ccr_eff, bcc_eff, scale_eff, returns_to_scale


# ============================================================
# 实际案例：8家社区医院效率评价
# ============================================================

if __name__ == "__main__":
    # 投入数据: 医护人员数, 床位数, 运营经费(万元)
    inputs = np.array([
        [105, 80, 1200],   # 医院A
        [90,  60,  900],   # 医院B
        [120, 100, 1500],  # 医院C
        [75,  50,  750],   # 医院D
        [130, 110, 1800],  # 医院E
        [85,  55,  800],   # 医院F
        [95,  70, 1000],   # 医院G
        [110, 90, 1400],   # 医院H
    ], dtype=float)
    
    # 产出数据: 门诊量(万人次), 住院人次, 满意度(分)
    outputs = np.array([
        [12.5, 3200, 92],  # 医院A
        [10.8, 2800, 88],  # 医院B
        [13.0, 3500, 85],  # 医院C
        [9.5,  2400, 90],  # 医院D
        [11.0, 3000, 82],  # 医院E
        [10.2, 2600, 91],  # 医院F
        [11.5, 3100, 87],  # 医院G
        [12.0, 3300, 86],  # 医院H
    ], dtype=float)
    
    dmu_names = ["医院A", "医院B", "医院C", "医院D",
                 "医院E", "医院F", "医院G", "医院H"]
    
    # 执行DEA分析
    ccr_eff, bcc_eff, scale_eff, rts = analyze_dea_results(
        inputs, outputs, dmu_names
    )
    
    # 绘制效率对比图
    try:
        import matplotlib.pyplot as plt
        
        plt.rcParams['font.sans-serif'] = ['SimHei']
        plt.rcParams['axes.unicode_minus'] = False
        
        x = range(len(dmu_names))
        width = 0.25
        
        fig, ax = plt.subplots(figsize=(10, 6))
        bars1 = ax.bar([i - width for i in x], ccr_eff, width, label='综合效率(CCR)')
        bars2 = ax.bar(x, bcc_eff, width, label='纯技术效率(BCC)')
        bars3 = ax.bar([i + width for i in x], scale_eff, width, label='规模效率')
        
        ax.set_xlabel('医院')
        ax.set_ylabel('效率值')
        ax.set_title('各医院DEA效率分析对比')
        ax.set_xticks(x)
        ax.set_xticklabels(dmu_names)
        ax.legend()
        ax.set_ylim(0, 1.1)
        ax.axhline(y=1.0, color='r', linestyle='--', alpha=0.5)
        
        plt.tight_layout()
        plt.savefig('dea_efficiency_comparison.png', dpi=150, bbox_inches='tight')
        plt.show()
        print("\n效率对比图已保存为 dea_efficiency_comparison.png")
    except ImportError:
        print("\n提示: 安装matplotlib后可生成效率对比图")

代码运行说明

确保安装 numpy 和 scipy，运行后将输出各医院的综合效率、纯技术效率和规模效率，并自动给出改进建议。如安装了 matplotlib，还会生成效率对比柱状图。

结果解读

• 综合效率为1的医院（如医院D、医院F）：投入产出配比最优，是其他医院的标杆
• 纯技术效率为1但规模效率小于1：管理水平最优，但规模不合适，需调整规模
• 纯技术效率小于1（如医院E）：管理水平有待提高，存在资源浪费
• 规模报酬递增：应适当扩大规模以提高效率
• 规模报酬递减：规模过大，应适当缩小或分流

应用注意事项与局限性

应用注意事项

1. DMU的同质性要求

所有被评价的DMU必须具有可比性：使用相同类型的投入和产出，在相似的环境下运营。例如，不应将三甲医院与社区诊所放在一起评价。

2. 指标选择原则

• 投入产出指标应具有明确的经济含义
• 指标之间不应存在严格的线性关系
• 投入指标应是“越少越好“，产出指标应是“越多越好“
• 建议DMU数量 \(\geq 2(m+s)\)，其中 \(m\) 为投入数，\(s\) 为产出数

3. 数据要求

• 所有投入产出数据必须为正数
• 如存在零值或负值，需进行适当的数据变换（如平移）
• 数据应具有一定的区分度，避免过于集中

4. 模型选择建议

若关注综合效率（含规模），使用CCR模型；若只关注管理效率、排除规模影响，使用BCC模型。实际中建议同时计算两种模型，综合分析。

5. 权重灵活性的双面性

DEA允许每个DMU选择对自己最有利的权重，这既是优点（客观性），也可能导致某些DMU通过极端权重获得高效率值，需结合实际判断。

局限性

1. 仅能评价相对效率

DEA评价的是相对效率，所有DMU都DEA有效并不意味着它们都达到了绝对最优水平，只说明在当前样本中没有更好的参照。

2. 对异常值敏感

个别异常DMU可能显著影响效率前沿面，建议在分析前进行异常值检测。

3. 无法进行统计检验

DEA是非参数方法，难以对效率差异进行显著性检验。可结合Bootstrap方法弥补。

4. 维度灾难

当投入产出指标数量相对DMU数量过多时，大部分DMU都会被评为有效，失去区分度。经验法则建议：

\[ n \geq \max{m \times s, ; 3(m+s)} \]

5. 不能处理随机误差

传统DEA将偏离前沿面的部分都归为效率损失，无法区分无效率与随机噪声。可结合随机前沿分析（SFA）弥补。

6. 投入产出方向性限制

标准DEA要求投入越少越好、产出越多越好。对于非期望产出（如污染），需使用方向性距离函数等特殊处理方法。

扩展方法

针对上述局限性，学者们提出了多种改进方法：超效率DEA（对有效DMU进一步排序）、Malmquist指数（动态效率变化分析）、网络DEA（考虑DMU内部子过程）、Bootstrap-DEA（统计推断）、非径向DEA/SBM模型（处理松弛问题）等。