ccc-sklearn-5-逻辑回归

1.逻辑回归概述

逻辑回归的本质是一个返回对数几率的，在线性数据上表现优异的分类器，它主要被应用在金融领域。即求解能够让模型对数据拟合程度最高的参数θ值，以此构建预测函数 y(x)，然后将特征矩阵输入预测函数来计算出逻辑回归的结果y。逻辑回归也可以做多分类。

逻辑回归的联系函数Sigmoid
用于将线性回归方程z变换为g(z)，并且g(z)的值分布与（0,1）之间，方程式如下：
$$g(z)=\frac{1}{1+e^{-z}}$$

说明：

• S型函数，能将任何实数映射到（0,1）区间
• 可以当做归一化的一种方法，与MinMaxSclaer的区别在于不能缩放到0和1只能趋近
• 为什么该模型实际被称为“对数几率回归”，下图可以看到y（x）的形似几率取对数本质为线性回归z
  
  几率说明：
  几率odds的本质是$\frac{p}{1-p}$，p为事件发生的概率

其他说明：

• 逻辑回归对线性关系的拟合效果非常好，计算速度很快。但是如果数据是非线性的，那么千万不要用逻辑回归
• 逻辑回归分类结果不是固定的0和1，而是小数形式呈现的类似概率数字

一直以来，人们都是以返回概率的方式来理解逻辑回归。可以说逻辑回归返回的数字，即便本质上不是概率，却有着概率的各种性质，可以被当成是概率来看待和使用。

2.sklearn中的逻辑回归

逻辑回归相关的类	说明
linear_model.LogisticRegression	逻辑回归分类器（又叫logit回归，最大熵分类器）
linear_model.LogisticRegressionCV	带交叉验证的逻辑回归分类器
linear_model.logistic_regression_path	计算Logistic回归模型以获得正则化参数的列表
linear_model.SGDClassifier	利用梯度下降求解的线性分类器（SVM，逻辑回归等等）
linear_model.SGDRegressor	利用梯度下降最小化正则化后的损失函数的线性回归模型
metrics.log_loss	对数损失，又称逻辑损失或交叉熵损失
其他会涉及的类	说明
metrics.confusion_matrix	混淆矩阵，模型评估指标之一
metrics.roc_auc_score	ROC曲线，模型评估指标之一
metrics.accuracy_score	精确性，模型评估指标之一

二元逻辑回归的损失函数
使用损失函数来衡量参数θ的模型拟合训练集时产生的信息损失的大小，从而衡量参数θ的优劣。损失函数越小，说明模型在训练集上表现越优异，拟合越充分，参数就优秀。我们也追求损失函数最小的参数θ。其极大似然估计的推导式为：
$$J(\theta )=-\sum _{i=1}^m(y_i\ast log(y_{\theta }(x_i))+(1-y_i)\ast log(1-y_{\theta }(x_i)))$$
其中θ表示求解出的一组参数，m是样本个数，$y_i$是样本i的真实标签，$y_{\theta }(x_i))$是样本i基于参数θ计算的逻辑回归返回值，$x_i$是样本i各个特征的取值。

说明：

• 由于追求损失函数的最小，可能导致过拟合；对于过拟合的控制可以通过正则化来实现
• 对于没有求解参数需求的模型没有损失函数，如KNN，决策树等

二元逻辑回归损失函数的数学解释，公式推导
二元逻辑回归标签服从伯努利分布（0-1分布），将一个特征向量为x，参数为θ的模型中一个样本的预测情况如下表示：


假设样本i的真实标签$y_i$为1，如果P1为1，就代表样本i的标签预测与真实值一致。此时对于单样本i来说，模型的预测就是完全准确；同样，如果P1此时为0 ，那么模型预测就完全错误。对于两种取值的概率整合可以定义如下等式：
$$P(\hat{y}|x_i,\theta )=P_1^{y_i}\ast P_0^{1-y_i}$$
这个式子同时表示了P1和P0，当真实标签yi为1,1-yi为0，P0的0次方就是1，结果就等于P1，如果P1确实为1则模型的拟合效果很好。同理，当yi为0时，结果为P0，如果P0非常接近1，则模型的效果很好，损失很小。所以，我们总是希望$P(\hat{y}|x_i,\theta )$的值等于1，即它的最大值。从而将模型拟合最小化损失问题转换成对函数求解极值的问题。对于一个训练集的m个样本，所有样本在特征矩阵X和参数θ组成的预测函数中，预测所有可能的$\hat{y}$的概率P为：

