友链

传统业务算法流程 V2.0

该博客主要面向算法工程师，意在经验总结和提升算法工作效率，非广义科普机器学习知识。

一、前期准备
  • 1. 业务评估
    • 明确甲方对算法模型的需求
      • 输入输出
      • 上线时间
    • 明确甲方对算法模型效果的期望
      • 准确率
      • precision
      • recall
      • 可解释性
      • 接口响应速度
    • 明确开发资源
      • 确认数据资源准备情况
      • 确认集群资源（CPU、MEM、GPU）准备情况
    • 按照数据评估和算法调研的结果，评估
      • 需求合理性
      • 所需数据量
      • 所需集群资源
      • 所需工时
      • 计划的结果返回形式
  • 2. 数据评估
    • 跟甲方确认各个字段的业务含义以及每个字段特征值的业务含义，询问是否有类别特征值码表。
    • 重点检查各个标签的数据量是否达标（一般按照单个标签5000条训练样本的标准）
    • 重点检查各个特征的缺失率情况
      • 是否符合业务经验
      • 确认缺失值出现的业务原因（为后续缺失值填充做准备）
      • 确认缺失率高的字段是否存在数据加工问题
    • 重点检查各个特征相在每个标签下的分布情况
      • 防止信息泄露
      • 针对分布显著不一致的特征，跟甲方进行确认
        • 是否符合业务经验
        • 是否存在数据加工问题
        • 是否存在信息泄露
    • 重点检查各个特征之间的相关性程度
      • 输出 correlation 矩阵
      • 对于高相关性的特征，跟甲方沟通确认
        • 是否符合业务经验
        • 是否可以删除部分特征
二、数据探索性分析（EDA）
  • 1. 为什么要做EDA
    • 使模型能够业务可解释
    • 初步判断特征重要性
    • 初步判断特征之间的关联关系，为构建新特征做准备
    • 初步进行异常值评估
      • 异常值一般为离群值
      • 大值和最小值不一定为异常值
    • EDA 的主要手段
      • 可视化画图
      • 量化统计
  • 2. 缺失值分析

    需要跟也无妨确认原因，可能是漏采，也可能因为数值为0

    连续型特征

    • 0填充：缺失值本身自带业务含义
    • 组内中位数、平均数填充：根据业务规则对样本进行恰当地分组，取组内统计量填充
    • 最邻近填充：一般在时间序列模型中使用最邻近值填充
    • 中位数、平均数填充

    非连续型特征

    • None填充：缺失值本身自带业务含义
    • 组内高频值填充
    • 高频值填充

    其它

    • 字段解析
      • 通过对其它字段的解析，映射到当前字段。
      • 举例：“Mr.Chen” → 性别：男。其中姓名字段或许没有作为特征的价值，但是可以从中提取有价值的信息。
  • 3. 异常值分析

    连续型特征

    • 异常值、离群值
      • 异常值在某些业务中应被理解为“离群值”，因为最大值和最小值不等同于离群值。
      • image.png
      • 如图所示
        • 红色标记样本远远偏离均值，可以被轻易判定为异常点。
        • 绿色标记样本并非极大值或者极小值，且比某些正常值还要小，但是因为偏离了当前“群体”，所以也需要被判定为异常点。
    • BoxPlot
      • 基本概念
        • 下四分位数(Q1) = 7
        • 中位数(Q2) = 8.5
        • 上四分位数(Q3) = 9
        • 四分位间距(ΔQ) = Q3 - Q1
        • 最小值 = Q1-1.5ΔQ
        • 最大值 = Q3+1.5ΔQ
        • 之外的就是 outliers
        • image.png
      • Python 画图
        ax = sns.boxplot("Age", data=data)
        image.png
        ax = data.boxplot("Age")
        image.png

    非连续型特征

    • 码表映射，是否存在该特征值。
    • 统计频率，极低频率特征值可视为异常值删除。

    异常值处理

    • 删除含有异常值的样本
      • 特别是针对数据采集导致的异常样本
    • 根据 boxplot 的最大值最小值进行数值截断
  • 4. 特征关联性分析

    统计工具（常用）

    • 交叉表（crosstab）
      • 交叉表是一种特殊的透视表，往往用来统计频次（分组），也可以使用参数【aggfunc】指定聚合函数实现其他功能。
      pd.crosstab(data['Pclass'],data['Survived'])

      image.png

    • 相关系数
      • Pearson相关系数（连续变量一对一）
        • Pearson相关系数优点：排除了特征间量纲的差异。
        • 输出范围为-1到+1, 0代表无相关性，负值为负相关，正值为正相关。
          image.png
        corr = data.corr()

        image.png

      • VIF（Variance Inflation Factor，方差膨胀因子）（连续变量多对一）
        • VIF：每个因子，用其他Ｎ个因子进行回归解释。
          • 将一个特征作为标签 y
          • 将剩余特征作为训练数据 X
          • 然后进行拟合，【拟合优度】高，说明 y 和其他特征 X 之间存在高度相关性
        • 拟合优度：测量的是线性回归中得到的模型对样本的拟合程度，使用 R^2 来表示
          image.png
          • SSR：预测值和均值的残差和平方
          • SSE：真实值和预测值的残差和平方
          • SST：真实值和均值的残差和平方
        • VIF = 1 / (1 - R^2)
        import pandas as pd
