3. PSI计算实战（Python实现）

本节将使用 pandas 与 numpy 实现群体稳定性指标（PSI）的完整计算流程。你将掌握从原始分布到最终 PSI 值的全链路编码，并学会处理实际业务中常见的缺失值与异常值问题。

3.1 PSI计算核心逻辑

PSI 的计算本质是衡量两个分布（预期分布 P 与实际分布 Q）在每个分箱上的偏差。其数学公式为：

PSI = Σ ( (Q_i - P_i) * ln(Q_i / P_i) )

其中 P_i 为预期分布在第 i 个分箱中的占比，Q_i 为实际分布在第 i 个分箱中的占比。计算流程包含三个关键步骤：

1. 分箱：将连续变量或离散变量划分为 k 个区间（通常为 10 等分）。
2. 统计占比：分别计算预期样本与实际样本在每个分箱中的数量占比。
3. 计算偏差：对每个分箱执行 (Q_i - P_i) * ln(Q_i / P_i) 并求和。

3.2 使用 pandas 与 numpy 实现 PSI

以下代码展示了从原始数据到 PSI 值的完整实现。我们使用 等宽分箱 策略，将变量分为 10 个区间。

import pandas as pd
import numpy as np

def calculate_psi(expected, actual, bins=10):
    """
    计算群体稳定性指标 (PSI)
    
    参数:
        expected: array-like, 预期分布（如训练集评分）
        actual: array-like, 实际分布（如当前生产评分）
        bins: int, 分箱数量（默认10）
    
    返回:
        psi: float, PSI值
        psi_df: DataFrame, 每个分箱的详细计算过程
    """
    # 1. 确定分箱边界（基于预期分布）
    min_val = min(expected.min(), actual.min())
    max_val = max(expected.max(), actual.max())
    # 添加微小偏移避免边界值落在区间外
    bin_edges = np.linspace(min_val - 1e-6, max_val + 1e-6, bins + 1)
    
    # 2. 统计预期分布占比
    expected_counts, _ = np.histogram(expected, bins=bin_edges)
    expected_pct = expected_counts / len(expected)
    
    # 3. 统计实际分布占比
    actual_counts, _ = np.histogram(actual, bins=bin_edges)
    actual_pct = actual_counts / len(actual)
    
    # 4. 处理零占比（避免除零和ln(0)）
    expected_pct = np.where(expected_pct == 0, 1e-6, expected_pct)
    actual_pct = np.where(actual_pct == 0, 1e-6, actual_pct)
    
    # 5. 计算每个分箱的PSI贡献
    psi_components = (actual_pct - expected_pct) * np.log(actual_pct / expected_pct)
    psi = np.sum(psi_components)
    
    # 6. 构建结果DataFrame
    psi_df = pd.DataFrame({
        'bin_start': bin_edges[:-1],
        'bin_end': bin_edges[1:],
        'expected_count': expected_counts,
        'expected_pct': expected_pct,
        'actual_count': actual_counts,
        'actual_pct': actual_pct,
        'psi_component': psi_components
    })
    
    return psi, psi_df

# 示例用法
np.random.seed(42)
expected_scores = np.random.normal(600, 100, 10000)
actual_scores = np.random.normal(620, 120, 10000)

psi_value, psi_detail = calculate_psi(expected_scores, actual_scores)
print(f"PSI值: {psi_value:.4f}")
print(psi_detail.head())

3.3 自定义 PSI 函数：支持多种分箱策略

实际业务中，等宽分箱可能不适用于长尾分布。我们扩展函数，支持 等频分箱 和 自定义分箱。

def calculate_psi_flexible(expected, actual, bins=10, strategy='equal_width', custom_edges=None):
    """
    支持多种分箱策略的PSI计算函数
    
    参数:
        strategy: str, 分箱策略 ('equal_width', 'equal_freq', 'custom')
        custom_edges: list, 自定义分箱边界（仅strategy='custom'时使用）
    """
    if strategy == 'equal_width':
        # 等宽分箱
        min_val = min(expected.min(), actual.min())
        max_val = max(expected.max(), actual.max())
        bin_edges = np.linspace(min_val - 1e-6, max_val + 1e-6, bins + 1)
    
    elif strategy == 'equal_freq':
        # 等频分箱（基于预期分布的分位数）
        percentiles = np.linspace(0, 100, bins + 1)
        bin_edges = np.percentile(expected, percentiles)
        bin_edges[0] -= 1e-6
        bin_edges[-1] += 1e-6
    
    elif strategy == 'custom':
        # 自定义分箱
        if custom_edges is None:
            raise ValueError("使用自定义分箱时必须提供custom_edges参数")
        bin_edges = np.array(custom_edges)
    
    else:
        raise ValueError(f"不支持的分箱策略: {strategy}")
    
    # 统计占比
    expected_counts, _ = np.histogram(expected, bins=bin_edges)
    actual_counts, _ = np.histogram(actual, bins=bin_edges)
    
    expected_pct = expected_counts / len(expected)
    actual_pct = actual_counts / len(actual)
    
    # 平滑处理
    expected_pct = np.where(expected_pct == 0, 1e-6, expected_pct)
    actual_pct = np.where(actual_pct == 0, 1e-6, actual_pct)
    
    # 计算PSI
    psi = np.sum((actual_pct - expected_pct) * np.log(actual_pct / expected_pct))
    
    return psi