对于概率P取对数：

结果为交叉熵函数，添加符号将极大值问题转为极小值问题，让θ称为函数自变量，即得到损失函数J(θ):

说明：

• 推导过程即为极大似然法
• 对于$P(\widehat{y_i}|x_i,\theta )$,如果θ已知，xi未知，便称P是在探索不同特征取值下获取所有可能$\widehat{y_i}$的可能性，即概率，研究自变量与因变量之间的关系。如果特征向量xi已知，参数θ未知，便称P是在探索不同参数下获取所有可能$\widehat{y_i}$的可能性，这种可能性被称为似然，研究参数取值与因变量之间的关系。

正则化参数penalty &C
常用L1正则化和L2正则化，通过在损失函数后加上参数向量θ的L1范式和L2范式的倍数来实现。增加的范式被称为“正则项”/“惩罚项”。损失函数改变，基于损失函数的最优化来求解的参数取值必然改变，以此调节模型的拟合程度，L1、L2范式表示如下：

J(θ)为损失函数，C为用来控制正则化程度的参数，n是方程中特征的总数，j代表每个参数，j>=1

也可以写作：

参数	说明
penalty	可以输入"l1"或"l2"来指定使用哪一种正则化方式，不填写默认"l2"。 注意，若选择"l1"正则化，参数solver仅能够使用求解方式”liblinear"和"saga“ 若使用“l2”正则化，参数solver中所有的求解方式都可以使用。
C	C正则化强度的倒数，必须是一个大于0的浮点数，不填写默认1.0，默认正则项与损失函数的比值是1：1。 C越小，损失函数会越小，模型对损失函数的惩罚越重，正则化的效力越强，参数θ会逐渐被压缩得越来越小

L1正则化会将参数压缩为0，L2正则只会让参数尽量小， 不会取到0。即L1正则化本质是特征选择，掌管参数的“稀疏性”，通过减少特征个数防止过拟合。L2正则化会尽量让每个特征都对模型有贡献，但携带信息少且对模型贡献不大的参数会非常接近0。通常要优先考虑L2正则化

L1与L2正则化曲线实例
步骤一：导入库和数据

from sklearn.linear_model import LogisticRegression as LR
from sklearn.datasets import load_breast_cancer
import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score

data = load_breast_cancer()
X = data.data
y = data.target
X.shape

步骤二：实例化模型，并查看L1、L2正则化后特征对应参数

lrl1 = LR(penalty="l1",solver="liblinear",C=0.5,max_iter=1000)
lrl2 = LR(penalty="l2",solver="liblinear",C=0.5,max_iter=1000)

lrl1 = lrl1.fit(X,y)
lrl1.coef_
lrl2 = lrl2.fit(X,y)
lrl2.coef_

步骤三：比较两者在训练集合测试集上的准确率

l1 = []
l2 = []
l1test = []
l2test = []

Xtrain, Xtest, Ytrain, Ytest = train_test_split(X,y,test_size=0.3,random_state=420)

for i in np.linspace(0.05,1,19):
    lrl1 = LR(penalty="l1",solver="liblinear",C=i,max_iter=1000)
    lrl2 = LR(penalty="l2",solver="liblinear",C=i,max_iter=1000)
    
    lrl1 = lrl1.fit(Xtrain, Ytrain)
    l1.append(accuracy_score(lrl1.predict(Xtrain),Ytrain))
    l1test.append(accuracy_score(lrl1.predict(Xtest),Ytest))
    lrl2 = lrl2.fit(Xtrain, Ytrain)
    l2.append(accuracy_score(lrl2.predict(Xtrain),Ytrain))
    l2test.append(accuracy_score(lrl2.predict(Xtest),Ytest))
    
graph = [l1,l2,l1test,l2test]
color = ['green','black','lightgreen','gray']
label = ['L1','L2','L1test','L2test']

plt.figure(figsize=(6,6))
for i in range(len(graph)):
    plt.plot(np.linspace(0.05,1,19),graph[i],color[i],label=label[i])
plt.legend(loc=4)
plt.show()

在这个数据集中，两则正则化结果区别不大。随着C的增大，正则化强度越来越小，模型在训练集和测试集的表现都变好，但C大于0.8后开始下降，也就是出现的过拟合。

3.逻辑回归中的特征工程

由于我们需要通过逻辑回归的结果来判断什么样的特征与分类结果相关，因此我们希望保留特征的原貌。PCA和SVD的降维结果是不可解释的，所以一般逻辑回归汇总不使用PCA和SVD

