逻辑回归 - 评分卡项目(五)

2,090 阅读25分钟

根据菜菜的课程进行整理,方便记忆理解

代码位置如下:

用逻辑回归制作评分卡

评分卡是一种以分数形式来衡量一个客户的信用风险大小的手段

  • 它衡量向别人借钱的人(受信人,需要融资的公司)不能如期履行合同中的还本付息责任
  • 借钱给别人的人(授信人,银行等金融机构)造成经济损失的可能性

一般来说,评分卡打出的分数越高,客户的信用越好,风险越小

一个完整的模型开发,需要有以下流程:

image.png

以个人消费类贷款数据,来为大家简单介绍A卡的建模和制作流程,核心会在”数据清洗“和“模型开发”上。

导库,获取数据

%matplotlib inline
import numpy as np
import pandas as pd
from sklearn.linear_model import LogisticRegression as LR
 
#其实日常在导库的时候,并不是一次性能够知道我们要用的所有库的。通常都是在建模过程中逐渐导入需要的库。

在银行系统中,这个数据通常使来自于其他部门的同事的收集,因此千万别忘记抓住给你数据的人,问问她/他各个项都是什么含义。通常来说,当特征非常多的时候(比如几百个),都会有一个附带的excel或pdf文档给到你,备注了各个特征都是什么含义。这种情况下其实要一个个去看还是非常困难,所以如果特征很多,建议先做降维.

探索数据与数据预处理

data = pd.read_csv(r".\rankingcard.csv",index_col=0)

#观察数据结构
data.shape   #(150000, 11)

data.info()

"""
<class 'pandas.core.frame.DataFrame'>
Int64Index: 150000 entries, 1 to 150000
Data columns (total 11 columns):
 #   Column                                Non-Null Count   Dtype  
---  ------                                --------------   -----  
 0   SeriousDlqin2yrs                      150000 non-null  int64  
 1   RevolvingUtilizationOfUnsecuredLines  150000 non-null  float64
 2   age                                   150000 non-null  int64  
 3   NumberOfTime30-59DaysPastDueNotWorse  150000 non-null  int64  
 4   DebtRatio                             150000 non-null  float64
 5   MonthlyIncome                         120269 non-null  float64
 6   NumberOfOpenCreditLinesAndLoans       150000 non-null  int64  
 7   NumberOfTimes90DaysLate               150000 non-null  int64  
 8   NumberRealEstateLoansOrLines          150000 non-null  int64  
 9   NumberOfTime60-89DaysPastDueNotWorse  150000 non-null  int64  
 10  NumberOfDependents                    146076 non-null  float64
dtypes: float64(4), int64(7)
memory usage: 13.7 MB
"""

#观察数据类型
data  #注意可以看到第一列为标签,剩下的10列为特征

了解样本总体的大概情况,比如查看缺失值,量纲是否统一,是否需要做哑变量等等。

特征/标签含义
SeriousDlqin2yrs出现 90 天或更长时间的逾期行为(即定义好坏客户)
RevolvingUtilizationOfUnsecuredLines贷款以及信用卡可用额度与总额度比例
age借款人借款年龄
NumberOfTime30-59DaysPastDueNotWorse过去两年内出现35-59天逾期但是没有发展得更坏的次数
DebtRatio每月偿还债务,赡养费,生活费用除以月总收入
MonthlyIncome月收入
NumberOfOpenCreditLinesAndLoans开放式贷款和信贷数量
NumberOfTimes90DaysLate过去两年内出现90天逾期或更坏的次数
NumberRealEstateLoansOrLines抵押贷款和房地产贷款数量,包括房屋净值信贷额度
NumberOfTime60-89DaysPastDueNotWorse过去两年内出现60-89天逾期但是没有发展得更坏的次数
NumberOfDependents家庭中不包括自身的家属人数(配偶,子女等)
去除重复值

现实数据,尤其是银行业数据,可能会存在的一个问题就是样本重复,即有超过一行的样本所显示的所有特征都一样。有时候可能时人为输入重复,有时候可能是系统录入重复,总而言之我们必须对数据进行去重处理。可能有人会说,难道不可能出现说两个样本的特征就是一模一样,但他们是两个样本吗?比如,两个人,一模一样的名字,年龄,性别,学历,工资……当特征量很少的时候,这的确是有可能的,但一些指标,比如说家属人数,月收入,已借有的房地产贷款数量等等,几乎不可能都出现一样。尤其是银行业数据经常是几百个特征,所有特征都一样的可能性是微乎其微的。即便真的出现了如此极端的情况,我们也可以当作是少量信息损失,将这条记录当作重复值除去。

#去除重复值
data.drop_duplicates(inplace=True)#inplace=True表示替换原数据

#删除之后千万不要忘记,恢复索引
data.index = range(len(data))
填补缺失值
data.isnull().sum()/data.shape[0]     #得到缺失值的比例
# data.isnull().mean()                #上一行代码的另一种形式书写

"""
SeriousDlqin2yrs                        0.000000
RevolvingUtilizationOfUnsecuredLines    0.000000
age                                     0.000000
NumberOfTime30-59DaysPastDueNotWorse    0.000000
DebtRatio                               0.000000
MonthlyIncome                           0.195601
NumberOfOpenCreditLinesAndLoans         0.000000
NumberOfTimes90DaysLate                 0.000000
NumberRealEstateLoansOrLines            0.000000
NumberOfTime60-89DaysPastDueNotWorse    0.000000
NumberOfDependents                      0.025624
dtype: float64
"""

第二个要面临的问题,就是缺失值。在这里我们需要填补的特征是“收入”和“家属人数”。“家属人数”缺失很少,仅缺失了大约2.5%,可以考虑直接删除,或者使用均值来填补。“收入”缺失了几乎20%,并且我们知道,“收入”必然是一个对信用评分来说很重要的因素,因此这个特征必须要进行填补。在这里,我们使用均值填补“家属人数”。

data["NumberOfDependents"].fillna(int(data["NumberOfDependents"].mean()),inplace=True)
# 这里用均值填补家庭人数这一项 
# 如果你选择的是删除那些缺失了2.5%的特征,千万记得恢复索引哟~
data.isnull().sum()/data.shape[0]

