[人工智能]评分卡：WOE、IV、PSI计算及ROC和KS曲线

公式定义和原理解释见：

风控模型—WOE与IV指标的深入理解应用 - 知乎

风控模型—群体稳定性指标(PSI)深入理解应用 - 知乎

1、WOE和IV

延伸：分箱后求 WOE 和 IV

1. WOE describes the?relationship?between a predictive variable and a binary target variable.
2. IV measures the?strength?of that relationship.

# Autor chenfeng
#!/usr/bin/env Python
# coding=utf-8

'''
学习demo：
皮尔逊相关系数，woe计算
'''

import pandas as pd
import numpy as np
from scipy import stats


# 从正太分布种抽取20个随机数
x = np.array(np.random.randn((20)))
print(type(x))
X = pd.DataFrame(x, columns=['x'])
# 生成范围在[0, 2)的20个整数
y = np.random.randint(0,2,size=[20])
Y = pd.DataFrame(y, columns=['y'])
print(type(Y))   # dataframe类型

r = 0
# n分箱个数
n=5

## 有sum和count函数，Y的类型必须是DataFrame
good = Y.sum()
bad = Y.count() - good
print(type(good))  # Series类型

while np.abs(r) < 1:
    # pandas中使用qcut()，边界易出现重复值，如果为了删除重复值设置 duplicates=‘drop'，则易出现于分片个数少于指定个数的问题
    # 这里组合DataFrame需要的X、Y的类型应该是array，不能是DataFrame。可以有2种写法：
    # 方式一：
    # d1 = pd.DataFrame({"X": X['x'], "Y": Y['y'], "Bucket": pd.qcut(X['x'], n, duplicates="drop")})
    # 方式二：
    d1 = pd.DataFrame({"X": x, "Y": y, "Bucket": pd.qcut(X['x'], n, duplicates="drop")})
    # as_index=True会记录每个分箱包含的数据index
    d2 = d1.groupby('Bucket', as_index=True)
    '''
    计算每个分箱的X列平均值d2.mean().X 与每个分箱的Y列平均值d2.mean().Y 的spearman秩相关性（pearson相关系数只可用来计算连续数据、正太分布的线性相关性）。
    1、返回值r-相关系数，p-不相关系数，p-values并不完全可靠，但对于大于500左右的数据集可能是合理的。
    2、和其他相关系数一样，r和p在-1和+1之间变化，其中0表示无相关。
      -1或+1的相关关系意味着确切的单调关系。正相关表明，随着x的增加，y也随之增加。负相关性表示随着x增加，y减小。
    '''
    # d2.mean()取每个分箱的平均值
    r, p = stats.spearmanr(d2.mean().X, d2.mean().Y)
    n = n - 1
    print('n={}, r={}，p={}'.format(n, r, p))

print("=" * 60)
'''
d3 = pd.DataFrame(d2.X.min(), columns=['min'])得到的是一个空的dataframe
'''
d3 = pd.DataFrame(d2.X.min(), columns=['min'])
'''
d3：记录每个分箱的min、max、sum、total、rate=sum/total、woe
'''
d3['min'] = d2.min().X
d3['max'] = d2.max().X
d3['sum'] = d2.sum().Y
d3['total'] = d2.count().Y
d3['rate'] = d2.mean().Y

'''
1、Dataframe类型和Series类型做加减乘除计算时，Series的行索引对应Dataframe的列索引，如果两者没有匹配的索引，则返回Nan。
这里good和bad的行索引都是y，但是d3是没有y列的，而且这里其实只是想除以一个实数，所以直接取Series的values就行了。
2、d3['sum']/(d3['total']-d3['sum']) = d3['rate']/(1 - d3['rate'])
'''
d3['woe'] = np.log((d3['rate']/(1 - d3['rate']))/(good / bad).values)

d3['goodrate'] = d3['sum']/good.values
d3['badrate'] = (d3['total']-d3['sum'])/bad.values
d3['iv'] = ((d3['goodrate']-d3['badrate'])*d3['woe']).sum()
'''
1、sort_values(axis=0, by=‘min’)表示沿第一维的方向，对min列排序，默认是递增。axis=1表示沿第二维排序。
   sort_index()只能操作索引index，不能操作其他行或列
2、reset_index(drop=True) 重置索引，drop=true表示删除旧索引，重置后的索引默认从0开始的递增数值。
'''
d4 = (d3.sort_values(axis=0, by='min')).reset_index(drop=True)
print(d4)

?2、稳定性PSI