embedded嵌入法SelectFromModel
步骤一：导入库和数据

from sklearn.linear_model import LogisticRegression as LR
from sklearn.datasets import load_breast_cancer
import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import cross_val_score
from sklearn.feature_selection import SelectFromModel

data = load_breast_cancer()
data.data.shape

步骤二：简单对比嵌入法提取的特征和全部特征的效果

LR_ = LR(solver="liblinear",C=0.9,random_state=420)
cross_val_score(LR_,data.data,data.target,cv=10).mean()

X_embedded = SelectFromModel(LR_,norm_order=1).fit_transform(data.data,data.target)
X_embedded.shape
cross_val_score(LR_,X_embedded,data.target,cv=10).mean()

说明：

• SelectFromModel类中的参数threshold是嵌入法的阈值，表示删除所有参数的绝对值低于这个阈值的特征。默认None，可以绘制threshold的学习曲线，观察不同threshold下模型的变化。可以使用模型属性.coef_（衡量特征的重要程度和贡献度）中生成的各个特征的系数来进行选择。

步骤三：绘制threshold的学习曲线

fullx = []
fsx = []
threshold = np.linspace(0,abs((LR_.fit(data.data,data.target).coef_)).max(),20)
k=0
for i in threshold:
    X_embedded = SelectFromModel(LR_,threshold=i).fit_transform(data.data,data.target)
    fullx.append(cross_val_score(LR_,data.data,data.target,cv=5).mean())
    fsx.append(cross_val_score(LR_,X_embedded,data.target,cv=5).mean())
    print((threshold[k],X_embedded.shape[1]))
    k+=1
    
plt.figure(figsize=(20,5))
plt.plot(threshold,fullx,label="full")
plt.plot(threshold,fsx,label="feature selection")
plt.xticks(threshold)
plt.legend()
plt.show()

随着threshold越来越大，被删除的特征越来越多，模型的效果也越来越差，说明调整该参数比较无效。

步骤四：绘制C的学习曲线

fullx = []
fsx = []
C = np.arange(0.01,10.01,0.5)

for i in C:
    LR_ = LR(solver="liblinear",C=i,random_state=420)
    fullx.append(cross_val_score(LR_,data.data,data.target,cv=10).mean())
    X_embedded = SelectFromModel(LR_,norm_order=1).fit_transform(data.data,data.target)
    fsx.append(cross_val_score(LR_,X_embedded,data.target,cv=10).mean())

print(max(fsx),C[fsx.index(max(fsx))])
    
plt.figure(figsize=(20,5))
plt.plot(C,fullx,label="full")
plt.plot(C,fsx,label="feature selection")
plt.xticks(C)
plt.legend()
plt.show()

继续细化最大值6.01左右的图像

fullx = []
fsx = []
C = np.arange(6.05,7.05,0.005)

for i in C:
    LR_ = LR(solver="liblinear",C=i,random_state=420)
    fullx.append(cross_val_score(LR_,data.data,data.target,cv=10).mean())
    X_embedded = SelectFromModel(LR_,norm_order=1).fit_transform(data.data,data.target)
    fsx.append(cross_val_score(LR_,X_embedded,data.target,cv=10).mean())

print(max(fsx),C[fsx.index(max(fsx))])
    
plt.figure(figsize=(20,5))
plt.plot(C,fullx,label="full")
plt.plot(C,fsx,label="feature selection")
plt.xticks(C)
plt.legend()
plt.show()

使用返回的效果最好的C值代入模型

LR_ = LR(solver="liblinear",C=6.079999999999999,random_state=420)
cross_val_score(LR_,data.data,data.target,cv=10).mean()

LR_ = LR(solver="liblinear",C=6.079999999999999,random_state=420)
X_embedded = SelectFromModel(LR_,norm_order=1).fit_transform(data.data,data.target)
cross_val_score(LR_,X_embedded,data.target,cv=10).mean()

X_embedded.shape

可以看到，在特征减少的情况下，准确率得到了提高。效果明显

其他的特征处理方法

• 系数累加法：找出曲线由锐利变平滑的转折点，转折点之前被累加的特征都是需要的，转折点之后的都不需要。相对麻烦，要提前对特征系数进行排序并确保排序后各系数对应的初始位置
• 包装法：直接设定需要的特征个数，具体参考数据预处理和特征工程

4.梯度下降

逻辑回归目的是求解能让损失函数$J(\theta )$最小化的θ值。对于二元逻辑回归，可以采用梯度下降、坐标下降、牛顿法求解。