"""
SeriousDlqin2yrs                        0.000000
RevolvingUtilizationOfUnsecuredLines    0.000000
age                                     0.000000
NumberOfTime30-59DaysPastDueNotWorse    0.000000
DebtRatio                               0.000000
MonthlyIncome                           0.195601
NumberOfOpenCreditLinesAndLoans         0.000000
NumberOfTimes90DaysLate                 0.000000
NumberRealEstateLoansOrLines            0.000000
NumberOfTime60-89DaysPastDueNotWorse    0.000000
NumberOfDependents                      0.000000
dtype: float64
"""

那字段"收入"怎么办呢?对于银行数据来说,我们甚至可以有这样的推断:

一个来借钱的人应该是会知道,“高收入”或者“稳定收入”于他/她自己而言会是申请贷款过程中的一个助力,因此如果收入稳定良好的人,肯定会倾向于写上自己的收入情况,那么这些“收入”栏缺失的人,更可能是收入状况不稳定或收入比较低的人。

基于这种判断,我们可以用比如说,四分位数来填补缺失值,把所有收入为空的客户都当成是低收入人群。当然了,也有可能这些缺失是银行数据收集过程中的失误,我们并无法判断为什么收入栏会有缺失,所以我们的推断也有可能是不正确的。具体采用什么样的手段填补缺失值,要和业务人员去沟通,观察缺失值是如何产生的。在这里,我们使用随机森林填补“收入”。

之前我们所做的随机森林填补缺失值的案例中,我们面临整个数据集中多个特征都有缺失的情况,因此要先对特征排序,遍历所有特征来进行填补。这次我们只需要填补“收入”一个特征,就无需循环那么麻烦了,可以直接对这一列进行填补。我们来写一个能够填补任何列的函数:

def fill_missing_rf(X,y,to_fill):
    """
    使用随机森林填补一个特征的缺失值的函数

    参数:
    X:要填补的特征矩阵
    y:完整的,没有缺失值的标签
    to_fill:字符串,要填补的那一列的名称
    """

    #构建我们的新特征矩阵和新标签
    df = X.copy()
    fill = df.loc[:,to_fill]
    df = pd.concat([df.loc[:,df.columns != to_fill],pd.DataFrame(y)],axis=1)

    # 找出我们的训练集和测试集
    Ytrain = fill[fill.notnull()]
    Ytest = fill[fill.isnull()]
    Xtrain = df.iloc[Ytrain.index,:]
    Xtest = df.iloc[Ytest.index,:]

    #用随机森林回归来填补缺失值
    from sklearn.ensemble import RandomForestRegressor as rfr
    rfr = rfr(n_estimators=100)
    rfr = rfr.fit(Xtrain, Ytrain)
    Ypredict = rfr.predict(Xtest)

    return Ypredict

接下来,我们来创造函数需要的参数,将参数导入函数,产出结果:

X = data.iloc[:,1:]
y = data["SeriousDlqin2yrs"]  # y = data.iloc[:,0]
X.shape                       # (149391, 10)

#=====[TIME WARNING:1 min]=====#
y_pred = fill_missing_rf(X,y,"MonthlyIncome")

#注意可以通过以下代码检验数据是否数量相同
y_pred.shape ==  data.loc[data.loc[:,"MonthlyIncome"].isnull(),"MonthlyIncome"].shape
# True

#确认我们的结果合理之后,我们就可以将数据覆盖了
data.loc[data.loc[:,"MonthlyIncome"].isnull(),"MonthlyIncome"] = y_pred
描述性统计处理异常值

现实数据永远都会有一些异常值,首先我们要去把他们捕捉出来,然后观察他们的性质。注意,我们并不是要排除掉所有异常值,相反很多时候,异常值是我们的重点研究对象,比如说,双十一中购买量超高的品牌,或课堂上让很多学生都兴奋的课题,这些是我们要重点研究观察的。

日常处理异常值,我们使用箱线图或者**3δ3\delta法则来找到异常值(千万不要说依赖于眼睛看,我们是数据挖掘工程师,除了业务理解,我们还要有方法)。但在银行数据中,我们希望排除的“异常值”不是一些超高或超低的数字,而是一些不符合常理**的数据:收入不能为负数,但是一个超高水平的收入却是合理的,可以存在的。

所以在银行业中,我们往往就使用普通的描述性统计来观察数据的异常与否与数据的分布情况。注意,这种方法只能在特征量有限的情况下进行,如果有几百个特征又无法成功降维或特征选择不管用,那还是用比较好。

#描述性统计
data.describe()

image.png

data.describe([0.01,0.1,0.25,.5,.75,.9,.99]).T

image.png

# 异常值也被我们观察到,年龄的最小值居然有0,这不符合银行的业务需求,即便是儿童账户也要至少8岁,我们可以
# 查看一下年龄为0的人有多少
(data["age"] == 0).sum()
# 1

#发现只有一个人年龄为0,可以判断这肯定是录入失误造成的,可以当成是缺失值来处理,直接删除掉这个样本
data = data[data["age"] != 0]

有三个指标看起来很奇怪

"""
另外,有三个指标看起来很奇怪:
 
"NumberOfTime30-59DaysPastDueNotWorse"
"NumberOfTime60-89DaysPastDueNotWorse"
"NumberOfTimes90DaysLate"
 
这三个指标分别是“过去两年内出现35-59天逾期但是没有发展的更坏的次数”,“过去两年内出现60-89天逾期但是没
有发展的更坏的次数”,“过去两年内出现90天逾期的次数”。这三个指标,在99%的分布的时候依然是2,最大值却是
98,看起来非常奇怪。一个人在过去两年内逾期35~59天98次,一年6个60天,两年内逾期98次这是怎么算出来的?
 
我们可以去咨询业务人员,请教他们这个逾期次数是如何计算的。如果这个指标是正常的,那这些两年内逾期了98次的
客户,应该都是坏客户。在我们无法询问他们情况下,我们查看一下有多少个样本存在这种异常:
 
"""
# 我们的常识是违约98次一定是一个坏客户,但是查看标签我们发现这个数据中有违约的人也有没有违约的人
# 并且我们通过下面的value_counts(),可以看到98次的是220,96次的是5个,0 - 17的人很少,且呈现递减的趋势,那么98/96一定都是异常值
data[data.loc[:,"NumberOfTimes90DaysLate"] > 90]