# Autor chenfeng
#!/usr/bin/env Python
# coding=utf-8

import math
import numpy as np
import pandas as pd


def calculate_psi(base_list, test_list, bins=20, min_sample=10):
    try:
        base_df = pd.DataFrame(base_list, columns=['score'])
        test_df = pd.DataFrame(test_list, columns=['score'])

        #---------------------------- 1.去除缺失值后，统计两个分布的样本量-----------------------------------------------
        base_notnull_cnt = len(list(base_df['score'].dropna()))
        test_notnull_cnt = len(list(test_df['score'].dropna()))

        # 空分箱
        base_null_cnt = len(base_df) - base_notnull_cnt
        test_null_cnt = len(test_df) - test_notnull_cnt

        #---------------------------- 2.最小分箱数：这里貌似是等距分箱？？每个分箱的距离是1/bin_num------------------------
        q_list = []
        if type(bins) == int:
            bin_num = min(bins, int(base_notnull_cnt / min_sample))
            q_list = [x / bin_num for x in range(1, bin_num)]
            break_list = []
            for q in q_list:
                '''
                df['score'].quantile(q=fraction, interpolation=linear)返回df['score']列的数据在q处的分位数值。
                其中，插值方式interpolation有以下几种：
                j-df['score']列中的最大值
                i-df['score']列中的最小值
                fraction-插值分位，当fraction=0.5, interpolation=linear时就是求中位数。
                * linear: `i + (j - i) * fraction`, where `fraction` is the
                  fractional part of the index surrounded by `i` and `j`.
                * lower: `i`.
                * higher: `j`.
                * nearest: `i` or `j` whichever is nearest.
                * midpoint: (`i` + `j`) / 2.
                '''
                bk = base_df['score'].quantile(q)
                break_list.append(bk)
            break_list = sorted(list(set(break_list)))  # 去重复后排序
            score_bin_list = [-np.inf] + break_list + [np.inf]
        else:
            score_bin_list = bins

        #----------------------------3.统计各分箱内的样本量----------------------------------------------------------------
        base_cnt_list = [base_null_cnt]
        test_cnt_list = [test_null_cnt]
        bucket_list = ["MISSING"]
        for i in range(len(score_bin_list) - 1):
            '''
            round()方法返回：按要求的小数位个数求四舍五入的结果。
            实际使用中发现round函数并不总是如上所说的四舍五入，这是因为该函数对于返回的浮点数并不是按照四舍五入的规则来计算，而会收到计算机表示精度的影响。
            当参数n不存在时，round()函数的输出为整数。
            当参数n存在时，即使为0，round()函数的输出也会是一个浮点数。
             此外，n的值可以是负数，表示在整数位部分四舍五入，但结果仍是浮点数。
            '''
            left = round(score_bin_list[i + 0], 4)
            right = round(score_bin_list[i + 1], 4)
            # 构建分箱格式，返回左开又闭的形式，如（0.1，0.2]
            bucket_list.append("(" + str(left) + ',' + str(right) + ']')

            base_cnt = base_df[(base_df.score > left) & (base_df.score <= right)].shape[0]
            base_cnt_list.append(base_cnt)

            test_cnt = test_df[(test_df.score > left) & (test_df.score <= right)].shape[0]
            test_cnt_list.append(test_cnt)

        #------------------------4.汇总统计结果:求每个分箱的占比率------------------------------------------------------------------------
        stat_df = pd.DataFrame({"bucket": bucket_list, "base_cnt": base_cnt_list, "test_cnt": test_cnt_list})
        stat_df['base_dist'] = stat_df['base_cnt'] / len(base_df)
        stat_df['test_dist'] = stat_df['test_cnt'] / len(test_df)

        def sub_psi(row):
            #------------------------ 5.计算每个分箱的不稳定性----------------------------------------------------------------------------
            base_list = row['base_dist']
            test_dist = row['test_dist']
            # 处理某分箱内样本量为0的情况
            if base_list == 0 and test_dist == 0:
                return 0
            elif base_list == 0 and test_dist > 0:
                base_list = 1 / base_notnull_cnt
            elif base_list > 0 and test_dist == 0:
                test_dist = 1 / test_notnull_cnt

            return (test_dist - base_list) * np.log(test_dist / base_list)

        '''
        对stat_df的每行应用sub_psi(...)函数
        '''
        stat_df['psi'] = stat_df.apply(lambda row: sub_psi(row), axis=1)
        '''
        从stat_df中取['bucket', 'base_cnt', 'base_dist', 'test_cnt', 'test_dist', 'psi']这些列的数据
        '''
        stat_df = stat_df[['bucket', 'base_cnt', 'base_dist', 'test_cnt', 'test_dist', 'psi']]
        # ------------------------6、PSI = 所有分箱的不稳定性之和------------------------------------------------------------------------
        psi = stat_df['psi'].sum()

    except:
        print('error!!!')
        psi = np.nan
        stat_df = None
    return psi, stat_df


