数据竞赛修炼笔记之工业化工生产预测

这段时间，会有系列真实的竞赛项目陪伴，我会通过修炼笔记的方式记录我这段时间学习数据竞赛的经历，希望每个竞赛都能给我们带来收获和成长！ 这个故事会很长，但我会坚持往下走，你看，天上太阳正晴，如果可以，我们一起吧…

1. 写在前面

终于下定决心涉足这个纠结很久的话题了，作为一个懵懂无知的竞赛小白，其实是非常渴望参加一场数据比赛的，因为数据比赛对于AIer来说真的很重要，不知道你是否遇到过这样的一些疑惑，就是涉足一个新领域的时候，比如数据挖掘，先非常努力的花时间学习python，numpy，pandas，scipy，sklearn这些基础知识，但是当后面真的遇到实际问题的时候，却不知道如何下手去分析，再回头想掏之前学会的工具时，发现这些工具依然在，但已记不清它们的使用说明。难不成之前学习过的这些都白学了？ 其实不是的，这些工具都在，只不过我们之前没有真正的去用过它， 所以，比赛就是一个让我们大展身手的舞台，通过解决实际问题，我们才能真正掌握之前的工具，也是学习知识的一个融合，毕竟解决实际问题，需要各个领域的知识交融和碰撞。在这个过程中，我们还可以认识一些志同道合的伙伴，一起进步和交流，进行思维的碰撞并得到成长。

所以这个数据竞赛修炼系列我会把我的所思所学都记录下来并分享，一是因为解决问题的思路可以迁移和变换，这样或许会帮助更多的人，二是通过整理和总结，可以使得知识和技能在自己脑海中逗留的时间长一些吧。

首先声明这个系列的每一篇都会很长，并且会有一大波代码来袭，毕竟每一篇都是一个完整的数据竞赛，每一篇都需要很长的时间消化整理，因为我想用最朴素，最详细的语言把每个比赛的思路和代码说给你听。

今天是数据竞赛修炼笔记的第二篇，2019年天池的一个比赛, 我们这次来做一个工业化生产预测的项目，任务目标是利用异烟酸生产过程中的各参数， 预测最终异烟酸的收率， 上一次做了一个二分类，而这次是一个回归问题。数据集包括生产过程中10个步骤的参数， 样本id、A1-A28、B1-B14包括原料，辅料，时间，温度，压强及收率等， 可以看一下生产过程的流程图：

这一次我们首先会学习数据检查和问题的修正，在这里面会用到一个强大的数据探索性分析工具pandas_profiling。然后学习如何根据原来的数据做特征工程， 在这里面会有时间特征应该怎么利用，时间段特征应该如何处理， 然后使用xgboost模型对异烟酸的收率进行预测，最后会对本次比赛的知识点进行总结。

大纲如下：

• 导入数据集并进行探索性分析
• 特征工程部分（时间和时间段预处理， 特征工程，筛选符合常规数据）
• 模型的建立和训练（xgboost）
• 模型的测试与结果
• 知识点和思路的总结（这一块重点是xgboost的使用和一些处理数据的tricks）

Ok, let’s go!

2. 导入数据集并进行探索性分析

我们使用的数据集是工业上直接进行记录的数据集，但是我们使用Excel观察一下数据集之后，就会发现一些格式上的问题，那么我们相应的就必须做出处理，关于数据集的检查这块，没有什么特别好的方式，毕竟不同的任务数据不同，只能是先打开数据大体上观察一遍，找到问题所在，当然也可以编写程序，遍历去观察等。

我们下面先导入数据集，然后进行数据的检查和问题的修正工作

"""导入数据集"""
df_trn = pd.read_csv('data/jinnan_round1_train_20181227.csv', encoding='GB2312')   # 如何没有后面的编码信息会报错
df_tst_a = pd.read_csv('data/jinnan_round1_testA_20181227.csv', encoding='GB2312')
df_tst_b = pd.read_csv('data/jinnan_round1_testB_20190121.csv', encoding='GB2312')  # 测试集在这会用不到，也可以不导入

# 观察数据
df_trn.head()
df_trn.info()

"""
1396个训练样本， 44列， 最后一列是label， 有很多都存在缺失值， 会发现数据量是比较少的，
数据量少也可以做一些人为的数据增强之类的。
"""