image.png

data[data.loc[:,"NumberOfTimes90DaysLate"] > 90].count()

"""
SeriousDlqin2yrs                        225
RevolvingUtilizationOfUnsecuredLines    225
age                                     225
NumberOfTime30-59DaysPastDueNotWorse    225
DebtRatio                               225
MonthlyIncome                           225
NumberOfOpenCreditLinesAndLoans         225
NumberOfTimes90DaysLate                 225
NumberRealEstateLoansOrLines            225
NumberOfTime60-89DaysPastDueNotWorse    225
NumberOfDependents                      225
dtype: int64
"""

data.loc[:,"NumberOfTimes90DaysLate"].value_counts()

# 有225个样本存在这样的情况,并且这些样本,我们观察一下,标签并不都是1,他们并不都是坏客户。因此,我们基
# 本可以判断,这些样本是某种异常,应该把它们删除。
 
data = data[data.loc[:,"NumberOfTimes90DaysLate"] < 90]

#一定要恢复索引
data.index = range(data.shape[0])
data.info()

为什么不统一量纲,也不标准化数据分布

在描述性统计结果中,我们可以观察到数据量纲明显不统一,而且存在一部分极偏的分布,虽然逻辑回归对于数据没有分布要求,但是我们知道如果数据服从正态分布的话梯度下降可以收敛得更快。但在这里,我们不对数据进行标准化处理,也不进行量纲统一,为什么?

无论算法有什么样的规定,无论统计学中有什么样的要求,我们的最终目的都是要为业务服务。现在我们要制作评分卡,评分卡是要给业务人员们使用的基于新客户填写的各种信息为客户打分的一张卡片,而为了制作这张卡片,我们需要对我们的数据进行一个“分档”,比如说,年龄20 ~ 30岁为一档,年龄30 ~ 50岁为一档,月收入1W以上为一档,5000~1W为一档,每档的分数不同。

一旦我们将数据统一量纲,或者标准化了之后,数据大小和范围都会改变,统计结果是漂亮了,但是对于业务人员来说,他们完全无法理解,标准化后的年龄在0.00328 ~ 0.00467之间为一档是什么含义。并且,新客户填写的信息,天生就是量纲不统一的,我们的确可以将所有的信息录入之后,统一进行标准化,然后导入算法计算,但是最终落到业务人员手上去判断的时候,他们会完全不理解为什么录入的信息变成了一串统计上很美但实际上根本看不懂的数字。由于业务要求,在制作评分卡的时候,我们要尽量保持数据的原貌,年龄就是8 ~ 110的数字,收入就是大于0,最大值可以无限的数字,即便量纲不统一,我们也不对数据进行标准化处理。

样本不均衡问题
#探索标签的分布
X = data.iloc[:,1:]
y = data.iloc[:,0]
 
y.value_counts()#查看每一类别值得数据量,查看样本是否均衡

n_sample = X.shape[0]

n_1_sample = y.value_counts()[1]
n_0_sample = y.value_counts()[0]
print('样本个数:{}; 1占{:.2%}; 0占{:.2%}'.format(n_sample,n_1_sample/n_sample,n_0_sample/n_sample))
# 样本个数:149165; 1占6.62%; 0占93.38%

可以看出,样本严重不均衡。虽然大家都在努力防范信用风险,但实际违约的人并不多。并且,银行并不会真的一棒子打死所有会违约的人,很多人是会还钱的,只是忘记了还款日,很多人是不愿意欠人钱的,但是当时真的很困难,资金周转不过来,所以发生逾期,但一旦他有了钱,他就会把钱换上。对于银行来说,只要你最后能够把钱还上,我都愿意借钱给你,因为我借给你就有收入(利息)。所以,对于银行来说,真正想要被判别出来的其实是”恶意违约“的人,而这部分人数非常非常少,样本就会不均衡。这一直是银行业建模的一个痛点:我们永远希望捕捉少数类。

之前提到过,逻辑回归中使用最多的是上采样方法来平衡样本。

#如果报错,就在prompt安装:pip install imblearn
import imblearn
#imblearn是专门用来处理不平衡数据集的库,在处理样本不均衡问题中性能高过sklearn很多
#imblearn里面也是一个个的类,也需要进行实例化,fit拟合,和sklearn用法相似
 
from imblearn.over_sampling import SMOTE
 
sm = SMOTE(random_state=42) #实例化
X,y = sm.fit_resample(X,y)

n_sample_ = X.shape[0]  # 278584

pd.Series(y).value_counts()
n_1_sample = pd.Series(y).value_counts()[1]
n_0_sample = pd.Series(y).value_counts()[0]
print('样本个数:{}; 1占{:.2%}; 0占{:.2%}'.format(n_sample_,n_1_sample/n_sample_,n_0_sample/n_sample_))
# 样本个数:278584; 1占50.00%; 0占50.00%

如此,我们就实现了样本平衡,样本量也增加了

分训练集和测试集
# 对中间数据进行保存
from sklearn.model_selection import train_test_split
X = pd.DataFrame(X)
y = pd.DataFrame(y)
 
X_train, X_vali, Y_train, Y_vali = train_test_split(X,y,test_size=0.3,random_state=420)
model_data = pd.concat([Y_train, X_train], axis=1)#训练数据构建模型
model_data.index = range(model_data.shape[0])
model_data.columns = data.columns
 
vali_data = pd.concat([Y_vali, X_vali], axis=1)#验证集
vali_data.index = range(vali_data.shape[0])
vali_data.columns = data.columns
 
model_data.to_csv(r".\model_data.csv")#训练数据
vali_data.to_csv(r".\vali_data.csv")#验证数据

分箱

前面提到过,我们要制作评分卡,是要给各个特征进行分档,以便业务人员能够根据新客户填写的信息为客户打分。因此在评分卡制作过程中,一个重要的步骤就是分箱。可以说,分箱是评分卡最难,也是最核心的思路,分箱的本质,其实就是离散化连续变量,好让拥有不同属性的人被分成不同的类别(打上不同的分数),其实本质比较类似于聚类。那我们在分箱中要回答几个问题:

要分多少个箱子才合适