import numpy as np
from statsmodels.stats.outliers_influence import variance_inflation_factor
 
# 宽表
data = pd.DataFrame([[15.9,16.4,19,19.1,18.8,20.4,22.7,26.5,28.1,27.6,26.3]
                     ,[149.3,161.2,171.5,175.5,180.8,190.7,202.1,212.1,226.1,231.9,239]
                     ,[4.2,4.1,3.1,3.1,1.1,2.2,2.1,5.6,5,5.1,0.7]
                     ,[108.1,114.8,123.2,126.9,132.1,137.7,146,154.1,162.3,164.3,167.6]]).T
  
# 上述宽表，第0列为因变量，1、2、3列为自变量
X = data[[1,2,3]]
  
# ✨✨✨务必注意✨✨✨，一定要加上常数项，如果没有常数项列，计算结果天差地别，可能VIF等于好几千
X[4]=1
print(X)
  
#计算第2个变量的(第二列)的方差膨胀因子，注意python的计数是从0开始的，所以第二个参数应该是1
variance_inflation_factor(X[[1,2,3,4]].values, 1)
        • 评判标准
          • VIF <= 1: 不存在多重共线性
          • 1 <= VIF <= 5: 存在一定程度共线性
          • VIF > 5: 存在高度共线性
        • 可以使用递归的方式，逐步筛选高相关性特征。
      • correlation_ratio（相关比）（连续型-非连续型）
        • 相关比：用来计算连续型和非连续型变量间的相关系数大小
          image.png
        • nx = 特征值为x的总样本量
          image.png
        • yxi = 非连续型特征值为x对应的连续型特征值
        from dython import nominal
nominal.correlation_ratio(data["Sex"], data["Age"])
# 0.12041709277881193
      • 卡方检验（非连续型变量）

        • 它主要是比较两个分类变量的关联性分析。其根本思想就是在于比较理论频数和实际频数的吻合程度或拟合优度问题。

        • 卡方检验的 p-value 会受样本量的影响，样本量越大卡方越大，p-value 越小，所以不建议使用。

        • 举例
          • 检验【职业类别】和【咖啡口味】的相关性
          • 假设显著性阈值为 α=0.01
            image.png
          • 假设 H0：职业与咖啡无关系(独立)
            image.png
            • e11 = 330/500 * 56 = 36.96
            • e12 = 125/500 * 21 = 14
              image.png
          • 在使用 chi2_contingency 前，需要先用 pd.crosstab 生成 contingency table
          from  scipy.stats import chi2_contingency
import numpy as np
import pandas as pd
data=[[25,21,10],[82,88,30],[223,16,5]]
 
df=pd.DataFrame(data,index=['美式咖啡','拿铁咖啡','卡布奇诺'],columns=['IT','行政','工程'])
kt=chi2_contingency(df)
 
print('卡方值=%.4f, p值=%.4f, 自由度=%i expected_frep=%s'%kt)

          image.png

          • p值小于0.01，则拒绝原假设，认为咖啡和职业不独立性。
      • Cramer's V（克莱姆相关系数）（非连续型变量）
        • 在卡方检验的基础上，克莱姆相关系数加入了对数据量的考量，规避了卡方检验的缺点，推荐使用。
          image.png
          image.png
        from dython import nominal
nominal.cramers_v(data["Sex"], data["Cabin"], bias_correction=True)
# 0.19804091740334265

    可视化工具（常用）

    • 特征间关联关系（A: 连续，B: 非连续）

      • A
        • df.plot(kind="kde")
        ax = data[["Age"]].plot(kind="kde")

        image.png

        • sns.distplot()
        plt.figure(figsize=[8,5])
ax = sns.distplot(data['Age'].dropna().values, bins=range(0, 81, 1), kde=False, color='blue')

        image.png

      • B
        • df.plot(kind="bar")
        ax = data['Pclass'].value_counts().plot(kind='bar')

        image.png

        • sns.barplot()
        vc = data['Pclass'].value_counts()