梯度下降求解逻辑回归
梯度是一个向量，既有大小又有方向。大小就是偏导数组成的向量的大小；方向是损失函数值变化最快的方向

想像小球从高处滚下到最低点的过程，小球坐标$({\theta }_1,{\theta }_2,J)$梯度向量的方向就是每次滚动的方向；由于每次位置都会发生变化，即每次滚动方向也不同。所以梯度下降，其实就是在众多$[{\theta }_1,{\theta }_2]$可能的值中遍历，一次次求解坐标点的梯度向量，不断让损失函数的取值逐渐逼近最小值，再返回这个最小值对应的参数取值$[{\theta }_1^{\ast },{\theta }_2^{\ast }]$的过程。

说明：

• 在多元函数上对各个自变量求偏导，并以向量的形式写出来，就是梯度。对于损失函数$J({\theta }_1,{\theta }_2)$来说，自变量是逻辑回归预测函数$y_{\theta }(x)$的参数θ1，θ2，在损失函数上对θ1，θ2求偏导，得到梯度向量$gradJ({\theta }_1，{\theta }_2)$
• 逻辑回归损失函数表示如下：
  $$J(\theta )=-\sum _{i=1}^m(y_i\ast log(y_{\theta }(x_i))+(1-y_i)\ast log(1-y_{\theta }(x_i)))$$
  对函数上自变量θ求偏导，可以得到梯度向量在第j组坐标点上的表示形式：
  $$\frac{\partial }{\partial {\theta }_j}J(\theta )=d_j=\sum _{i=i}^m(y_{\theta }(x_i)-y_i)x_{ij}$$

只要给定一组θ的取值再代入特征矩阵，就能求出一组θ下的预测结果，结合真实标签向量y，可以获得θj取值下的梯度向量，大小表示为dj。遍历θ的过程表示为：
$${\theta }_{j+1}={\theta }_j-\alpha \ast d_j={\theta }_j-\alpha \ast \sum _{i=i}^m(y_{\theta }(x_i)-y_i)x_{ij}$$
$\alpha$为步长，控制每走一步后θ的变化，从而影响迭代后的梯度大小和方向。

步长概念
步长不是物理距离，也不是梯度下降过程中距离的直接变化，它是梯度向量的大小d上的一个比例，影响参数θ每次迭代后改变的部分。

从A运动到B，参数向量θ的变化为${\theta }_a-{\theta }_b$，根据梯度向量迭代公式为步长*梯度向量的大小，为二维平面三角形中的邻边。损失函数减少的量，$J({\theta }_b)-J({\theta }_a)$是二维平面三角形中的对边。
步长可以调节损失函数下降的速率，在损失函数降低的方向上，步长越长，θ的变动就越大。

逻辑回归中用参数max_iter代替步长，它表示最大迭代次数

max_iter实例

l2 = []
l2test = []

Xtrain, Xtest, Ytrain, Ytest = train_test_split(X,y,test_size=0.3,random_state=420)
for i in np.arange(1,201,10):
    lrl2 = LR(penalty="l2",solver="liblinear",C=0.9,max_iter=i)
    lrl2 = lrl2.fit(Xtrain,Ytrain)
    l2.append(accuracy_score(lrl2.predict(Xtrain),Ytrain))
    l2test.append(accuracy_score(lrl2.predict(Xtest),Ytest))

graph = [l2,l2test]
color = ["black","gray"]
label = ["L2","L2test"]
plt.figure(figsize=(20,5))
for i in range(len(graph)):
    plt.plot(np.arange(1,201,10),graph[i],color[i],label=label[i])
plt.legend(loc=4)
plt.xticks(np.arange(1,201,10))
plt.show()

warinning是指函数没有收敛，可以增大max_iter。实际情况中以预测效果为基准。

二元回归与多元回归：solver & multi_class

首先理解一对多和多对多
OvR：某种分类看作1，其他都为0
MvM：某些分类看作1，其他都为0

multi_class
auto：表示根据分类情况和其他参数来确定分类问题的类型。默认
multinomial：表示处理多分类问题，solver是“liblinear”时不可用

solver

简单对比“multinomial”和“over”

from sklearn.datasets import load_iris
iris = load_iris()

for multi_class in ('multinomial','auto'):
    clf = LR(solver='sag',max_iter=100,random_state=42,
             multi_class=multi_class).fit(iris.data,iris.target)
    print("training score: %.3f (%s)"%(clf.score(iris.data,iris.target),multi_class))