最开始我们并不知道,但是既然是将连续型变量离散化,想也知道箱子个数必然不能太多,最好控制在十个以下。而用来制作评分卡,最好能在4 ~ 5个为最佳。我们知道,离散化连续变量必然伴随着信息的损失,并且箱子越少,信息损失越大。为了衡量特征上的信息量以及特征对预测函数的贡献,银行业定义了概念Information value(IV):

image.png

其中N是这个特征上箱子的个数,i代表每个箱子,good% 是这个箱内的优质客户(标签为0的客户)占整个特征中所有优质客户的比例, bad%是这个箱子里的坏客户(就是那些会违约,标签为1的那些客户)占整个特征中所有坏客户的比例,而WOEiWOE_i则写作:

image.png

这是我们在银行业中用来衡量违约概率的指标,中文叫做证据权重(weight of Evidence),本质其实就是优质客户比上坏客户的比例的对数。WOE是对一个箱子来说的,WOE越大,代表了这个箱子里的优质客户越多。而IV是对整个特征来说的,IV代表的意义是我们特征上的信息量以及这个特征对模型的贡献,由下表来控制:

IV特征对预测函数的贡献
< 0.03特征几乎不带有效信息,对模型没有贡献,这种特征可以被删除
0.03 ~ 0.09有效信息很少,对模型的贡献度低
0.1 ~ 0.29有效信息一般,对模型的贡献度中等
0.3 ~ 0.49有效信息较多,对模型的贡献度较高
>=0.5有效信息非常多,对模型的贡献超高并且可疑

可见,IV并非越大越好,我们想要找到IV的大小和箱子个数的平衡点。箱子越多,IV必然越小,因为信息损失会非常多,所以,我们会对特征进行分箱,然后计算每个特征在每个箱子数目下的WOE值,利用IV值的曲线,找出合适的分箱个数.

分箱要达成什么样的效果

我们希望不同属性的人有不同的分数,因此我们希望在同一个箱子内的人的属性是尽量相似的,而不同箱子的人的属性是尽量不同的,即业界常说的”组间差异大,组内差异小“。对于评分卡来说,就是说我们希望一个箱子内的人违约概率是类似的,而不同箱子的人的违约概率差距很大,即WOE差距要大,并且每个箱子中坏客户所占的比重(bad%)也要不同。那我们,可以使用卡方检验来对比两个箱子之间的相似性,如果两个箱子之间卡方检验的P值很大,则说明他们非常相似,那我们就可以将这两个箱子合并为一个箱子

基于这样的思想,我们总结出我们对一个特征进行分箱的步骤:

  • 我们首先把连续型变量分成一组数量较多的分类型变量,比如,将几万个样本分成100组,或50组
  • 确保每一组中都要包含两种类别的样本,否则IV值会无法计算
  • 我们对相邻的组进行卡方检验,卡方检验的P值很大的组进行合并,直到数据中的组数小于设定的N箱为止
  • 我们让一个特征分别分成[2,3,4.....20]箱,观察每个分箱个数下的IV值如何变化,找出最适合的分箱个数
  • 分箱完毕后,我们计算每个箱的WOE值,bad% ,观察分箱效果

这些步骤都完成后,我们可以对各个特征都进行分箱,然后观察每个特征的IV值,以此来挑选特征。

接下来,我们就以"age"为例子,来看看分箱如何完成。注意,分箱代码的版权属于Hsiaofei Tsien,我已获得授权在这门课中使用和讲解他的代码。

等频分箱
#按照等频对需要分箱的列进行分箱

#“age”为例子
model_data["qcut"], updown = pd.qcut(model_data["age"], retbins=True, q=20)#等频分箱

model_data["qcut"]

"""
0         (52.0, 54.0]
1         (61.0, 64.0]
2         (36.0, 39.0]
3         (68.0, 74.0]
4         (52.0, 54.0]
              ...     
195003    (31.0, 34.0]
195004    (48.0, 50.0]
195005    (45.0, 46.0]
195006    (61.0, 64.0]
195007    (52.0, 54.0]
Name: qcut, Length: 195008, dtype: category
Categories (20, interval[float64]): [(20.999, 28.0] < (28.0, 31.0] < (31.0, 34.0] < (34.0, 36.0] ... (61.0, 64.0] < (64.0, 68.0] < (68.0, 74.0] < (74.0, 107.0]]
"""

# 在最后的位置拼接一列特征,所所属的箱体
model_data.head()

image.png

"""
pd.qcut,基于分位数的分箱函数,本质是将连续型变量离散化
只能够处理一维数据。返回箱子的上限和下限
参数q:要分箱的个数
参数retbins=True来要求同时返回结构为索引为样本索引,元素为分到的箱子的Series
现在返回两个值:每个样本属于哪个箱子,以及所有箱子的上限和下限
"""
#在这里时让model_data新添加一列叫做“分箱”,这一列其实就是每个样本所对应的箱子
model_data["qcut"].value_counts()

"""
(36.0, 39.0]      12701
(20.999, 28.0]    11788
(58.0, 61.0]      11347
(48.0, 50.0]      11093
(46.0, 48.0]      10978
(31.0, 34.0]      10809
(50.0, 52.0]      10553
(43.0, 45.0]      10390
(61.0, 64.0]      10186
(39.0, 41.0]       9798
(52.0, 54.0]       9720
(41.0, 43.0]       9666
(28.0, 31.0]       9528
(74.0, 107.0]      9143
(64.0, 68.0]       8928
(54.0, 56.0]       8688
(68.0, 74.0]       8665
(56.0, 58.0]       7883
(34.0, 36.0]       7467
(45.0, 46.0]       5677
Name: qcut, dtype: int64
"""

#所有箱子的上限和下限
updown

# array([ 21.,  28.,  31.,  34.,  36.,  39.,  41.,  43.,  45.,  46.,  48.,50.,  52.,  54.,  56.,  58.,  61.,  64.,  68.,  74., 107.])

model_data[model_data["SeriousDlqin2yrs"] == 0].groupby(by="qcut").count()

image.png

确保每个箱中都有0和1
# 统计每个分箱中0和1的数量
# 这里使用了数据透视表的功能groupby
coount_y0 = model_data[model_data["SeriousDlqin2yrs"] == 0].groupby(by="qcut").count()["SeriousDlqin2yrs"]