ax = sns.barplot(vc.index, vc.values, ci=None)

        image.png

        • dataframe.plot(kind="pie")
        ax = data['Pclass'].value_counts().plot(kind='pie', explode=[0,0.1,0.3], autopct='%1.1f%%', shadow=True)

        image.png

      • AA
        • df.plot(kind="scatter",x,y)
        ax = data[["Age","Fare"]].plot(kind="scatter",x="Age",y="Fare")

        image.png

        • sns.scatterplot()
        ax = sns.scatterplot(data=data, x=data["Age"], y=data["Fare"])

        image.png

      • BB
        • sns.countplot()
        ax = sns.countplot(data["Pclass"], edgecolor='black', hue="Survived", linewidth=1, data=data)

        image.png

        • sns.histplot()
        ax = sns.histplot(data=data,x='Initial', hue='Survived', multiple='stack')

        image.png

      • AB
        • sns.violinplot()
        ax = sns.violinplot(x="Pclass", y="Age", data=train, split=True)

        image.png

        • sns.boxplot()
        ax = sns.boxplot(x="Pclass", y="Age", data=data)

        image.png

      • AAB
        • sns.scatterplot()
        ax = sns.scatterplot(data=data, x=data["Age"], y=data["Fare"], hue="Survived", palette=['red','blue'])

        image.png

      • ABB
        • sns.violinplot
        ax = sns.violinplot(x="Pclass", y="Age", hue="Survived", data=train, split=True)

        image.png

        • sns.boxplot()
        sns.boxplot(x="Pclass", y="Age", hue="Survived", data=data)

        image.png

        • sns.factorplot
        ax = sns.factorplot(x="Pclass", y="Age", hue="Survived", data=data, aspect=1.5, size=3.5, ci=95.0)

        image.png

      • BBB
        • sns.factorplot()
        ax = sns.factorplot(x="Pclass", y="Sex", hue="Survived", data=data, aspect=1.5, size=3.5, ci=95.0)

        image.png

      • ABBB
        • sns.factorplot()
        ax = sns.factorplot(x="Pclass", y="Age", hue="Sex", col="Embarked", data=data, aspect=0.9, size=3.5, ci=95.0)

        image.png

      • BBBB
        • sns.factorplot()
        sns.factorplot(x="Pclass", y="Survived", hue="Sex", col="Embarked", data=train, aspect=0.9, size=3.5, ci=95.0)

        image.png

    • Dython
      • dython.nominal.associations
        • 可以统计非连续型特征和连续型特征间的关联关系
      from dython import nominal
nominal.associations(data,figsize=(20,10),mark_columns=True);

      image.png

    • HeatMap 
      plt.figure(figsize=[7,7])
ax = sns.heatmap(data.corr(), vmax=0.6, square=True, annot=True)

      image.png

    • PairPlot 
      sns.pairplot(data=train, vars=['Survived','Pclass','Age','SibSp','Parch','Fare'], hue='Survived', palette=['red','blue'], size=1.5)

      image.png

    • Secondary_y Plot 
      # 默认顺序为 【x,y1,y2】
df.columns = ['ip', 'click_count', 'prop_downloaded']
# 必须要先排序
df = df.sort_values('is_attributed', ascending=False)
# 画图
ax = df.plot(secondary_y='prop_downloaded')
plt.show()

      image.png

  • 5. 新特征挖掘
    • 分组归纳
      • 基于业务，但是需要Python进行统计分析，所以由算法完成
      • 检查非连续型特征的特征值唯一性，若有大量重复特征值则可进行分组分析。
      • 基于业务，观察分组后的目标特征是否具有可归纳的组内特性。
      • 举例：经典Titanic项目中的【共用船票】特征
        • 通过对【船票】特征的唯一性分析，发现有多个重复的船票。
        • 根据业务，确定多个需要分析的目标特征，其中包括性别、年龄、【消费金额】等等。
        • 通过对【船票】特征的特征值进行 groupBy 分组统计后发现，使用相同船票的乘客都有着几乎相同的【消费金额】。
        • 故将【共用船票】作为新特征加入到模型当中。
    • 信息抽取
      • 基于业务，从已有特征当中提取新的重要信息，尤其是字符串类别的特征。
      • 基于业务，检查特征值中是否存在值得提取的信息。
      • 举例：经典Titanic项目中的【姓名】特征
        • 【姓名】特征中包含了一些【头衔】信息，例如：Mr、Mrs、Doctor 等等。
        • 我们可以通过解析这些【头衔】，并结合业务，将其转换为【性别】和【年龄】。
    • 特征平铺
      • 基于算法，部分重要且特征值数量较小的非连续型特征可以进行平铺。
      • 把类别特征的每个特征值都以 0/1 OneHot 的形式展开成为新特征。
      • 举例：经典Titanic项目中的【客舱等级】特征
        • 【客舱等级】特征对于旅客是否能够幸存十分重要，且一共只有3个特征值：A，B，C。
        • 将【客舱等级】特征进行 OneHot 平铺，生成新特征：Pclass_A, Pclass_B 和 Pclass_C
        dummy_variables = pd.get_dummies(data["Pclass"])
三、特征工程
  • 1. 数据采样
    • 数据不均衡

      一般情况下，正负样本比例小于 1:100 则可以认为有显著的标签分布不均衡。

      产生的问题

      • 无法评估特征关联性（以 Pearson系数为例）
        • 正负样本比 1:1000
          image.png
        • 正负样本比 1:1
          image.png
      • 模型训练一边倒
        • 在训练标签极度不均衡的情况下，以正负样本比 1:1000为例，正样本所产生的的 loss 占总 loss 比例过小，以至于被模型完全忽略。
        • classification_report 结果中，recall趋向于0，precision趋向于1，auc趋向于 0.5 。
    • 解决方案
      • 上采样

        SMOTE

        from imblearn.over_sampling import SMOTE