2.1 数据的探索性分析

这一块介绍一个好的工具： pandas_profiling
这个包可以直接利用ProfileReport生成一份数据探索性报告， 在这里面会看到：

• 总体的数据信息（首先是数据集信息：变量数(列)、观察数(行)、数据缺失率、内存；数据类型的分布情况），
• 警告信息
  • 要点：类型，唯一值，缺失值
  • 分位数统计量，如最小值，Q1，中位数，Q3，最大值，范围，四分位数范围
  • 描述性统计数据，如均值，模式，标准差，总和，中位数绝对偏差，变异系数，峰度，偏度
• 单变量描述（对每一个变量进行描述）
• 相关性分析（皮尔逊系数和斯皮尔曼系数）
• 采样查看等

"""数据探索性分析"""
import pandas_profiling as ppf
ppf.ProfileReport(df_trn)

就这两句话，可是效果很强大：

这里只显示了一部分，如果数据有什么问题，这里会给出一些提示信息， 通过查看，我们初步获得的信息： 数据缺失问题，字段相关问题， 时间字段的问题，高基数警告（列中有很多唯一值）， 44列1396个样本,数据量也是比较少的，这种情况下也可以考虑人为做一些数据的增强。

2.2 数据的检查与问题的修正

拿到数据之后，首先应该对数据进行怀疑，这些数据中有没有一些有问题的量，每一列常规情况下应该是一些什么样的值，正常情况下可以画一个分布图或者箱线图等这种图，这样可以检测出一些离群点，这样就可以指定筛选的范围，去过滤不同的离群点。所以拿到数据之后，先对数据做一个检查，看看有没有什么样的问题，有问题，应该及时的做出修改。尤其是这种工业上的问题，很可能会由于记录时不仔细会出现数据的错误问题，要进行一个检查。

那怎么检查呢？ 数据量小的时候可以通过Excel简单的看一下，因为Excel表里面有个套用表格样式，套用上之后就会发现每一列的取值情况，这时候可以具体看到取得值是不是合法。 数据量比较大的时候就需要使用python写一个函数了， 遍历所有的列，然后每一列输出取值的情况，这样来检查是不是取值合法，不合法的进行处理。

所以下面我写一个python函数来检查每一列的值是不是合法。尤其注意那种时间段的列

# 这里可以遍历列，然后输出一下每一列取得是什么值，从这里面看看有什么不合适的值
for col in df_trn.columns[1:]:
    print(col, ": ", df_trn[col].value_counts().index)

找到问题之后，下面我们就先修正数据，去掉一些不合法的值。这里只给出一部分代码，这块没有什么固定的模板，毕竟不同的数据集不一样， 没有技术含量，但是需要提前做。

# 这一块没有什么固定的模板，需要自己去观察修正， 没有什么技术含量，但是必须去做的事情。

def train_abnormal_revise(data):
    df_trn = data.copy()              # 记得要复制一份处理，不要改变原来的数据
    df_trn.loc[(df_trn['A1'] == 200) & (df_trn['A3'] == 405), 'A1'] = 300   # 这样的就1个，可能是写错了，因为300对应405, 毕竟这都是配料
    # A5这一列，会发现有三个不合法的值，要替换掉
    df_trn['A5'] = df_trn['A5'].replace('1900/1/21 0:00', '21:00:00')
    df_trn['A5'] = df_trn['A5'].replace('1900/1/29 0:00', '14:00:00')
    df_trn['A9'] = df_trn['A9'].replace('1900/1/9 7:00', '23:00:00')
    # A9有两个不合法的值
    df_trn['A9'] = df_trn['A9'].replace('1900/1/9 7:00', '23:00:00')
    df_trn['A9'] = df_trn['A9'].replace('700', '7:00:00')
    # A11有一个不合法的值
    df_trn['A11'] = df_trn['A11'].replace('1900/1/1 2:30', '2:30:00')
    df_trn['A11'] = df_trn['A11'].replace(':30:00', '00:30:00')
    df_trn['A16'] = df_trn['A16'].replace('1900/1/12 0:00', '12:00:00')
    df_trn['A20'] = df_trn['A20'].replace('6:00-6:30分', '6:00-6:30')
    df_trn['A20'] = df_trn['A20'].replace('18:30-15:00', '18:30-19:00')
    # A22有个不合法的值，记错了应该是
    df_trn['A22'] = df_trn['A22'].replace(3.5, np.nan)
    df_trn['A25'] = df_trn['A25'].replace('1900/3/10 0:00', 70).astype(int)
    df_trn['A26'] = df_trn['A26'].replace('1900/3/13 0:00', '13:00:00')
    df_trn['B1'] = df_trn['B1'].replace(3.5, np.nan)
    df_trn['B4'] = df_trn['B4'].replace('15:00-1600', '15:00-16:00')
    df_trn['B4'] = df_trn['B4'].replace('18:00-17:00', '16:00-17:00')
    df_trn['B4'] = df_trn['B4'].replace('19:-20:05', '19:05-20:05')
    df_trn['B9'] = df_trn['B9'].replace('23:00-7:30', '23:00-00:30')
    df_trn['B14'] = df_trn['B14'].replace(40, 400)
    
    return df_trn