coount_y1 = model_data[model_data["SeriousDlqin2yrs"] == 1].groupby(by="qcut").count()["SeriousDlqin2yrs"]

#num_bins值分别为每个区间的上界,下界,0出现的次数,1出现的次数
num_bins = [*zip(updown,updown[1:],coount_y0,coount_y1)]
 
#注意zip会按照最短列来进行结合
num_bins
"""
[(21.0, 28.0, 4243, 7545),
 (28.0, 31.0, 3571, 5957),
 (31.0, 34.0, 4075, 6734),
 (34.0, 36.0, 2908, 4559),
 (36.0, 39.0, 5182, 7519),
 (39.0, 41.0, 3956, 5842),
 (41.0, 43.0, 4002, 5664),
 (43.0, 45.0, 4389, 6001),
 (45.0, 46.0, 2419, 3258),
 (46.0, 48.0, 4813, 6165),
 (48.0, 50.0, 4900, 6193),
 (50.0, 52.0, 4728, 5825),
 (52.0, 54.0, 4681, 5039),
 (54.0, 56.0, 4677, 4011),
 (56.0, 58.0, 4483, 3400),
 (58.0, 61.0, 6583, 4764),
 (61.0, 64.0, 6968, 3218),
 (64.0, 68.0, 6623, 2305),
 (68.0, 74.0, 6753, 1912),
 (74.0, 107.0, 7737, 1406)]
"""

# 确保每个箱子中都含正样本或负样本
for i in range(20):
    #如果第一个组没有包含正样本或负样本,向后合并
    if 0 in num_bins[0][2:]:
        num_bins[0:2] = [(
            num_bins[0][0],
            num_bins[1][1],
            num_bins[0][2]+num_bins[1][2],
            num_bins[0][3]+num_bins[1][3])]
        continue

    """
    合并了之后,第一行的组是否一定有两种样本了呢?不一定
    如果原本的第一组和第二组都没有包含正样本,或者都没有包含负样本,那即便合并之后,第一行的组也还是没有
    包含两种样本
    所以我们在每次合并完毕之后,还需要再检查,第一组是否已经包含了两种样本
    这里使用continue跳出了本次循环,开始下一次循环,所以回到了最开始的for i in range(20), 让i+1
    这就跳过了下面的代码,又从头开始检查,第一组是否包含了两种样本
    如果第一组中依然没有包含两种样本,则if通过,继续合并,每合并一次就会循环检查一次,最多合并20次
    如果第一组中已经包含两种样本,则if不通过,就开始执行下面的代码
    """
    #已经确认第一组中肯定包含两种样本了,如果其他组没有包含两种样本,就向前合并
    #此时的num_bins已经被上面的代码处理过,可能被合并过,也可能没有被合并
    #但无论如何,我们要在num_bins中遍历,所以写成in range(len(num_bins))
    for i in range(len(num_bins)):
        if 0 in num_bins[i][2:]:
            num_bins[i-1:i+1] = [(
                num_bins[i-1][0],
                num_bins[i][1],
                num_bins[i-1][2]+num_bins[i][2],
                num_bins[i-1][3]+num_bins[i][3])]
        break
        #如果对第一组和对后面所有组的判断中,都没有进入if去合并,则提前结束所有的循环
    else:
        break

    """
    这个break,只有在if被满足的条件下才会被触发
    也就是说,只有发生了合并,才会打断for i in range(len(num_bins))这个循环
    为什么要打断这个循环?因为我们是在range(len(num_bins))中遍历
    但合并发生后,len(num_bins)发生了改变,但循环却不会重新开始
    举个例子,本来num_bins是5组,for i in range(len(num_bins))在第一次运行的时候就等于for i in 
    range(5)
    range中输入的变量会被转换为数字,不会跟着num_bins的变化而变化,所以i会永远在[0,1,2,3,4]中遍历
    进行合并后,num_bins变成了4组,已经不存在=4的索引了,但i却依然会取到4,循环就会报错
    因此在这里,一旦if被触发,即一旦合并发生,我们就让循环被破坏,使用break跳出当前循环
    循环就会回到最开始的for i in range(20)中
    此时判断第一组是否有两种标签的代码不会被触发,但for i in range(len(num_bins))却会被重新运行
    这样就更新了i的取值,循环就不会报错了
    """

分箱实战

定义WOE和IV函数
#计算WOE和BAD RATE
#BAD RATE与bad%不是一个东西
#BAD RATE是一个箱中,坏的样本所占的比例 (bad/total)
#而bad%是一个箱中的坏样本占整个特征中的坏样本的比例
 
def get_woe(num_bins):
    # 通过 num_bins 数据计算 woe
    columns = ["min","max","count_0","count_1"]
    df = pd.DataFrame(num_bins,columns=columns)

    df["total"] = df.count_0 + df.count_1#一个箱子当中所有的样本数
    df["percentage"] = df.total / df.total.sum()#一个箱子里的样本数,占所有样本的比例
    df["bad_rate"] = df.count_1 / df.total#一个箱子坏样本的数量占一个箱子里边所有样本数的比例
    df["good%"] = df.count_0/df.count_0.sum()
    df["bad%"] = df.count_1/df.count_1.sum()
    df["woe"] = np.log(df["good%"] / df["bad%"])
    return df

#计算IV值
def get_iv(df):
    rate = df["good%"] - df["bad%"]
    iv = np.sum(rate * df.woe)
    return iv
卡方检验,合并箱体,画出IV曲线
num_bins_ = num_bins.copy()
 
import matplotlib.pyplot as plt
import scipy
 
IV = []
axisx = []
 
