【数据探索分析】

提取数据后，对数据进行探索分析是非常重要的，通过探索分析，我们可以对变量的缺失、异常值、变量分布、相关性及稳定性等情况有一个整体的认识。

代码如下：

1.数据探索分析

######################### 数据探索分析
​def data_summary(data,null_rate=0.8,nunique_max_rate=0.9,nmiss_unique_num=1):
    '''
    args:
        data:数据源
        null_rate：缺失率阈值
        nunique_max_rate：单值率阈值
        nmiss_unique_num：非空值个数阈值
    return：
 返回
    '''
    import numpy as np
    import pandas as pd
    from pandas import DataFrame as df
    import time
    
    # 缺失情况
    df_info = df(data.isnull().sum(),columns = ['null_num']) # 缺失数量    
    df_info['data_length'] = len(data) # 数据的长度
    df_info['null_rate'] = df_info['null_num']/df_info['data_length'] # 缺失率
​
    # 单值率、非空值个数、数据类型等信息
    df_info['nunique_max_rate'] = data.apply(lambda x:1 if x.isnull().sum()==len(data) else max(x.value_counts())/len(data)) # 单值率
    df_info['nmiss_unique_num'] = data.nunique()-(df_info['null_rate']==0) # 非空值个数
    df_info['type'] = df(data.dtypes).applymap(str) # 数据类型
    df_info = df_info.reset_index().rename(columns = {'index':'feature_name'})
    
    # 描述性统计
    des_info = data.describe().T.reset_index().rename(columns={'index':'feature_name'})
       
    # 合并数据
    df_info = pd.merge(df_info,des_info,on='feature_name',how='left')
    
    # 按照阈值标记删除字段
    df_info['drop_col'] = np.where((df_info['null_rate']>=null_rate)|
        (df_info['nunique_max_rate']>=nunique_max_rate)|
        (df_info['nmiss_unique_num']<=nmiss_unique_num)|
        ((df_info['type']<='object')&(df_info['nmiss_unique_num']>=50)),1,np.nan)
    
    # 导出数据
    t = str(time.strftime("%Y-%m-%d %H:%M:%S", time.localtime()).replace(' ','_').replace(':','_'))
    df_info.to_excel(r'df_info_'+ t +'.xlsx',index=False)
     
    return df_info

2.变量分布分析

import seaborn as sns
​
# 法1
from matplotlib import pyplot as plt 
nrow = 2
ncol = 2
fig,axs = plt.subplots(nrows=nrow,ncols=ncol)
​
z=0
for i in range(nrow):
    print(i)
    for j in range(ncol):
        print(j)
        sns.distplot(data.iloc[:,z],ax=axs[i,j])
        z = z+1  
        
# 法2
sns.pairplot(data) # 对角线上为变量的分布，非对角线上为变量间的关系

3.变量相关性分析

- 自变量间相关性（热力图）

coeff = data.corr() # 相关系数
sns.heatmap(coeff, cmap='YlGnBu', annot=True, annot_kws={'size':5,'weight':'bold', 'color':'white'})  # cmap：颜色设置，annot：单元格内写入值，annot_kws：单元格内写入值的设置

- 自变量与目标变量的关系

sns.pairplot(data, hue='label', palette='YlGnBu') # hue：针对某一字段进行分类，palette：设置分类字段的颜色

具体作图方法可见以下链接：

【python画图_变量相关性(heatmap、pairplot)】

4.多重共线性

如果变量间的相关性较高，需要考虑去除多重共线性的问题。

SST：总平方和（真实值与均值差的平方和）
SSR：回归平方（预测值与均值差的平方和）
SSE：残差平方和（真实值与预测值差的平方和）

具体过程:
假设数据集共有a,b,c三个自变量，且变量间存在共线性，那么vif的具体过程为：

• 把自变量 a 作为因变量，b 和 c 为自变量做回归，得出回归方程的R^2，根据计算公式得到变量 a 的VIF
• 把自变量 b 作为因变量，a 和 c 为自变量做回归，得出回归方程的R^2，根据计算公式得到变量 b 的VIF
• 把自变量 c 作为因变量，a 和 b 为自变量做回归，得出回归方程的R^2，根据计算公式得到变量 c 的VIF

如果变量 c 的VIF较大，则说明变量 c 与a,b之间存在较强的共线性。

阈值范围：

• 一般样本量在10W+时,VIF>10就存在多重共线性
• 样本量在10W以下时，VIF>5就出现多重共线性
• 几千的小样本时需要VIF<2

# 计算每个变量的vif    
def cal_vif(data):
    '''
    atgs:
        data：数据源
    return:
        vif：每个变量对应的vif
    '''    
    from statsmodels.stats.outliers_influence import variance_inflation_factor
    vif = df(data.columns,columns=['feature_name'])
    vif['VIF'] =[variance_inflation_factor(data.values,i) for i in range(len(data.columns))]
    vif.sort_values('VIF',ascending=False,ignore_index=True,inplace=True)
    return vif