if __name__=='__main__':
    # 从正太分布种抽取20个随机数
    base_x = np.array(np.random.randn(100))
    base_list = pd.DataFrame({"score": base_x})

    test_x = np.array(np.random.randn(20))
    test_list = pd.DataFrame(test_x, columns=['score'])

    psi, stat_df = calculate_psi(base_list=base_list,
                                 test_list=test_list,
                                 bins=20, min_sample=10)
    print("psi:", psi)
    print(stat_df)

3、ROC和KS曲线

roc曲线涉：
x轴--假正率(negative的召回率)fpr
y轴--真正率(positive的召回率)tpr

常用：sklearn.metrics.roc_curve使用简要说明_sun91019718的博客-CSDN博客_sklearn.metrics.roc_curve

tpr,fpr,threshold = sklearn.metrics.roc_curve(y_label,y_pred,pos_label=1)
备注：pos_label=positive label,当y_label是{0，1}或{-1，1}内的值，可以不写；否则必须明确。

# Autor chenfeng
#!/usr/bin/env Python
# coding=utf-8
import numpy as np
from sklearn import metrics
import matplotlib.pyplot as plt



def plot_roc(y_label,y_pred):
    """
    绘制roc曲线
    param:
        y_label -- 真实的y值 list/array
        y_pred -- 预测的y值 list/array
    return:
        roc曲线
    """
    '''
    1、tpr（真正率）和fpr（假正率）
    2、当y的值不在{0, 1} or {-1, 1}范围内时，需要设置positive label参数：pos_label
    '''
    tpr,fpr,threshold = metrics.roc_curve(y_label, y_pred, pos_label=2)
    print("真正率tpr：", tpr)
    print("假正率fpr：", fpr)
    print("threshold：", threshold)

    '''
    AUC是roc曲线下方的面积
    '''
    AUC = metrics.roc_auc_score(y_label,y_pred)

    fig = plt.figure(figsize=(6,4))
    ax = fig.add_subplot(1,1,1)
    ax.plot(fpr, tpr, color='blue', label='AUC=%.3f'%AUC)
    ax.plot([0,1],[0,1],'r--')
    ax.set_ylim(0,1)
    ax.set_xlim(0,1)
    ax.set_title('ROC')
    ax.legend(loc='best')
    plt.show()


def plot_model_ks(y_label, y_pred):
    """
    绘制ks曲线
    param:
        y_label -- 真实的y值 list/array
        y_pred -- 预测的y值 list/array
    return:
        ks曲线
    """
    pred_list = list(y_pred)
    label_list = list(y_label)
    total_bad = sum(label_list)
    total_good = len(label_list) - total_bad
    items = sorted(zip(pred_list, label_list), key=lambda x: x[0])
    step = (max(pred_list) - min(pred_list)) / 200

    pred_bin = []
    good_rates = []
    bad_rates = []
    ks_list = []
    for i in range(1, 201):
        idx = min(pred_list) + i * step
        pred_bin.append(idx)
        label_bin = [x[1] for x in items if x[0] < idx]
        bad_num = sum(label_bin)
        good_num = len(label_bin) - bad_num
        goodrate = good_num / total_good
        badrate = bad_num / total_bad
        ks = abs(goodrate - badrate)
        good_rates.append(goodrate)
        bad_rates.append(badrate)
        ks_list.append(ks)

    fig = plt.figure(figsize=(6, 4))
    ax = fig.add_subplot(1, 1, 1)
    ax.plot(pred_bin, good_rates, color='green', label='good_rate')
    ax.plot(pred_bin, bad_rates, color='red', label='bad_rate')
    ax.plot(pred_bin, ks_list, color='blue', label='good-bad')
    ax.set_title('KS:{:.3f}'.format(max(ks_list)))
    ax.legend(loc='best')
    plt.show()



y_label = np.array([1, 1, 2, 2])
y_pred = np.array([0.1, 0.4, 0.35, 0.8])
plot_roc(y_label, y_pred)
plot_model_ks(y_label, y_pred)