while len(num_bins_) > 2:#大于设置的最低分箱个数
    pvs = []
    #获取 num_bins_两两之间的卡方检验的置信度(或卡方值)
    for i in range(len(num_bins_)-1):
        x1 = num_bins_[i][2:]
        x2 = num_bins_[i+1][2: ]
        # 0 返回 chi2 值,1 返回 p 值。
        pv = scipy.stats.chi2_contingency([x1,x2])[1]#p值
        # chi2 = scipy.stats.chi2_contingency([x1,x2])[0]#计算卡方值
        pvs.append(pv)
        
    # 通过 p 值进行处理。合并 p 值最大的两组
    # 找到最大的chi2计算得到的pv值,我们通过下标可以找到对应的两组数据,相邻的,所以合并的时候,min = 第一个数据的小年龄
    # max = 第二个数据的大年龄,然后将0/1的数量进行相加即可
    i = pvs.index(max(pvs))
    num_bins_[i:i+2] = [(
            num_bins_[i][0],
            num_bins_[i+1][1],
            num_bins_[i][2]+num_bins_[i+1][2],
            num_bins_[i][3]+num_bins_[i+1][3])]
    
    bins_df = get_woe(num_bins_)
    axisx.append(len(num_bins_))
    IV.append(get_iv(bins_df))
    
plt.figure()
plt.plot(axisx,IV)
plt.xticks(axisx)
plt.xlabel("number of box")
plt.ylabel("IV")
plt.show()
# 选择转折点处,也就是下坠最快的折线点,所以这里对于age来说选择箱数为6

image.png

用最佳分箱个数分箱,并验证分箱结果
def get_bin(num_bins_,n):
    while len(num_bins_) > n:
        pvs = []
        for i in range(len(num_bins_)-1):
            x1 = num_bins_[i][2:]
            x2 = num_bins_[i+1][2:]
            pv = scipy.stats.chi2_contingency([x1,x2])[1]
            # chi2 = scipy.stats.chi2_contingency([x1,x2])[0]
            pvs.append(pv)

        i = pvs.index(max(pvs))
        num_bins_[i:i+2] = [(
                num_bins_[i][0],
                num_bins_[i+1][1],
                num_bins_[i][2]+num_bins_[i+1][2],
                num_bins_[i][3]+num_bins_[i+1][3])]
    return num_bins_
 
afterbins = get_bin(num_bins,6)
 
afterbins

"""
[(21.0, 36.0, 14797, 24795),
 (36.0, 54.0, 39070, 51506),
 (54.0, 61.0, 15743, 12175),
 (61.0, 64.0, 6968, 3218),
 (64.0, 74.0, 13376, 4217),
 (74.0, 107.0, 7737, 1406)]
"""

bins_df = get_woe(num_bins)
 
bins_df
#希望每组的bad_rate相差越大越好;
# woe差异越大越好,应该具有单调性,随着箱的增加,要么由正到负,要么由负到正,只能有一个转折过程;
# 如果woe值大小变化是有两个转折,比如呈现w型,证明分箱过程有问题
# num_bins保留的信息越多越好
将选取最佳分箱个数的过程包装为函数
def graphforbestbin(DF, X, Y, n=5,q=20,graph=True):
    '''
    自动最优分箱函数,基于卡方检验的分箱

    参数:
    DF: 需要输入的数据
    X: 需要分箱的列名
    Y: 分箱数据对应的标签 Y 列名
    n: 保留分箱个数
    q: 初始分箱的个数
    graph: 是否要画出IV图像

    区间为前开后闭 (]

    '''
    
    DF = DF[[X,Y]].copy()

    DF["qcut"],bins = pd.qcut(DF[X], retbins=True, q=q,duplicates="drop")
    coount_y0 = DF.loc[DF[Y]==0].groupby(by="qcut").count()[Y]
    coount_y1 = DF.loc[DF[Y]==1].groupby(by="qcut").count()[Y]
    num_bins = [*zip(bins,bins[1:],coount_y0,coount_y1)]

    for i in range(q):
        if 0 in num_bins[0][2:]:
            num_bins[0:2] = [(
                num_bins[0][0],
                num_bins[1][1],
                num_bins[0][2]+num_bins[1][2],
                num_bins[0][3]+num_bins[1][3])]
            continue

        for i in range(len(num_bins)):
            if 0 in num_bins[i][2:]:
                num_bins[i-1:i+1] = [(
                    num_bins[i-1][0],
                    num_bins[i][1],
                    num_bins[i-1][2]+num_bins[i][2],
                    num_bins[i-1][3]+num_bins[i][3])]
                break
        else:
            break

    def get_woe(num_bins):
        columns = ["min","max","count_0","count_1"]
        df = pd.DataFrame(num_bins,columns=columns)
        df["total"] = df.count_0 + df.count_1
        df["percentage"] = df.total / df.total.sum()
        df["bad_rate"] = df.count_1 / df.total
        df["good%"] = df.count_0/df.count_0.sum()
        df["bad%"] = df.count_1/df.count_1.sum()
        df["woe"] = np.log(df["good%"] / df["bad%"])
        return df

    def get_iv(df):
        rate = df["good%"] - df["bad%"]
        iv = np.sum(rate * df.woe)
        return iv

    IV = []
    axisx = []
    while len(num_bins) > n:
        pvs = []
        for i in range(len(num_bins)-1):
            x1 = num_bins[i][2:]
            x2 = num_bins[i+1][2:]
            pv = scipy.stats.chi2_contingency([x1,x2])[1]
            pvs.append(pv)

        i = pvs.index(max(pvs))
        num_bins[i:i+2] = [(
            num_bins[i][0],
            num_bins[i+1][1],
            num_bins[i][2]+num_bins[i+1][2],
            num_bins[i][3]+num_bins[i+1][3])]

        bins_df = pd.DataFrame(get_woe(num_bins))
        axisx.append(len(num_bins))
        IV.append(get_iv(bins_df))
        
    if graph:
        plt.figure()
        plt.plot(axisx,IV)
        plt.xticks(axisx)
        plt.xlabel("number of box")
        plt.ylabel("IV")
        plt.show()
    return bins_df
对所有特征进行分箱选择
model_data.columns
"""
Index(['SeriousDlqin2yrs', 'RevolvingUtilizationOfUnsecuredLines', 'age',
       'NumberOfTime30-59DaysPastDueNotWorse', 'DebtRatio', 'MonthlyIncome',
       'NumberOfOpenCreditLinesAndLoans', 'NumberOfTimes90DaysLate',
       'NumberRealEstateLoansOrLines', 'NumberOfTime60-89DaysPastDueNotWorse',
       'NumberOfDependents', 'qcut'],
      dtype='object')
"""

for i in model_data.columns[1:-1]:
    print(i)
    graphforbestbin(model_data,i,"SeriousDlqin2yrs",n=2,q=20)