2.3 标签与数据集整合

这里的一个技巧就是训练集和测试集放在一块进行处理会比较方便

df_trn, df_tst = df_trn.copy(), df_tst_a.copy()
df_target = df_trn['收率']   # 单独把收率拿出来
del df_trn['收率']

# 下面也是类似于一个数据处理的技巧，就是训练集和测试集放在一块预处理操作会比较方便
df_trn_tst = df_trn.append(df_tst, ignore_index=False).reset_index(drop=True)

2.4 某些列的缺失值处理

对于缺失值的处理，可以使用众数，中位数，平均值等进行填充， 对于年龄的这种离散的情况，一般众数比较合理，这里A3列，也是使用众数进行的填充。

# df_trn['A3'].isnull().sum()  42

# A3列有缺失值，所以可以进行一个填充，使用众数
for _df in [df_trn, df_tst, df_trn_tst]:
    _df['A3'] = _df['A3'].fillna(405)

3. 特征工程

这一部分，要提取特征出来，虽然我们收到的数据是结构化好了的数据，但是大部分时候，这些数据不能直接用，第一会有冗余性，第二特征不足，模型可能无法学习，所以我们需要特征工程部分的处理，这一部分很重要，数据和特征决定了机器学习的上限，而模型和算法只是去逼近这个上限而已，所以模型的好坏直接取决于特征工程部分。 这一块主要可以分为下面几部分：

• 时间段和时间特征预处理
• 温度相关特征， 温度相关统计特征
• 时间相关特征
• 水耗相关特征
• 所有特征进行合并

3.1 时间和时间段特征预处理

这一块，主要是针对时间特征的一个处理，首先把所有时间相关的列拿出来，然后同时对测试集，训练集进行处理，对于每一列时间特征，我们先改名字（带上时间），这样后面处理的时候方便查找操作，对于时间段的话，可以进行拆分成起止时间和终止时间，也可以求出持续时间，换成三个特征。 这些都是tricks。对于拆分好的时间，我们转换成时长，这样计算机才认识。

# 所有时间相关列
cols_timer = ['A5', 'A7', 'A9', 'A11', 'A14', 'A16', 'A24', 'A26', 'B5', 'B7']

# 同时对训练和测试集进行相同处理, 不管用得上用不上，先处理完再说
for _df in [df_trn_tst, df_trn, df_tst]:
    # 添加列名标记
    _df.rename(columns={_col:_col+'_t' for _col in cols_timer}, inplace=True)
    # 遍历所有持续时间相关列  如21:00-21:30
    for _col in ['A20', 'A28', 'B4', 'B9', 'B10', 'B11']:
        # 获取当前列的索引
        _idx_col = _df.columns.tolist().index(_col)
        # 添加新的一列，表示起始时间， split表示分别取开始和结束时间， 用索引来指定
        _df.insert(_idx_col+1, _col+'_at', _df[_col].str.split('-').str[0])
        #添加新的一列，表示终止时间
        _df.insert(_idx_col+2, _col+'_bt', _df[_col].str.split('-').str[1])
        
        # 删除持续时间
        del _df[_col]
        cols_timer = cols_timer + [_col + '_at', _col+'_bt']

这样，我们就很容易把带时间的列筛选出来：

