3、PSI计算实战：使用Python计算PSI值、pandas/numpy实现、处理缺失值与异常值

在AML模型监控中，PSI（Population Stability Index）是衡量模型评分分布稳定性的核心指标。本节将带领你从零开始，使用Python的pandas和numpy库实现PSI计算，并重点处理实际业务中常见的缺失值与异常值问题。

3.1 PSI计算原理回顾

PSI衡量的是预期分布（通常为模型开发时的基准分布）与实际分布（当前监控期的分布）之间的差异。其计算公式为：

PSI = Σ (实际占比_i - 预期占比_i) * ln(实际占比_i / 预期占比_i)

其中，i 表示第 i 个分箱（通常为10等分或20等分）。PSI值越小，分布越稳定；PSI值越大，漂移越严重。

行业经验阈值：
  • PSI < 0.1：分布稳定，无需告警
  • 0.1 ≤ PSI < 0.25：轻度漂移，需关注
  • PSI ≥ 0.25：严重漂移，触发告警

3.2 基础PSI计算函数实现

我们首先实现一个基础的PSI计算函数，使用numpy进行向量化运算，确保计算效率。

import numpy as np
import pandas as pd

def calculate_psi(expected, actual, bins=10):
    """
    计算PSI值
    
    Parameters:
    expected (array-like): 预期分布（基准分布）
    actual (array-like): 实际分布（当前分布）
    bins (int): 分箱数量，默认10
    
    Returns:
    float: PSI值
    """
    # 基于预期分布确定分箱边界
    percentiles = np.linspace(0, 100, bins + 1)
    bin_edges = np.percentile(expected, percentiles)
    
    # 处理边界重复情况（当数据中存在大量重复值时）
    bin_edges = np.unique(bin_edges)
    
    # 分箱计数
    expected_counts, _ = np.histogram(expected, bins=bin_edges)
    actual_counts, _ = np.histogram(actual, bins=bin_edges)
    
    # 计算占比（添加平滑因子避免除零）
    epsilon = 1e-10
    expected_pct = expected_counts / len(expected) + epsilon
    actual_pct = actual_counts / len(actual) + epsilon
    
    # 计算PSI
    psi = np.sum((actual_pct - expected_pct) * np.log(actual_pct / expected_pct))
    
    return psi
关键点说明：
  • 使用 np.percentile 基于预期分布确定分箱边界，确保分箱具有业务可比性
  • 使用 np.unique 处理重复边界，避免空箱
  • 添加 epsilon 平滑因子，防止对数运算中出现0值

3.3 处理缺失值与异常值

实际业务数据中，缺失值和异常值普遍存在。直接计算PSI会导致结果失真。我们需要在计算前进行预处理。

3.3.1 缺失值处理策略

策略 方法 适用场景 注意事项
删除 直接移除缺失值 缺失比例 < 5% 可能引入样本偏差
填充中位数 用预期分布的中位数填充 缺失比例 5%-20% 适用于对称分布
填充众数 用预期分布的众数填充 缺失比例 5%-20% 适用于离散型评分
模型预测 用回归模型预测缺失值 缺失比例 > 20% 计算成本高，需谨慎

3.3.2 异常值处理策略

异常值通常指评分超出合理范围（如[0,1]区间外的概率值）或极端离群点。我们采用Winsorization（截尾处理）方法。

def handle_outliers(data, lower_percentile=1, upper_percentile=99):
    """
    使用Winsorization处理异常值
    
    Parameters:
    data (array-like): 输入数据
    lower_percentile (float): 下限百分位
    upper_percentile (float): 上限百分位
    
    Returns:
    array: 处理后的数据
    """
    lower_bound = np.percentile(data, lower_percentile)
    upper_bound = np.percentile(data, upper_percentile)
    
    # 截尾处理
    data_clipped = np.clip(data, lower_bound, upper_bound)
    
    return data_clipped

3.4 完整的PSI计算流程

将上述步骤整合为一个完整的PSI计算函数，包含缺失值处理和异常值处理。

def robust_psi_calculation(expected, actual, bins=10, 
                           handle_missing='median', 
                           handle_outlier=True,
                           outlier_percentiles=(1, 99)):
    """
    鲁棒的PSI计算函数
    
    Parameters:
    expected (array-like): 预期分布
    actual (array-like): 实际分布
    bins (int): 分箱数量
    handle_missing (str): 缺失值处理策略：'drop', 'median', 'mode'
    handle_outlier (bool): 是否处理异常值
    outlier_percentiles (tuple): 异常值截尾百分位
    
    Returns:
    float: PSI值
    dict: 包含中间统计信息的字典
    """
    # 转换为numpy数组
    expected = np.array(expected, dtype=float)
    actual = np.array(actual, dtype=float)
    
    # 记录原始统计信息
    stats = {
        'expected_original_count': len(expected),
        'actual_original_count': len(actual),
        'expected_missing_count': np.isnan(expected).sum(),
        'actual_missing_count': np.isnan(actual).sum()
    }
    
    # 处理缺失值
    if handle_missing == 'drop':
        expected = expected[~np.isnan(expected)]
        actual = actual[~np.isnan(actual)]
    elif handle_missing == 'median':
        expected_median = np.nanmedian(expected)
        expected = np.where(np.isnan(expected), expected_median, expected)
        actual_median = np.nanmedian(actual)
        actual = np.where(np.isnan(actual), actual_median, actual)
    elif handle_missing == 'mode':
        # 计算众数（离散评分适用）
        expected_mode = pd.Series(expected).mode().iloc[0]
        expected = np.where(np.isnan(expected), expected_mode, expected)
        actual_mode = pd.Series(actual).mode().iloc[0]
        actual = np.where(np.isnan(actual), actual_mode, actual)
    
    # 处理异常值
    if handle_outlier:
        expected = handle_outliers(expected, 
                                   outlier_percentiles[0], 
                                   outlier_percentiles[1])
        actual = handle_outliers(actual, 
                                 outlier_percentiles[0], 
                                 outlier_percentiles[1])
    
    # 计算PSI
    psi_value = calculate_psi(expected, actual, bins)
    
    # 更新统计信息
    stats['expected_processed_count'] = len(expected)
    stats['actual_processed_count'] = len(actual)
    stats['psi_value'] = psi_value
    
    return psi_value, stats