image.png

计算各箱的WOE并映射到数据中

auto_col_bins = {"RevolvingUtilizationOfUnsecuredLines":6,
                "age":5,
                "DebtRatio":4,
                "MonthlyIncome":3,
                "NumberOfOpenCreditLinesAndLoans":5}
 
#不能使用自动分箱的变量
hand_bins = {"NumberOfTime30-59DaysPastDueNotWorse":[0,1,2,13]
            ,"NumberOfTimes90DaysLate":[0,1,2,17]
            ,"NumberRealEstateLoansOrLines":[0,1,2,4,54]
            ,"NumberOfTime60-89DaysPastDueNotWorse":[0,1,2,8]
            ,"NumberOfDependents":[0,1,2,3]}
 
#保证区间覆盖使用 np.inf替换最大值,用-np.inf替换最小值 
#原因:比如一些新的值出现,例如家庭人数为30,以前没出现过,改成范围为极大值之后,这些新值就都能分到箱里边了
hand_bins = {k:[-np.inf,*v[:-1],np.inf] for k,v in hand_bins.items()}
bins_of_col = {}
 
# 生成自动分箱的分箱区间和分箱后的 IV 值
 
for col in auto_col_bins:
    bins_df = graphforbestbin(model_data,col
                             ,"SeriousDlqin2yrs"
                             ,n=auto_col_bins[col]
                             #使用字典的性质来取出每个特征所对应的箱的数量
                             ,q=20
                             ,graph=False)
    bins_list = sorted(set(bins_df["min"]).union(bins_df["max"]))
    #保证区间覆盖使用 np.inf 替换最大值 -np.inf 替换最小值
    bins_list[0],bins_list[-1] = -np.inf,np.inf
    bins_of_col[col] = bins_list
    
#合并手动分箱数据    
bins_of_col.update(hand_bins)
 
bins_of_col
"""
{'RevolvingUtilizationOfUnsecuredLines': [-inf,
  0.09905982831820331,
  0.2978753766811289,
  0.4653785849632207,
  0.9820244866500002,
  0.9999998999999999,
  inf],
 'age': [-inf, 36.0, 54.0, 61.0, 74.0, inf],
 'DebtRatio': [-inf,
  0.0175061132,
  0.4015951649678128,
  1.4688993235185674,
  inf],
 'MonthlyIncome': [-inf, 0.1, 6167.0, inf],
 'NumberOfOpenCreditLinesAndLoans': [-inf, 1.0, 3.0, 5.0, 17.0, inf],
 'NumberOfTime30-59DaysPastDueNotWorse': [-inf, 0, 1, 2, inf],
 'NumberOfTimes90DaysLate': [-inf, 0, 1, 2, inf],
 'NumberRealEstateLoansOrLines': [-inf, 0, 1, 2, 4, inf],
 'NumberOfTime60-89DaysPastDueNotWorse': [-inf, 0, 1, 2, inf],
 'NumberOfDependents': [-inf, 0, 1, 2, inf]}
"""

data = model_data.copy()
 
#函数pd.cut,可以根据已知的分箱间隔把数据分箱
#参数为 pd.cut(数据,以列表表示的分箱间隔)
data = data[["age","SeriousDlqin2yrs"]].copy()
 
data["cut"] = pd.cut(data["age"],[-np.inf, 36.0, 54.0, 61.0, 74.0, np.inf])
data

image.png

#将数据按分箱结果聚合,并取出其中的标签值
data.groupby("cut")["SeriousDlqin2yrs"].value_counts()

"""
cut           SeriousDlqin2yrs
(-inf, 36.0]  1                   24795
              0                   14797
(36.0, 54.0]  1                   51506
              0                   39070
(54.0, 61.0]  0                   15743
              1                   12175
(61.0, 74.0]  0                   20344
              1                    7435
(74.0, inf]   0                    7737
              1                    1406
Name: SeriousDlqin2yrs, dtype: int64
"""

# 使用unstack()来将树状结构变成表状结构
data.groupby("cut")["SeriousDlqin2yrs"].value_counts().unstack()

image.png

bins_df = data.groupby("cut")["SeriousDlqin2yrs"].value_counts().unstack()

bins_df["woe"] = np.log((bins_df[0]/bins_df[0].sum())/(bins_df[1]/bins_df[1].sum()))

bins_df

image.png

def get_woe(df,col,y,bins):
    df = df[[col,y]].copy()
    df["cut"] = pd.cut(df[col],bins)
    bins_df = df.groupby("cut")[y].value_counts().unstack()
    woe = bins_df["woe"] = np.log((bins_df[0]/bins_df[0].sum())/(bins_df[1]/bins_df[1].sum()))
    return woe

#将所有特征的WOE存储到字典当中
woeall = {}
for col in bins_of_col:
    woeall[col] = get_woe(model_data,col,"SeriousDlqin2yrs",bins_of_col[col])
    
woeall