"""我们把带有时间特征的列筛选出来"""
cols_timer = list(filter(lambda x: x.endswith('t'), df_trn_tst.columns))  
cols_timer

下面把时间全部换成分钟的形式：

"""将时间全部全换成分钟的形式"""
def time_to_min(x):
    if x is np.nan:
        return np.nan
    else:
        x = x.replace(';', ':').replace('；', ':')
        x = x.replace('::', ':').replace('"', ':')
        h, m = x.split(':')[:2]
        h = 0 if not h else h
        m = 0 if not m else m
        return int(h)*60 + int(m)   


for _df in [df_trn_tst, df_trn, df_tst]:
    for _col in cols_timer:
        _df[_col] = _df[_col].map(time_to_min)     

3.2 正式特征工程部分

这一块就是进行特征的组合，提取，转换等操作了，主要分为温度相关特征，时间相关特征和水耗相关特征的提取等。

3.2.1 创建一个新的df来准备添加特征

之所以要创建一个新的df，是因为我们不想在原来的数据集上进一步特征工程，这样拼接太大了，毕竟我们要从不同的角度构建特征出来，所以我们建立一个新的df之后，我们就可以把新造出来的特征先加到这个空的df上，然后我们把这个和原来的进行拼接就是最后的特征了，这也是做特征工程的一种trick吧， 我们只需要统一的样本id。

raw = df_trn_tst.copy()
df = pd.DataFrame(raw['样本id'])
df.head()

3.2.2 温度相关特征

温度特征这块，也是分为四个过程，加热，水解，脱色，结晶。 在每一个过程当中，我们首先提取原始的温度特征，然后再尝试进行一些加减组合等。特征工程中，我们需要发挥想象力，尽可能多的创造特征，不用先考虑哪些特征可能好，可能不好，先弥补这个广度，所以我们不仅只做一些简单的加减组合，我们还要提取温度的统计特征。

# 加热过程
df['P1_S1_A6_0C'] = raw['A6']   # 容器初始温度
df['P1_S2_A8_1C'] = raw['A8']  # 首次测温温度
df['P1_S3_A10_2C'] = raw['A10']  # 准备水解温度
df['P1_C1_C0_D'] = raw['A8'] - raw['A6']   # 测温温差
df['P1_C2_C0_D'] = raw['A10'] - raw['A6']  # 初次沸腾温差

# 水解过程
df['P2_S1_A12_3C'] = raw['A12']  # 水解开始温度
df['P2_S2_A15_4C'] = raw['A15']  # 水解过程测温温度
df['P2_S3_A17_5C'] = raw['A17'] # 水解结束温度
df['P2_C3_C0_D'] = raw['A12'] - raw['A6']  # 水解开始与初始温度温差
df['P2_C3_C2_D'] = raw['A12'] - raw['A10']  # 水解开始前恒温温差
df['P2_C4_C3_D'] = raw['A15'] - raw['A12']  # 水解过程中途温差
df['P2_C5_C4_D'] = raw['A17'] - raw['A15']  # 水解结束中途温差
df['P2_C5_C3_KD'] = raw['A17'] - raw['A12']  # 水解起止温差

# 脱色过程
df['P3_S2_A25_7C'] = raw['A25']  # 脱色保温开始温度
df['P3_S3_A27_8C'] = raw['A27']  # 脱色保温结束温度
df['P3_C7_C5_D'] = raw['A25'] - raw['A17']  # 降温温差
df['P3_C8_C7_KD'] = raw['A27'] - raw['A25']  # 保温温差

# 结晶过程
df['P4_S2_B6_11C'] = raw['B6']  # 结晶开始温度
df['P4_S3_B8_12C'] = raw['B8']  # 结晶结束温度
df['P4_C11_C8_D'] = raw['B6'] - raw['A27']  # 脱色结束到结晶温差
df['P4_C12_C11_KD'] = raw['B8'] - raw['B6']  # 结晶温差

温度相关的统计特征， 比如不同工序里面类似过程每个样本做一些统计的特征等。