# 三种策略对比
psi_equal_width = calculate_psi_flexible(expected_scores, actual_scores, strategy='equal_width')
psi_equal_freq = calculate_psi_flexible(expected_scores, actual_scores, strategy='equal_freq')
psi_custom = calculate_psi_flexible(expected_scores, actual_scores, strategy='custom', 
                                    custom_edges=[0, 300, 400, 500, 600, 700, 800, 1000])

print(f"等宽分箱 PSI: {psi_equal_width:.4f}")
print(f"等频分箱 PSI: {psi_equal_freq:.4f}")
print(f"自定义分箱 PSI: {psi_custom:.4f}")

3.4 处理缺失值与异常值

实际数据中缺失值和异常值会严重影响 PSI 计算的准确性。以下是三种常见处理策略：

def calculate_psi_robust(expected, actual, bins=10, missing_strategy='separate_bin', 
                         clip_outliers=True):
    """
    处理缺失值和异常值的鲁棒PSI计算
    
    参数:
        missing_strategy: str, 缺失值处理策略 ('drop', 'separate_bin', 'impute_mean')
        clip_outliers: bool, 是否截断异常值
    """
    # 复制数据避免修改原始数据
    exp = expected.copy()
    act = actual.copy()
    
    # 处理缺失值
    if missing_strategy == 'drop':
        exp = exp.dropna()
        act = act.dropna()
    elif missing_strategy == 'separate_bin':
        # 记录缺失标记，后续单独成箱
        exp_missing = exp.isna().sum()
        act_missing = act.isna().sum()
        exp = exp.dropna()
        act = act.dropna()
    elif missing_strategy == 'impute_mean':
        mean_val = pd.concat([exp, act]).mean()
        exp = exp.fillna(mean_val)
        act = act.fillna(mean_val)
    
    # 处理异常值（截断到1%和99%分位数）
    if clip_outliers:
        combined = pd.concat([exp, act])
        lower_bound = combined.quantile(0.01)
        upper_bound = combined.quantile(0.99)
        exp = exp.clip(lower_bound, upper_bound)
        act = act.clip(lower_bound, upper_bound)
    
    # 计算分箱边界
    min_val = min(exp.min(), act.min())
    max_val = max(exp.max(), act.max())
    bin_edges = np.linspace(min_val - 1e-6, max_val + 1e-6, bins + 1)
    
    # 统计占比
    exp_counts, _ = np.histogram(exp, bins=bin_edges)
    act_counts, _ = np.histogram(act, bins=bin_edges)
    
    # 如果使用separate_bin策略，将缺失值作为最后一个分箱
    if missing_strategy == 'separate_bin':
        exp_counts = np.append(exp_counts, exp_missing)
        act_counts = np.append(act_counts, act_missing)
    
    exp_pct = exp_counts / len(expected)
    act_pct = act_counts / len(actual)
    
    # 平滑处理
    exp_pct = np.where(exp_pct == 0, 1e-6, exp_pct)
    act_pct = np.where(act_pct == 0, 1e-6, act_pct)
    
    psi = np.sum((act_pct - exp_pct) * np.log(act_pct / exp_pct))
    
    return psi

# 测试含缺失值的数据
np.random.seed(42)
exp_with_na = pd.Series(np.random.normal(600, 100, 10000))
act_with_na = pd.Series(np.random.normal(620, 120, 10000))
# 随机引入5%缺失值
exp_with_na[np.random.choice(10000, 500, replace=False)] = np.nan
act_with_na[np.random.choice(10000, 500, replace=False)] = np.nan

psi_drop = calculate_psi_robust(exp_with_na, act_with_na, missing_strategy='drop')
psi_separate = calculate_psi_robust(exp_with_na, act_with_na, missing_strategy='separate_bin')
print(f"剔除缺失值 PSI: {psi_drop:.4f}")
print(f"单独成箱 PSI: {psi_separate:.4f}")

3.5 实战案例：信贷评分卡 PSI 监控

以下是一个完整的信贷评分卡 PSI 监控案例，展示从数据加载到报警的全流程。

# 模拟信贷评分数据
np.random.seed(42)
train_scores = np.random.normal(650, 80, 50000)  # 训练集评分
prod_scores = np.random.normal(670, 95, 50000)   # 生产环境评分

# 计算PSI
psi_value, psi_detail = calculate_psi(train_scores, prod_scores, bins=10)

# 输出监控报告
print("="*50)
print("信贷评分卡 PSI 监控报告")
print("="*50)
print(f"PSI值: {psi_value:.4f}")

if psi_value < 0.1:
    status = "✅ 稳定"
elif psi_value < 0.25:
    status = "⚠️ 需关注"
else:
    status = "🚨 显著漂移"

print(f"状态: {status}")
print(f"\n分箱详情（前5箱）:")
print(psi_detail[['bin_start', 'bin_end', 'expected_pct', 'actual_pct', 'psi_component']].head())

3.6 性能优化与最佳实践

在大规模数据场景下，PSI 计算需要关注性能。以下是一些优化建议：

• 向量化计算：使用 numpy 数组操作替代 Python 循环，可提升 10-100 倍性能。
• 分箱数量选择：推荐 10-20 箱。箱数过少会丢失细节，箱数过多会导致部分箱占比为 0。
• 平滑处理：对零占比使用 1e-6 替代，避免除零错误和 ln(0) 问题。
• 缓存分箱边界：在连续监控场景中，固定分箱边界可确保不同时间点的 PSI 可比。