"""
{'RevolvingUtilizationOfUnsecuredLines': cut
 (-inf, 0.0991]     2.202126
 (0.0991, 0.298]    0.668970
 (0.298, 0.465]    -0.125724
 (0.465, 0.982]    -1.075822
 (0.982, 1.0]      -0.475662
 (1.0, inf]        -2.032914
 dtype: float64,
 'age': cut
 (-inf, 36.0]   -0.520053
 (36.0, 54.0]   -0.280179
 (54.0, 61.0]    0.253175
 (61.0, 74.0]    1.002752
 (74.0, inf]     1.701429
 dtype: float64,
 'DebtRatio': cut
 (-inf, 0.0175]     1.523794
 (0.0175, 0.402]    0.035165
 (0.402, 1.469]    -0.388177
 (1.469, inf]       0.175920
 dtype: float64,
 'MonthlyIncome': cut
 (-inf, 0.1]      1.330940
 (0.1, 6167.0]   -0.218048
 (6167.0, inf]    0.270218
 dtype: float64,
 'NumberOfOpenCreditLinesAndLoans': cut
 (-inf, 1.0]   -0.836758
 (1.0, 3.0]    -0.332920
 (3.0, 5.0]    -0.057622
 (5.0, 17.0]    0.123964
 (17.0, inf]    0.464595
 dtype: float64,
 'NumberOfTime30-59DaysPastDueNotWorse': cut
 (-inf, 0.0]    0.352558
 (0.0, 1.0]    -0.871069
 (1.0, 2.0]    -1.379556
 (2.0, inf]    -1.544336
 dtype: float64,
 'NumberOfTimes90DaysLate': cut
 (-inf, 0.0]    0.234569
 (0.0, 1.0]    -1.752176
 (1.0, 2.0]    -2.258120
 (2.0, inf]    -2.401458
 dtype: float64,
 'NumberRealEstateLoansOrLines': cut
 (-inf, 0.0]   -0.393291
 (0.0, 1.0]     0.194089
 (1.0, 2.0]     0.618576
 (2.0, 4.0]     0.386065
 (4.0, inf]    -0.292362
 dtype: float64,
 'NumberOfTime60-89DaysPastDueNotWorse': cut
 (-inf, 0.0]    0.124547
 (0.0, 1.0]    -1.387000
 (1.0, 2.0]    -1.764356
 (2.0, inf]    -1.811287
 dtype: float64,
 'NumberOfDependents': cut
 (-inf, 0.0]    0.627636
 (0.0, 1.0]    -0.581528
 (1.0, 2.0]    -0.528492
 (2.0, inf]    -0.477878
 dtype: float64}
"""
#不希望覆盖掉原本的数据,创建一个新的DataFrame,索引和原始数据model_data一模一样
model_woe = pd.DataFrame(index=model_data.index)

#将原数据分箱后,按箱的结果把WOE结构用map函数映射到数据中
model_woe["age"] = pd.cut(model_data["age"],bins_of_col["age"]).map(woeall["age"])

#对所有特征操作可以写成:
for col in bins_of_col:
    model_woe[col] = pd.cut(model_data[col],bins_of_col[col]).map(woeall[col])
    
#将标签补充到数据中
model_woe["SeriousDlqin2yrs"] = model_data["SeriousDlqin2yrs"]
model_woe.head()

image.png

vali_woe = pd.DataFrame(index=vali_data.index)
 
for col in bins_of_col:
    vali_woe[col] = pd.cut(vali_data[col],bins_of_col[col]).map(woeall[col])
vali_woe["SeriousDlqin2yrs"] = vali_data["SeriousDlqin2yrs"]
 
vali_X = vali_woe.iloc[:,:-1]
vali_y = vali_woe.iloc[:,-1]

建模与模型验证

X = model_woe.iloc[:,:-1]
y = model_woe.iloc[:,-1]
 
from sklearn.linear_model import LogisticRegression as LR
 
lr = LR().fit(X,y)
lr.score(vali_X,vali_y)

c_1 = np.linspace(0.01,1,20)
c_2 = np.linspace(0.01,0.2,20)
 
score = []
for i in c_1: 
    lr = LR(solver='liblinear',C=i).fit(X,y)
    score.append(lr.score(vali_X,vali_y))
plt.figure()
plt.plot(c_1,score)
plt.show()

image.png

lr.n_iter_  #array([5], dtype=int32)

score = []
for i in [1,2,3,4,5,6]: 
    lr = LR(solver='liblinear',C=0.025,max_iter=i).fit(X,y)
    score.append(lr.score(vali_X,vali_y))
plt.figure()
plt.plot([1,2,3,4,5,6],score)
plt.show()

image.png

import scikitplot as skplt
 
#%%cmd
#pip install scikit-plot
 
vali_proba_df = pd.DataFrame(lr.predict_proba(vali_X))
skplt.metrics.plot_roc(vali_y, vali_proba_df,
                        plot_micro=False,figsize=(6,6),
                        plot_macro=False)

image.png

制作评分卡

建模完毕,我们使用准确率和ROC曲线验证了模型的预测能力。接下来就是要讲逻辑回归转换为标准评分卡了。评分卡中的分数,由以下公式计算

image.png

其中A与B是常数,A叫做“补偿”,B叫做“刻度”, 代表了一个人违约的可能性。其实逻辑回归的结果取对数几率形式会得到,即我们的参数*特征矩阵,所以其实就是我们的参数。两个常数可以通过两个假设的分值带入公式求出,这两个假设分别是:

  • 某个特定的违约概率下的预期分值

  • 指定的违约概率翻倍的分数(PDO)

    • 例如,假设对数几率为时设定的特定分数为600,PDO=20,那么对数几率为时的分数就是620。带入以上线性表达式,可以得到:

      image.png

B = 20/np.log(2)
A = 600 + B*np.log(1/60)
 
B,A
# (28.85390081777927, 481.8621880878296)

有了A和B,分数就很容易得到了。其中不受评分卡中各特征影响的基础分,就是将截距作为带入公式进 行计算,而其他各个特征各个分档的分数,也是将系数带入进行计算:

base_score = A - B*lr.intercept_    # lr.intercept_:截距
base_score                          # array([481.56390143])

score_age = woeall["age"] * (-B*lr.coef_[0][1])#lr.coef_:每一个特征建模之后得出的系数
score_age#"age"特征中每个箱对应的分数

"""
cut
(-inf, 36.0]   -11.263674
(36.0, 54.0]    -6.068313
(54.0, 61.0]     5.483448
(61.0, 74.0]    21.718289
(74.0, inf]     36.850735
dtype: float64
"""
  • 保存评分卡文件
file = "./ScoreData.csv"
 
#open是用来打开文件的python命令,第一个参数是文件的路径+文件名,如果你的文件是放在根目录下,则你只需要文件名就好
#第二个参数是打开文件后的用途,"w"表示用于写入,通常使用的是"r",表示打开来阅读
#首先写入基准分数
#之后使用循环,每次生成一组score_age类似的分档和分数,不断写入文件之中
 
with open(file,"w") as fdata:
    fdata.write("base_score,{}\n".format(base_score))
for i,col in enumerate(X.columns):#[*enumerate(X.columns)]
    score = woeall[col] * (-B*lr.coef_[0][i])
    score.name = "Score"
    score.index.name = col
    score.to_csv(file,header=True,mode="a")

image.png

可以发现,真正建模的部分不多,更多是我们如何处理数据,如何利用统计和机器学习的方法将数据调整成我们希望的样子,所以除了算法,更加重要的是我们能够达成数据目的的工程能力.