label_dict = {0:517900, 1:327100} #指定每个标签上采样后的个数
smote = SMOTE(sampling_strategy=label_dict, random_state=2022)
X_sm, y_sm = smote.fit_resample(df_features, df_label)

        image.png

      • 下采样

        随机下采样

        • 建议做多次随机采样后分别训练，预测时取多个模型的平均概率。
        from sklearn.utils import resample
data_undersampled = resample(data,
                          n_samples=n_samples, # 下采样后的期望样本数量
                          random_state=2022)

        NearMiss

        • sampling_strategy
          • version=1: 选取正例样本中与N个最近邻负样本平均距离最短的样本
          • version=2: 选取正例样本中与N个最远邻负样本平均距离最短的样本
          • version=3: 2-steps 算法,首先对于每个负类样本,他们的M近邻正样本保留。然后，选取到N近邻负平均聚类最大
          • dict: 指定每个类别采样后的数量
          from imblearn.under_sampling import NearMiss

X = df[feature_names]
y = df[label_name]
resample_dict = {"label1":10, "label2": 100, ...}
 
nm1 = NearMiss(version=1, random_state=2022, sampling_strategy=resample_dict)
X_resampled_num1, y_resampled = nm1.fit_resample(X, y)
  • 2. 特征降维

    特征降维一般可以解决多重共线性和维度灾难问题

    • 导致特征系数不再准确，产生偏差，进而导致模型可解释性变差。

    • 导致模型预测结果方差变大，更不稳定。

    • 过高纬度的特征可能会导致内存溢出、训练速度慢和预测接口响应慢等问题

    • 特征选择
      • filter（过滤法）
        • 进行单特征分析，例如：分布分析、频率统计等等
        • 进行多特征间关联性分析，例如：VIF、Cramers'v、相关比系数等等
        • 选择最佳特征子集
      • wrapper（包装法）
        • 选择不同特征子集训练模型
        • 对比 feature_importance 和模型指标
        • 选择最佳的特征子集
      • embedded（嵌入法）
        • 将特征选择的过程嵌入到模型本身当中
        • 例如：L1 norm、L2 norm
    • 降维算法

      主要对比最常用到的 PCA 和 t-SNE

      • 均为无监督算法
      • t-SNE 的结果相比 PCA 更具有代表性。
      • 因为 t-SNE 的速度极慢且极消耗资源，在处理大量数据（100000+）时，推荐先使用 PCA 进行降维后再使用 t-SNE。
      • t-SNE 以 KL散度函数作为 loss函数，而 KL散度函数不是凸函数，所以不同的 random_seed 会让算法收敛在不同的局部最优点，因此推荐用不同的seed多尝试几次。
        • KL散度：一般用于度量两个概率分布函数之间的“距离”
          image.png
      • image.png
      • PCA
        from sklearn.decomposition import PCA

pca = PCA(n_components=2, svd_solver='full', random_state=2022)
data_pca = pca.fit_transform(data[feature_con])
      • t-SNE
        from sklearn.manifold import TSNE

tsne = TSNE(n_components=2, init='pca', random_state=2022, perplexity=30, method='barnes_hut', n_iter=1000, verbose=1)
data_tsne = tsne.fit_transform(data[feature_con])
        • n_component：降维后的期望维度
        • perplexity：困惑度。在进行流行学习（manifold learning）时，用来计算的最邻近 N 个样本点，数据总量越多，perplexity 应越大。
        • init：default=pca; t-SNE未保留明确的全局结构，所以通过用 PCA 进行初始化来缓解。
        • method：default=barnes_hut
          • Barnes-Hut：时间复杂度为 O(NlogN)，n_component 只能为 2或者3。
          • exact：时间复杂度为 O(N^2)，不适用于 100w+ 数据量。
        • exact：时间复杂度为 O(N^2)，不适用于 100w+ 数据量。
          n_jobs：决定并行运行的线程数量。
      • 二维结果绘制
        import matplotlib.pyplot as plt

plt.figure(1, figsize=(10, 10))
colors = ['blue', 'red']
for color, i, target_name in zip(colors, [0, 1], [0, 1]):
    plt.scatter(data_pca[y == i, 0], data_pca[y == i, 1], color=color, s=1,
                alpha=.8, label=target_name, marker='.')
  • 3. 类别特征数值化
    • OneHotEncoding
      • 基于算法，部分重要且特征值数量较小的非连续型特征可以进行平铺。
      • 把类别特征的每个特征值都以 0/1 OneHot 的形式展开成为新特征。
      • 举例：经典Titanic项目中的【客舱等级】特征
        • 【客舱等级】特征对于旅客是否能够幸存十分重要，且一共只有3个特征值：A，B，C。
        • 将【客舱等级】特征进行 OneHot 平铺，生成新特征：Pclass_A, Pclass_B 和 Pclass_C
          dummy_variables = pd.get_dummies(data["Pclass"])
    • LabelEncoding
      • LabelEncoding 将类别特征转换为了数值特征，缺点是在建模过程当中会带入数值的大小信息。
        df["feature_name"] = pd.factorize(df["feature_name"])[0]
  • 4. 特征标准化
    • 业务经验总结
      • 所有 TreeBased 模型不需要使用特征标准化，例如：Xgboost、LightGBM、RandomForest等等。
      • 虽然 RobustScaler 适用于异常值较多的情况，但是也要尝试其它两种标准化算法，在实际应用当中，往往 MinMaxScaler 效果最佳
      • 常规流程会存在一定程度的数据泄露，参考：Pipeline+CV
    • MinMaxScaler
      • 适用于所有场景。
      • 根据在项目上的多次使用，MinMaxScaler 的效果一直较好。
        from sklearn import preprocessing