5.银行评分卡实例

项目介绍
评分卡在银行的借贷中经常使用，通过对客户的信用大小评分来区分风险用户与和正常用户，对于评分卡的完整模型的开发，需要有以下流程：

本次项目核心在数据清洗与模型开发！

导入库与数据

%matplotlib inline
import numpy as np
import pandas as pd
from sklearn.linear_model import LogisticRegression as LR

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

步骤一.数据的初步探索

data.shape

data.info()

11个特征的解释说明：

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

步骤二.数据预处理

2.1去重并恢复索引

data.drop_duplicates(inplace=True)
data.info()

data.index = range(data.shape[0])
data.info()

第一张图为没有恢复索引的效果
2.2平均值填补NumberOfDependents
由于NumberOfDependents缺的很少，所以直接用平均值来进行填充处理

data.isnull().sum()/data.shape[0]
data["NumberOfDependents"].fillna(int(data["NumberOfDependents"].mean()),inplace=True)
data.info()

2.3逻辑森林填充MonthlyIncome
在实际情况中如此大量的数据缺失肯定需要我们与相关业务进行沟通解决。这里采用其他特征值预测的思想来填补这个缺失的特征。

#X 要填充的特征矩阵 y 完整的无缺失值的标签  to_fill 要填补的那一列缺失值
def fill_missing_rf(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"]
X.shape

y_pred =fill_missing_rf(X,y,"MonthlyIncome")
data.loc[data.loc[:,"MonthlyIncome"].isnull(),"MonthlyIncome"]=y_pred
data.info()

2.4处理描述性错误的异常值

上面数据可以看出，年龄最小值出现了问题。还有NumberOfTimes90DaysLate这类特征中也出现了不符合常理的数字，这些都需要进行处理

(data["age"] == 0).sum()
data = data[data["age"]!=0]

data[data.loc[:,"NumberOfTime60-89DaysPastDueNotWorse"]> 90].count()
data  = data[data.loc[:,"NumberOfTime60-89DaysPastDueNotWorse"] < 90]
data.index=range(data.shape[0])
data.info()

由于异常特征值数量都很少，所以这里就直接删除了

2.5处理样本不均衡的问题

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))

违约的人毕竟只是极少数人，样本统计肯定会出现不均衡的情况，这里采用上采样（这里不进行解释）的方法来平衡样本

import imblearn
from imblearn.over_sampling import  SMOTE

sm = SMOTE(random_state=42)
X,y = sm.fit_resample(X,y)
n_sample_ =X.shape[0]
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))

2.6 划分训练集和数据集并保存方便后续使用

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"./data/model_data.csv")
vali_data.to_csv(r"./data/vali_data.csv")

到这里，数据预处理就做完了

步骤三：分箱操作

分箱是制作评分卡最核心的步骤，本质就是离散化连续变量，让拥有不同属性的人被分为不同的类别，类似于聚类

如何判断合适的分箱个数

离散化连续变量必然伴随信息的损失，IV(information value)定义特征信息量对预测函数贡献的衡量值
$$IV=\sum _{i=1}^N(good\%-bad\%)\ast {WOE}_i$$
N表示特征分箱个数，good%代表标签为0的优质客户占整个特征中所有优质客户的比例，bad%是这个箱子中的坏客户占整个特征中所有坏客户的比例，而$WO{E}_i$(weight of Evidence),写作：
$$WO{E}_i=ln(\frac{good\%}{bad\%})$$

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

箱子越多，信息损失越多，IV必然越小；所以分箱时应该计算每个特征在每个箱子数目下的WOE值，利用IV值的曲线，找出合适的分箱个数

说明：

• 期望的分箱效果：组间差异大，组内差异小。使用卡方检验两个箱子的相似性，并合并差异性很小的箱子
• 特征分箱的基本过程：
  1.首先将连续变量分成一组数量较多的分类型变量
  2.确保每一组中都包含有2种类别的样本，这是为了IV值的计算
  3.对相邻组卡方检验，对检验结果P值很大的组合并，知道数据中组的数量小于设定的N箱为止
  4.让一个特征分别分成[2,3,4…20]箱，观察每个分享个数下的IV值如何变化，找出最合适的分享个数
  5.分享完毕后，计算每个箱的WOE值，观察分箱效果
  最后对各个特征进行分箱，观察每个特征的IV值，以此挑选特征

3.1 等频分箱（针对age）