​
# 循环删除vif较高的变量，直至全部变量的vif小于设定阈值
vif = cal_vif(data)
while vif['VIF'].max()>2: # vif的阈值为2
    drop_col = vif['feature_name'][0] # vif按照降序排列，所以vif值最大的索引是0
    data.drop(drop_col,axis=1,inplace=True)
    vif = cal_vif(data)

5.变量/评分 稳定性

- 计算psi

# 计算psi（传入数据可为变量、分数）
def cal_psi(actual,expected,cut_method='equidistance',cut_points=[],bins=10):
    '''
    args:
        actual：实际数据（线上），Series
        expected：期望数据（线下），Series
        cut_method：分箱方法，'auto'：自动分箱，'equidistance'：等距分箱，'equalfrequency'：等频分箱
        cut_points：切点列表
        bins：分箱数
    return:
        psi_info：各区间的分布及psi
        psi：psi值
    '''
    # 1.确定分箱方法
    if cut_method=='auto':  # 自动分箱
        cut_points = [-np.inf]+cut_points+[np.inf]
    elif cut_method=='equidistance': # 等距分箱
        bins = min(bins, expected.nunique()) # 去掉空值后值的个数与给定箱数取小
        cut_points = list(np.linspace(expected.min(),expected.max(),num=bins))
        cut_points = [-np.inf]+cut_points+[np.inf]
    elif cut_method=='equalfrequency': # 等频分箱
        bins = min(bins, expected.nunique()) # 去掉空值后值的个数与给定箱数取小
        cut_points = list(np.nanquantile(expected,[i/bins for i in range(1,bins)]))
        cut_points = [-np.inf]+cut_points+[np.inf]
    
    # 2.计算各区间的样本量
    Range = ['missing']
    expected_bins_cnt = [expected.isnull().sum()]
    actual_bins_cnt = [actual.isnull().sum()]
    for i in range(len(cut_points)-1):
        print(i)
        # 区间
        left = round(cut_points[i],2)
        right = round(cut_points[i+1],2)
        Range.append('('+str(left)+','+str(right)+']')
        # 计算expected样本量
        ebins_cnt = ((expected>left)&(expected<=right)).sum()
        expected_bins_cnt.append(ebins_cnt)
        # 计算actual样本量
        abins_cnt = ((actual>left)&(actual<=right)).sum()
        actual_bins_cnt.append(abins_cnt)
    psi_info = df({'Range':Range,'Expected_cnt':expected_bins_cnt,'Actual_cnt':actual_bins_cnt})
    
    # 3.计算psi
    psi_info['Expected_pct'] = psi_info['Expected_cnt'] / len(expected)
    psi_info['Expected_pct'] = np.where(psi_info['Expected_pct']==0,0.0001,psi_info['Expected_pct'])
    
    psi_info['Actual_pct'] = psi_info['Actual_cnt'] / len(actual)
    psi_info['Actual_pct'] = np.where(psi_info['Actual_pct']==0,0.0001,psi_info['Actual_pct'])
    
    psi_info['A-E'] = psi_info['Actual_pct'] - psi_info['Expected_pct']
    psi_info['A/E'] = psi_info['Actual_pct'] / psi_info['Expected_pct']
    psi_info['ln(A/E)'] = np.log(psi_info['A/E'])
    
    psi_info['PSI'] = psi_info['A-E'] * psi_info['ln(A/E)']
    psi = psi_info['PSI'].sum()
    
    return psi_info,psi

- psi画图

# psi画图    
def plot_psi(psi_info):
    '''
    args:
        psi_info：各区间的分布及psi            
    '''
    from matplotlib import pyplot as plt
    plt.rcParams['font.sans-serif']=['SimHei'] # 显示中文
    # 画图
    plt.bar(psi_info.index-0.1,psi_info['Expected_pct'],width=0.2,color='#007d65')    
    plt.bar(psi_info.index+0.1,psi_info['Actual_pct'],width=0.2,color='#918597') 
    
    # 添加数据标签
    exp_label = psi_info['Expected_pct'].apply(lambda x:round(x,4))
    for a,b in zip(psi_info.index,exp_label):
        plt.text(a-0.2,b+0.001,b,fontsize=8)
        
    act_label = psi_info['Actual_pct'].apply(lambda x:round(x,4))
    for a,b in zip(psi_info.index,act_label):
        plt.text(a+0.001,b+0.001,b,fontsize=8)
        
    # 坐标轴、标题及图例  
    plt.xticks(psi_info.index,psi_info['Range']) # 横轴标签，索引对应标签
    plt.xlabel('Range')   # 横轴标签
    plt.ylabel('A/E_pct') # 纵轴标签
    plt.title('A/E分布图') # 图表标题
    plt.legend(labels=['Expected_pct','Actual_pct'],fontsize=8) # 图例（当图例有多个的时候需要制定图例名称）
    
    # 保存图形文件
    try:
        import os
        plt.savefig(os.getcwd() + '\\PSI.png', format='png', dpi=80)
    except:
        print('保存图片失败！')