import numpy as np
 
x = np.array([[   3000.,   -1000.,    2000.,  613000.],
       [   2000.,       0.,       0.,   50000.],
       [      0.,    1000.,   -1000., -113000.],
       [   1000.,    2000.,   -3000.,    6000.]])
 
 
minmax_scaler = preprocessing.MinMaxScaler()
x_mm = minmax_scaler.fit_transform(x)
print(x_mm)

        [[1.         0.         1.         1.        ]
 [0.66666667 0.33333333 0.6        0.22451791]
 [0.         0.66666667 0.4        0.        ]
 [0.33333333 1.         0.         0.16391185]]
    • StandardScaler
      • 适用于异常值较少的场景。
      • 如果数据中，有异常数据，如绝对值特别大的数据，如果把他们作为正常数据，会导致其他数据放缩后，非常非常的小，接近于0。
        from sklearn import preprocessing
import numpy as np
 
x = np.array([[   3000.,   -1000.,    2000.,  613000.],
       [   2000.,       0.,       0.,   50000.],
       [      0.,    1000.,   -1000., -113000.],
       [   1000.,    2000.,   -3000.,    6000.]])
 
 
standard_scaler = preprocessing.StandardScaler()
x_ss = standard_scaler.fit_transform(x)
print(x_ss)

        [[ 1.34164079 -1.34164079  1.38675049  1.69234455]
 [ 0.4472136  -0.4472136   0.2773501  -0.3177609 ]
 [-1.34164079  0.4472136  -0.2773501  -0.89972748]
 [-0.4472136   1.34164079 -1.38675049 -0.47485617]]
    • RobustScaler
      • 对异常值鲁棒，适用于异常值较大或者较多的特征缩放场景。
      • 使用中位数和四分位区间进行缩放。Xi = (xi - median) / (75th quantile - 25th quantile)
        from sklearn import preprocessing
import numpy as np
 
# 定义数据集
x = np.array([[   3000.,   -1000.,    2000.,  613000.],
       [   2000.,       0.,       0.,   50000.],
       [      0.,    1000.,   -1000., -113000.],
       [   1000.,    2000.,   -3000.,    6000.]])
 
# 定义放缩对象
robust_scaler = preprocessing.RobustScaler(with_centering=True, with_scaling=True, quantile_range=(25.0, 75.0), copy=True)
 
# 开始放缩
x_robust = robust_scaler.fit_transform(x)
 
# 打印放缩后的结果
print(x_robust)

        [[ 1.         -1.          1.25        2.72727273]
 [ 0.33333333 -0.33333333  0.25        0.1025641 ]
 [-1.          0.33333333 -0.25       -0.65734266]
 [-0.33333333  1.         -1.25       -0.1025641 ]]
  • 5. 数据压缩

    • 几乎所有 Kaggle top solutions 的必备环节

    • 能够有50%~80%的内存释放

    • 加速模型训练

      data['Age'] = data['Age'].astype(np.int8)
    • 案例 
      def reduce_mem_usage2(df):
    """ iterate through all the columns of a dataframe and modify the data type
        to reduce memory usage.       
    """
    start_mem = df.memory_usage().sum() / 1024**2
    print('Memory usage of dataframe is {:.2f} MB'.format(start_mem))
      
    for col in df.columns:
        col_type = df[col].dtype
          
        if col_type != object:
            c_min = df[col].min()
            c_max = df[col].max()
            if str(col_type)[:3] == 'int':
                if c_min > np.iinfo(np.int8).min and c_max < np.iinfo(np.int8).max:
                    df[col] = df[col].astype(np.int8)
                elif c_min > np.iinfo(np.int16).min and c_max < np.iinfo(np.int16).max:
                    df[col] = df[col].astype(np.int16)
                elif c_min > np.iinfo(np.int32).min and c_max < np.iinfo(np.int32).max:
                    df[col] = df[col].astype(np.int32)
                elif c_min > np.iinfo(np.int64).min and c_max < np.iinfo(np.int64).max:
                    df[col] = df[col].astype(np.int64) 
            else:
                if c_min > np.finfo(np.float16).min and c_max < np.finfo(np.float16).max:
                    df[col] = df[col].astype(np.float16)
                elif c_min > np.finfo(np.float32).min and c_max < np.finfo(np.float32).max:
                    df[col] = df[col].astype(np.float32)
                else:
                    df[col] = df[col].astype(np.float64)
        else:
            df[col] = df[col].astype('category')
  
    end_mem = df.memory_usage().sum() / 1024**2
    print('Memory usage after optimization is: {:.2f} MB'.format(end_mem))
    print('Decreased by {:.1f}%'.format(100 * (start_mem - end_mem) / start_mem))
      
    return df
 