model_data["qcut"], updown = pd.qcut(model_data["age"], retbins=True, q=20)
#retbins：True，同时返回索引和箱子的Series
#model_data:每个样本所分的箱子  updown：所有箱子的上下限

#统计每个分箱中0,1数量
cnt_y0 = model_data[model_data["SeriousDlqin2yrs"] == 0].groupby(by="qcut").count()["SeriousDlqin2yrs"]
cnt_y1 = model_data[model_data["SeriousDlqin2yrs"] == 1].groupby(by="qcut").count()["SeriousDlqin2yrs"]
cnt_y0
#分别为每个区间上界，下界，0出现次数，1出现次数
num_bins = [*zip(updown,updown[1:],cnt_y0,cnt_y1)]
num_bins
#zip按照最短列进行结合

3.2确保IV值计算合并分箱

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
    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

这个案例不需要，基本思想就是如果第一次不行，就往后合并；之后的向前合并即可
3.3定义WOE和IV函数

def get_woe(nums_bins):
    columns = ["min","max","count_0","count_1"]
    df = pd.DataFrame(nums_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

注意区分bad_rate和bad%，分别是一个箱中坏样本占的比例以及一个箱中坏样本占整个特征中坏样本的比例

3.4卡方检验，合并箱体，绘制IV曲线

num_bins_ = num_bins.copy()
import matplotlib.pyplot as plt
import scipy

IV = []
axisx = []

while len(num_bins_) > 2:
    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 = 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()

一般选择下降趋势变化最大的，这张图可以选择4箱作为分箱结果
3.5 箱子合并函数

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]
            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,4)

3.6将整个选择最佳分箱个数过程封装为函数

def graphforbestbin(DF, X, Y, n=5, q=20, graph=True):
    DF = DF[[X,Y]].copy()
    DF["qcut"],bins = pd.qcut(DF[X], retbins=True, q=q, duplicates="drop")
    cnt_y0 = DF.loc[DF[Y]==0].groupby(by="qcut").count()[Y]
    cnt_y1 = DF.loc[DF[Y]==1].groupby(by="qcut").count()[Y]
    num_bins = [*zip(bins,bins[1:],cnt_y0,cnt_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(nums_bins):
        columns = ["min","max","count_0","count_1"]
        df = pd.DataFrame(nums_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_df1 = pd.DataFrame(get_woe(num_bins))
        axisx.append(len(num_bins))
        IV.append(get_iv(bins_df1))
        
    if graph:
        plt.figure()
        plt.plot(axisx,IV)
        plt.xticks(axisx)
        plt.xlabel("number of box")
        plt.ylabel("IV")
        plt.show()
        
    return -1

3.7对所有特征都绘制IV曲线

model_data.columns

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

选取特征绘制IV曲线正常：

选取特征绘制IV曲线不正常：

3.8对所有特征进行分箱的选择
分别可以自动分箱和观察后手动分箱

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]} 
hand_bins = {k:[-np.inf,*v[:-1],np.inf] for k,v in hand_bins.items()}

bins_of_col ={}
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"]))
    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

步骤四：计算各箱的WOE并映射到数据中

4.1计算并存储每箱的woe值

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

woeall = {}
for col in bins_of_col:
    woeall[col] = get_woe(model_data,col,"SeriousDlqin2yrs",bins_of_col[col])

4.2将WOE映射到原始数据中

model_woe = pd.DataFrame(index=model_data.index)

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"]

步骤五：建模与模型验证

5.1处理测试集数据并建模

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)

5.2绘制参数C的学习曲线

c_1 = np.linspace(0.01,1,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_2,score)
plt.show()

5.3 查看ROC曲线

import scikitplot as skplt

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)

步骤六：制作评分卡

评分卡公式如下：A、B是常数，分别为补偿和刻度，log（odds）表示一个人违约的可能性，也是参数
$$Score=A-B\ast log(odds)$$
A与B可以通过两个假设的分值代入求得：

1. 某个特定违约概率下的预期分值
2. 指定的违约概率翻倍的分数（PDO）

将所有特征的评分卡内容全部写入本地文件
6.1通过方程组（自己假设的）解AB

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

6.2计算基础分base_score

base_score = A - B*lr.intercept_
base_score

6.3计算全部特征的评分卡并保存

file = "scoredata.csv"
with open(file,"w") as fdata:
    fdata.write("base_score,{}\n".format(base_score))
for i,col in 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")

终于做完了，呜呜呜呜呜~~~~~~~~~~~~~

其他附录表

逻辑回归参数表：

逻辑回归属性表：

逻辑回归接口列表