3.5 实战案例：监控模型评分漂移

假设我们有一个信贷评分模型，开发时的基准评分分布（预期）和当前月份的评分分布（实际）如下：

# 模拟数据
np.random.seed(42)
expected_scores = np.random.beta(2, 5, 10000)  # 基准分布
actual_scores = np.random.beta(2.5, 4.5, 10000)  # 当前分布（轻度漂移）

# 引入缺失值和异常值
actual_scores[::100] = np.nan  # 1%缺失
actual_scores[::500] = 1.5     # 异常值（超出[0,1]范围）

# 计算PSI
psi_value, stats = robust_psi_calculation(
    expected_scores, 
    actual_scores, 
    bins=10,
    handle_missing='median',
    handle_outlier=True
)

print(f"PSI值: {psi_value:.4f}")
print(f"预期样本数: {stats['expected_original_count']} -> {stats['expected_processed_count']}")
print(f"实际样本数: {stats['actual_original_count']} -> {stats['actual_processed_count']}")
print(f"实际缺失值数: {stats['actual_missing_count']}")
输出示例：
PSI值: 0.1523
预期样本数: 10000 -> 10000
实际样本数: 10000 -> 10000
实际缺失值数: 100
注意事项：
  • 缺失值处理策略的选择应基于业务理解，避免机械套用
  • 异常值截尾的百分位选择需结合评分分布特征，通常1%-99%是安全范围
  • 分箱数量（bins）建议保持与模型开发时一致，通常为10或20
  • 当预期分布和实际分布样本量差异较大时，需关注统计显著性

3.6 PSI计算流程可视化

下图展示了完整的PSI计算流程，从原始数据输入到最终PSI值输出，包含缺失值处理和异常值处理两个关键预处理步骤。

PSI计算完整流程图 原始数据输入 步骤1：缺失值处理 删除 / 中位数填充 / 众数填充 步骤2：异常值处理 Winsorization截尾（1%-99%） 步骤3：等频分箱（10箱） 步骤4：PSI公式计算 输出PSI值 预处理关键点 缺失比例<5%：删除 缺失比例5-20%：填充 异常值：截尾处理 分箱注意事项 基于预期分布分箱 避免空箱（unique处理） 添加平滑因子ε

3.7 批量监控与告警集成

在实际AML系统中，我们需要对多个模型、多个时间窗口进行批量PSI计算，并集成告警逻辑。

def batch_psi_monitoring(reference_data, current_data_dict, 
                         model_name, bins=10, threshold_warning=0.1,
                         threshold_alarm=0.25):
    """
    批量PSI监控函数
    
    Parameters:
    reference_data (array): 基准数据
    current_data_dict (dict): 键为时间窗口，值为当前数据
    model_name (str): 模型名称
    bins (int): 分箱数量
    threshold_warning (float): 警告阈值
    threshold_alarm (float): 告警阈值
    
    Returns:
    pd.DataFrame: 监控结果表
    """
    results = []
    
    for window_name, current_data in current_data_dict.items():
        psi_value, stats = robust_psi_calculation(
            reference_data, current_data, bins=bins
        )
        
        # 判断告警级别
        if psi_value >= threshold_alarm:
            alert_level = '严重告警'
        elif psi_value >= threshold_warning:
            alert_level = '轻度警告'
        else:
            alert_level = '正常'
        
        results.append({
            '模型名称': model_name,
            '时间窗口': window_name,
            'PSI值': round(psi_value, 4),
            '告警级别': alert_level,
            '基准样本数': stats['expected_processed_count'],
            '当前样本数': stats['actual_processed_count']
        })
    
    return pd.DataFrame(results)

# 使用示例
reference = expected_scores
current_windows = {
    '2024-01': actual_scores,
    '2024-02': np.random.beta(3, 4, 10000),  # 更明显的漂移
    '2024-03': np.random.beta(2, 5, 10000)   # 稳定
}

monitor_df = batch_psi_monitoring(reference, current_windows, '信贷评分模型_v2')
print(monitor_df.to_string(index=False))
输出示例：
模型名称        时间窗口   PSI值  告警级别  基准样本数  当前样本数
信贷评分模型_v2  2024-01  0.1523  轻度警告    10000    10000
信贷评分模型_v2  2024-02  0.3127  严重告警    10000    10000
信贷评分模型_v2  2024-03  0.0421    正常      10000    10000
最佳实践总结：
  • 数据预处理优先：缺失值和异常值处理必须在PSI计算之前完成，否则结果不可靠
  • 分箱策略一致：基准分布和当前分布使用相同的分箱边界，确保可比性
  • 阈值动态调整：根据模型特性和业务容忍度，可适当调整0.1和0.25的阈值
  • 批量监控自动化：将PSI计算集成到定时任务中，实现自动化漂移检测
← 上一章 📖 返回目录 下一章 →