_funcs = ['mean', 'std', 'sum']
# 遍历每一种统计指标
for _func in _funcs:
    # 对每一个样本计算各项指标
    df[f'P2_C2-C5_{_func}'] = raw[['A10', 'A12', 'A15', 'A17']].agg(_func, axis=1)   # 沸腾过程温度
    df[f'P2_D3-D5_{_func}'] = df[[f'P2_C{i}_C{i-1}_D' for i in range(3, 6)]].abs().agg(_func, axis=1) # 沸腾过程绝对温差
    df[f'P2_C1-C12_KD_ABS_{_func}'] = df[[_f for _f in df.columns if _f.endswith('KD')]].abs().agg(_func, axis=1)  # 关键过程绝对温差
    df[f'P2_C1-C12_D_{_func}'] = df[[_f for _f in df.columns if _f.endswith('D')]].abs().agg(_func, axis=1)  # 所有过程绝对温差
    df[f'P2_LARGE_KD_{_func}'] = df[['P2_C3_C0_D', 'P3_C7_C5_D', 'P4_C12_C11_KD']].abs().agg(_func, axis=1)  # 大温差绝对温差

最后得到温度的相关特征：

"""得到温度相关特征"""
df_temperature = df.set_index('样本id')

df_temperature.head()   # 这样后期利于拼接

看一下做的特征：

关于时间特征的提取，水耗特征的提取思路和温度的一样，也是先单独加减组合，然后提取统计特征，在这里就不描述了， 代码太多， 具体的notebook可以见后面的链接。

3.3 合并所有特征

由于上面提取特征时都把索引重置成一样的了，这里合并的时候，直接相连，然后再重置索引即可。

df_feature = pd.concat([df_materials, df_duration, df_temperature, df_interact], axis=1).reset_index()

# 这样总共是做了143个特征出来
df_feature.head()

然后，我们记得把训练集和测试集分开：

df_trn = df_feature.iloc[:len(df_trn)].reset_index(drop=True)
df_trn['收率'] = df_target
df_tst = df_feature.iloc[len(df_trn):].reset_index(drop=True)
df_tst['收率'] = np.nan

3.4 筛选常规数据

import matplotlib.pyplot as plt

df_trn['收率'].plot(kind='hist')
plt.xlim(0.8, 1)
plt.show()

这个条形图可能看不出来：

我们下面换箱型图就会一目了然：

df_trn.boxplot(['收率'])

就会发现有几个离群点，我们使用query过滤掉。

# 删除离群点
df_trn = df_trn.query('收率 > 0.8671').reset_index(drop=True)
df_trn = df_trn.query('收率 < 0.9861').reset_index(drop=True)

df_trn.boxplot(['收率'])

结果如下：

4. 模型的建立和训练

本次模型使用xgboost， 关于xgboost的详细原理和应用，可以参考我的另一篇博客：白话机器学习算法理论+实战番外篇之Xgboost算法

def xgb_cv(train, test, params, fit_params, feature_names, nfold, seed):
    # 创建结果df
    train_pred = pd.DataFrame({
        'id': train['样本id'],
        'true': train['收率'],
        'pred': np.zeros(len(train))
    })
    
    # 测试提交结果
    test_pred = pd.DataFrame({'id':test['样本id'], 'pred':np.zeros(len(test))})
    # 交叉验证
    kfolder = KFold(n_splits=nfold, shuffle=True, random_state=seed)
    # 构造测试DMatrix
    xgb_tst = xgb.DMatrix(data=test[feature_names])
    print('\n')
    
    # 遍历cv的每一折数据，通过索引来指定
    for fold_id, (trn_idx, val_idx) in enumerate(kfolder.split(train['收率'])):
        # 构造当前训练的DMatrix
        xgb_trn = xgb.DMatrix(
            train.iloc[trn_idx][feature_names],
            train.iloc[trn_idx]['收率']
        )
        # 构造当前验证的DMatrix
        xgb_val = xgb.DMatrix(
            train.iloc[val_idx][feature_names],
            train.iloc[val_idx]['收率']
        )
        # 训练回归模型
        xgb_reg = xgb.train(params=params, dtrain=xgb_trn, **fit_params,
                           evals=[(xgb_trn, 'train'), (xgb_val, 'valid')])
        # 得到验证结果
        val_pred = xgb_reg.predict(
            xgb.DMatrix(train.iloc[val_idx][feature_names]),
            ntree_limit = xgb_reg.best_ntree_limit
        )
        train_pred.loc[val_idx, 'pred'] = val_pred
        # 这里是k折求了个平均
        test_pred['pred'] += xgb_reg.predict(xgb_tst, ntree_limit=xgb_reg.best_ntree_limit) / nfold
    print('\nCV LOSS: ', mse(train_pred['true'], train_pred['pred']), '\n')
    return test_pred