reduce_mem_usage2(df)

      Memory usage of dataframe is 1970.87 MB
Memory usage after optimization is: 547.14 MB
Decreased by 72.2%
Memory usage of dataframe is 1673.87 MB
Memory usage after optimization is: 460.02 MB
Decreased by 72.5%
CPU times: user 1min 53s, sys: 5min 48s, total: 7min 41s
Wall time: 7min 42s
四、模型构建
  • 1. 模型选择

    预留，该部分为经验积累和案例展示

  • 2. 模型训练
    • 数据集切分
      • 建议将总样本 6:2:2 分割成 Train Set、Valid Set 和 Test Set，并固定 random_state。
      • 只有 Train Set 可以进行采样等预处理，Valid Set 和 Train Set 不做采样处理！！！
      from sklearn.model_selection import StratifiedKFold

sss = StratifiedKFold(n_splits=5, random_state=None, shuffle=False)
for train_index, test_index in sss.split(X, y):
    original_Xtrain, original_Xtest = X.iloc[train_index], X.iloc[test_index]
    original_ytrain, original_ytest = y.iloc[train_index], y.iloc[test_index]
    • Xgboost
      • 导入关键依赖包（xgboost==1.5.2）
        import xgboost as xgb
from xgboost import rabit
from xgboost.callback import TrainingCallback
from typing import Callable, List, Optional, Union, Dict, Tuple, TypeVar, cast, Sequence, Any
      • 定义 callbacks
        • 支持根据 data_valid 的指标进行 early_stop
        • 支持将模型输出同时打印在屏幕上和输入到 logger 中
                    def get_logger(log_file_path, level=logging.DEBUG, when='midnight', back_count=0):
    """
    :brief  日志记录
    :param log_filename: 日志名称
    :param level: 日志等级
    :param when: 间隔时间:
        S:秒
        M:分
        H:小时
        D:天
        W:每星期（interval==0时代表星期一）
        midnight: 每天凌晨
    :param back_count: 备份文件的个数，若超过该值，就会自动删除
    :return: logger
    """
    logger = logging.getLogger(log_file_path)
    logger.setLevel(level)
    # log输出格式
    formatter = logging.Formatter('%(asctime)s - %(levelname)s: %(message)s')
    # 输出到控制台
    ch = logging.StreamHandler()
    ch.setLevel(level)
    # 输出到文件
    fh = logging.FileHandler(
        filename=log_file_path,
        encoding='utf-8')
    fh.setLevel(level)
    # 设置日志输出格式
    fh.setFormatter(formatter)
    ch.setFormatter(formatter)
    # 添加到logger对象里
    logger.addHandler(fh)
    logger.addHandler(ch)
    return logger
 
_Model = Any
class EvaluationMonitor(TrainingCallback):
    '''Print the evaluation result at each iteration.
    .. versionadded:: 1.3.0
    Parameters
    ----------
    metric :
        Extra user defined metric.
    rank :
        Which worker should be used for printing the result.
    period :
        How many epoches between printing.
    show_stdv :
        Used in cv to show standard deviation.  Users should not specify it.
    '''
    def __init__(self, rank: int = 0, period: int = 1, show_stdv: bool = False) -> None:
        self.printer_rank = rank
        self.show_stdv = show_stdv
        self.period = period
        assert period > 0
        # last error message, useful when early stopping and period are used together.
        self._latest: Optional[str] = None
        super().__init__()
 
    def _fmt_metric(
        self, data: str, metric: str, score: float, std: Optional[float]
    ) -> str:
        if std is not None and self.show_stdv:
            msg = f"\t{data + '-' + metric}:{score:.5f}+{std:.5f}"
        else:
            msg = f"\t{data + '-' + metric}:{score:.5f}"
        return msg
 
    def after_iteration(self, model: _Model, epoch: int,
                        evals_log: TrainingCallback.EvalsLog) -> bool:
        if not evals_log:
            return False
 
        msg: str = f'[{epoch}]'
        if rabit.get_rank() == self.printer_rank:
            for data, metric in evals_log.items():
                for metric_name, log in metric.items():
                    stdv: Optional[float] = None
                    if isinstance(log[-1], tuple):
                        score = log[-1][0]
                        stdv = log[-1][1]
                    else:
                        score = log[-1]
                    msg += self._fmt_metric(data, metric_name, score, stdv)
            msg += '\n'
 
            if (epoch % self.period) == 0 or self.period == 1:
                rabit.tracker_print(msg)
                logger.info(msg)
                self._latest = None
            else:
                # There is skipped message
                self._latest = msg
        return False
 
    def after_training(self, model: _Model) -> _Model:
        if rabit.get_rank() == self.printer_rank and self._latest is not None:
            rabit.tracker_print(self._latest)
            logger.info(msg)
        return model
 
log_file_path = "./model.log"
logger = get_logger(log_file_path)
 
es = xgb.callback.EarlyStopping(
    rounds=20,
    save_best=True,
    data_name="validation_0",
    metric_name="auc",
)
em = EvaluationMonitor()
      • 模型训练
        n_jobs, n_estimators, max_depth, reg_alpha, reg_lambda, scale_pos_weight=100, 500, 10, 1, 1, 1