设置训练参数：

fit_params = {'num_boost_round':10800,
              'verbose_eval': 300,
              'early_stopping_rounds': 360
             }
params_xgb = {'eta':0.01, 'max_depth': 7, 'subsample':0.8,
              'booster':'gbtree', 'colsample_bytree':0.8,
              'objective':'reg:linear', 'silent':True, 'nthread':4
             }

开始训练：

# 开始训练
pred_xgb_a = xgb_cv(df_trn, df_tst, params_xgb, fit_params, df_trn.columns.tolist()[1:-1], 5, 0)

由于这是一个冠军的解决方案，所以可以看到效果还是很好的：

其实，感觉做比赛大部分时间都花在了特征工程上面， 模型这块估计都是xgboost这些主流的算法模型，之所以人家做的好，因为人家特征这块下了功夫。

5. 模型的测试和结果：

# 得到预测结果
df_tst_a['收率'] = pred_xgb_a['pred'].values
df_tst_a.head()

这个就是最后提交的形式， 最后结果看一下：

最后，每个样本对应着一个预测收益。

6. 知识点和思路的总结

关于知识点这块，这一次主要是整理学习pandas的数值访问，包括.loc, .iloc, .ix的区别，然后是pandas的set_index和reset_index的学习，然后是pandas的query过滤。这些我都整理到了另一篇博客：Pandas的数值访问(.loc, .iloc, .ix访问数据的区别)+query+set_index和reset_index

其次是xgboost的使用：白话机器学习算法理论+实战番外篇之Xgboost算法

最后是这个比赛的思路总结和技巧总结，这个比赛的思路其实并不是太复杂，相比较于上一次的快手用户活跃度预测来说，这个工业化生产预测可能更加符合常规，但是从这里面我们还是需要学习一些数据处理的tricks， 我整理如下：

1. 拿到数据之后，一定要从质疑的态度去审视数据，看看有没有个别的数值不合法（这里可以借助Excel的套用表格样式，或者自己写一个python函数遍历进行检查审视）
   
2. 数据预处理的时候，标准化归一化这种，最好是把测试集和训练集放在一块进行处理，这样可以更方便，前提得先统一维度，即把训练集的标签单独存放，然后再两者进行合并处理，合并的时候要记录训练集样本个数。处理完之后，再分开。
   
3. 进行数据处理的时候，最好是不要在原数据上进行处理，写成函数，然后进行好备份操作。
   
4. 做特征工程的时候，如果是从不同的角度构建特征，我们最好是每个角度写成函数创建一个空的DataFrame单独构建， 最后利用统一index进行合并
   
5. 关于数据的时间字段，处理的时候最好是给重命名一下，这样方便后面特殊处理的时候好找，或者还可以利用to_datetime转成标准的时间，这样这个字段就可以做很多事情了，也可以设置成索引等好多处理方式。
   
6. 时间段字段，可以从上面重新构建新特征，比如起始时间，终止时间，持续时间等。如果计算持续时间的时候，普通的表示方式计算机并不认识，这时候我们可以考虑转成分钟或者小时的方式。
   
7. 做特征工程的时候，我们需要发挥想象力，尽可能多的创造特征，不用先考虑哪些特征可能好，可能不好，先从弥补广度出发，做出来尽可能多的，不好的特征我们可以再后期删除嘛。 所以我们不仅只做一些简单的加减组合，我们还要考虑提取一些统计的特征。当然，这是连续的数值，如果是离散的数据，后面也会有相应的方式。
   
8. 筛选不常规的数据，我们观察数据的时候，也可以进行数据可视化，筛除掉一些离群的数据样本先。
   

关于这个比赛的详细代码和笔记，我已经放到了GitHub上，可以自行查阅， 如果感觉有帮助，欢迎star 😉
https://github.com/zhongqiangwu960812/AIGame/tree/master/GameOfTianChi/工业化工生产预测