logger.info("model parameters: \nn_jobs={}\nn_estimators={}\nmax_depth={}\nreg_alpha={}\nreg_lambda={}\nscale_pos_weight={}\n".format(n_jobs, n_estimators, max_depth, reg_alpha, reg_lambda, scale_pos_weight))
model_xgb = xgb.XGBClassifier(n_jobs=n_jobs,n_estimators=n_estimators,max_depth=max_depth,reg_alpha=reg_alpha,reg_lambda=reg_lambda,scale_pos_weight=scale_pos_weight)
model_xgb.fit(df_yly_train, label_yly_train, verbose=False, eval_metric=["auc","logloss","rmse"],
              eval_set=[(df_yly_test, label_yly_test)], callbacks=[es,em])
model_xgb.save_model(".model.json")
logger.handlers = [] # 每次结束清空 handlers
  • 3. 模型优化
    • 参数调优
      • 方法介绍

        • 对于GridSearch、RandomizedSearch和BayesianSearch，都是用于调优机器学习模型中的超参数。

        • 那么，应该用哪个？如何使用它？

        GridSearch VS RandomizedSearch VS BayesSearch

        这里我直接给出结论：（实验过程和使用细节可参考GridSearch VS RandomizedSearch VS BayesSearch）

        结论1: RandomizedSearch 永远优于 GridSearch

        对于GridSearch和RandomizedSearch，始终使用RandomizedSearch，因为RandomizedSearch性能更好（两者搜索时间相差较小）；

        原因在于RandomizedSearch可以轻松高效地在重要参数中找到优化值。

        image.png

        如上图所示，有一个重要参数和一个不重要参数。对于GridSearch的3x3 = 9个组合，实际上它在9次迭代中只搜索了3个不同的值来寻找重要参数。然而，对于RandomizedSearch，它在9次迭代中搜索9个不同的重要参数值。因此，RandomizedSearch更容易搜索重要参数。

        结论2: BayesSearch VS RandomizedSearch

        当RandomizedSearch随机搜索参数时，如果我们有意地、有方向地使用类似Gradient Descent思想去搜索会怎么样呢？

        由于BayesSearch包含了很少的条件概率计算，对于一个小型数据集，BayesSearch将比RandomizedSearch花费更多的时间。对于一个庞大的数据集（≥100000 样本），在给定复杂的搜索分布参数和足够迭代的情况下，BayesSearch比RandomizedSearch有更好的性能。

        image.png

        • 我们可以从上面的测试结果中看到
          • BayesSearchCV 的效果一直优于 RandomSearchCV
          • 随着数据量的增加，BayesSearchCV 的效率也要高于 RandomSearch
      • 样例代码

        • 加载依赖包，生成随机训练数据，定义模型

          from sklearn.datasets import make_classification
from xgboost import XGBClassifier
from sklearn.model_selection import GridSearchCV, RandomizedSearchCV, StratifiedShuffleSplit
from scipy.stats import randint, uniform
 
from skopt import BayesSearchCV
from skopt.space import Real, Categorical, Integer
 
X2, y2 = make_classification(n_samples=200000,
                             n_features=120,
                             n_redundant=50,
                             n_informative=70,
                             random_state=1)
 
model = XGBClassifier(random_state=42, n_jobs=-1)  # input hyperparams without tuning

        • RandomSearchCV

          param = {
    'learning_rate': uniform(.05, .1),
    'subsample': uniform(.2, .3),
    'n_estimators': randint(20, 70),
    'min_child_weight': randint(20, 40),
    'reg_alpha': uniform(0, .7),
    'reg_lambda': uniform(0, .7),
    'colsample_bytree': uniform(.1, .7),
    'max_depth': randint(2, 6)
}
 
randomsearch = RandomizedSearchCV(model,
                                  param_distributions=param,
                                  n_iter=30,  # specify how many iterations
                                  cv=StratifiedShuffleSplit(n_splits=1, test_size=.25, random_state=1),
                                  n_jobs=-1,
                                  scoring='accuracy',
                                  return_train_score=True)
randomsearch.fit(X2, y2)
print('best score of Randomized Search over 10 iterations:', randomsearch.best_score_)

          best score of Randomized Search over 10 iterations: 0.93264
Wall time: 9min 45s

        • BayesSearchCV

          param = {
    'learning_rate': Real(.05, .1+.05),
    'subsample': Real(.2, .5),
    'n_estimators': Integer(20, 70),
    'min_child_weight': Integer(20, 40),
    'reg_alpha': Real(0, 0+.7),
    'reg_lambda': Real(0, 0+.7),
    'colsample_bytree': Real(.1, .1+.7),
    'max_depth': Integer(2, 6)
}
 
bayessearch = BayesSearchCV(model,
                            param,
                            n_iter=30, # specify how many iterations
                            scoring='accuracy',
                            n_jobs=-1,
                            cv=StratifiedShuffleSplit(n_splits=1, test_size=.25, random_state=1)
                           )
bayessearch.fit(X2, y2)

          best score of Bayes Search over 10 iterations: 0.959
Wall time: 9min 50s
    • 阈值选择
      • Dython

        • Kolmogorov–Smirnov test

          from dython.model_utils import ks_abc

y_pred_proba = model_xgb.predict_proba(X_test)
ks_abc(y_test.values.reshape(y_test.shape[0]), y_pred_proba[:,1])

          'eopt': 0.5052752494812012

          image.png

        • ROC

          from dython.model_utils import ks_abc

y_pred_proba = model_xgb.predict_proba(X_test)
metric_graph(y_test.values.reshape(y_test.shape[0]), y_pred_proba[:,1], metric='roc', class_names=['0','1'])

          这里的 AUC 指的是 ROC_curve 所对应的 AUC

          {'1': {'auc': {'val': 0.8404166666666667, 'naive': 0.5},
  'eopt': {'val': 0.3808262, 'x': 0.17857142857142858, 'y': 0.72}},
 'ax': <AxesSubplot:title={'center':'Receiver Operating Characteristic'}, xlabel='False Positive Rate', ylabel='True Positive Rate'>}

          image.png

        • Precision-Recall

          from dython.model_utils import ks_abc

y_pred_proba = model_xgb.predict_proba(X_test)
metric_graph(y_test.values.reshape(y_test.shape[0]), y_pred_proba[:,1], metric='pr', class_names=['0','1'])

          这里的 AUC 指的是 Precision-Recall Curve 所对应的 AUC

          {'1': {'auc': {'val': 0.775124260854837, 'naive': 0.373134328358209},
  'eopt': {'val': 0.53355354, 'x': 0.7, 'y': 0.7291666666666666}},
 'ax': <AxesSubplot:title={'center':'Precision-Recall Curve'}, xlabel='Recall', ylabel='Precision'>}

          image.png

  • 4. 模型评估
    • classification_report 
      from sklearn.metrics import classification_report

classification_report(y_test, y_pred)

      image.png

    • AUC、ROC、Precision-Recall：见模型优化部分

  • 5. 模型可解释性 + EDA
    • 特征权重

      参考项目：https://www.kaggle.com/code/headsortails/pytanic

      • 不同模型建模效果对比
        image.png
      • 不同模型特征权重对比
        image.png
      • 画图展示
        image.png
        image.png
      • 总结：
        • 在模型效果基本一致的前提下，特征权重会有着很大的差异。所以在业务上，很难根据一个模型的特征权重去定义一个模型的重要性。
        • 因为特征间存在着或多或少的共线性，共线特征的重要性会被相互稀释，所以我们不能直接根据特征权重分数直接定义一个特征的重要程度。
        • 在上面的对比中：
          • RandomForest 的特征权重跟 df.corr() 中跟标签”Survived“的关联关系基本相符。
          • Adaboost 的特征权重中，部分高线性相关但特征重要性偏低的特征【”Alone“，”Parch“，”SibSp“】的权重直接变为了0。
    • SHAP
      • 全量样本
        • 图像解释
          • 每个点都是一个样本
          • 颜色越红表示数值越大
          • 在竖轴左边表示对标签“0”的贡献程度，竖轴右边表示对标签“1”的贡献程度
          • 以“Pclass”举例，“Pclass”越小的样本，被判为标签“1”的概率越高
        • SHAP
          import shap

explainer = shap.TreeExplainer(model_xgb)
feature_names = list(X_test.columns)
 
shap_values = explainer.shap_values(X_test.values)
shap.summary_plot(shap_values, X_test.values, feature_names=feature_names, alpha=0.6)
          image.png
      • 单条样本
        • SHAP 可以分析特征在单条样本预测上的重要性程度。
          import shap

explainer = shap.TreeExplainer(model_xgb)
# 单条样本的 np.array
feature_names = list(X_test.columns)
Xd = X_test.values[44].reshape([1, len(feature_names)])
shap_values = explainer.shap_values(Xd)
 
shap.summary_plot(shap_values, Xd, feature_names=feature_names)
          image.png
      • 单个特征
        • 原非连续型特征
          import shap

explainer = shap.TreeExplainer(model_xgb)
feature_names = list(X_test.columns)
shap_values = explainer.shap_values(X_test.values)
 
shap.dependence_plot('Pclass', shap_values, X_test, interaction_index="Pclass",alpha=0.4)
          image.png
      • 双特征 
        import shap

explainer = shap.TreeExplainer(model_xgb)
feature_names = list(X_test.columns)
shap_values = explainer.shap_values(X_test.values)
 
shap.dependence_plot('Pclass', shap_values, X_test, interaction_index="Sex",alpha=0.4)

        image.png

五、现存问题

数据加工（以虚开走逃为例）

  • 算法需获得原始数据，目前只有业务老师根据经验加工过后的特征数据
    • 缺陷一：业务老师做的这一步骤等同于特征工程中的新特征挖掘，但是无法确保其重要性以及全面性。
    • 缺陷二：无法保证特征加工的质量，目前已经发现了很多特征质量问题。
posted @ 2022-07-01 09:27  爱吃碳水想减肥  阅读(111)  评论(0)    收藏  举报