阿里天池--工业蒸汽预测

赛题描述：

经脱敏后的锅炉传感器采集的数据（采集频率是分钟级别），根据锅炉的工况，预测产生的蒸汽量。

数据说明:
数据分成训练数据（train.txt）和测试数据（test.txt），其中字段 V0-V37，这38个字段是作为特征变量，target作为目标变量。 选手利用训练数据训练出模型，预测测试数据的目标变量，排名结果依据预测结果的 MSE（mean square error）。

结果提交：
选手需要提交测试数据的预测结果（txt格式，只有1列预测结果）。

结果评估：
预测结果以mean square error作为评判标准。

1  数据分析

导入相关库

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
sns.set(style='ticks')

from scipy import stats
from sklearn.decomposition import PCA
from sklearn.metrics import mean_squared_error, r2_score
from sklearn.model_selection import train_test_split, GridSearchCV, learning_curve, RepeatedKFold

from sklearn.linear_model import LinearRegression, RidgeCV, LassoCV
from sklearn.svm import SVR
from xgboost import XGBRegressor

import warnings
warnings.filterwarnings('ignore')

导入数据并打印

dtrain = pd.read_csv('zhengqi_train.txt', sep='\t')
dtest = pd.read_csv('zhengqi_test.txt', sep='\t')
dfull = pd.concat([dtrain, dtest], ignore_index=True, sort=False)
print('训练集大小: ', np.shape(dtrain))
print('测试集大小: ', np.shape(dtest))

训练集大小: (2888, 39)
测试集大小: (1925, 38)

 分析数据的基本分布情况：

plt.figure(figsize=(18,8),dpi=100)
dfull.boxplot(sym='r^', patch_artist=True, notch=True)
plt.title('DATA-FULL')

• 数据共有38个特征变量和1目标变量
• 所有特征变量及蒸汽量目标变量均为连续型数据
• 所有数据均分布在0附近，且无量级差异

  # 查看数据缺失情况
  
  print('缺失值统计：')
  print(dfull.apply(lambda x: sum(x.isnull())))
  

  缺失值统计：
  V0 0
  V1 0
  V2 0
  V3 0
  V4 0
  …
  V34 0
  V35 0
  V36 0
  V37 0
  target 1925
  dtype: int64

 并没有确实的数据。

接下来绘制数据相关性热力图，查看数据相关性情况：

def heatmap(df):
    plt.figure(figsize=(20,16), dpi=100)
    cols = df.columns.tolist()
    mcorr = df[cols].corr(method = 'spearman')
    cmap = sns.diverging_palette(220, 10, as_cmap=True)
    mask = np.zeros_like(mcorr, dtype = np.bool)
    mask[np.triu_indices_from(mask)] = True              # 角分线右侧为True
    g = sns.heatmap(mcorr, mask=mask, cmap=cmap, square=True, annot=True, fmt='0.2f')
    plt.xticks(rotation=45)
    return mcorr
    
dtrain_mcorr = heatmap(dtrain)

 热力图中发现，特征之间存在多重共线性，针对这一问题处理方法有以下三种：

• 主成分分析（PCA）：通常情况下可能导致模型精度降低
• 正则化：许多较为成熟模型中都带有正则项，同时能够防止模型过拟合
• 逐步回归法：较为繁琐，不建议

2 特征工程

 特征工程往往是影响预测结果好坏的关键所在，在此项目中我主要对数据进行了以下处理：

• 剔除训练集与测试集分布不均匀的特征变量
• 查看特征与目标变量相关系数，剔除相关系数低于0.1的无关变量
• 查看特征分布偏态系数，将部分明显左偏/右偏的特征数据进行正态化
• 对剩余特征数据进行标准化

创建特征剔除函数，方便对数据进行统一处理

# 删除dfull表中指定的列
def drop_var(var_lst):
    dfull.drop(var_lst, axis=1, inplace=True)
 
# 将dfull重新分割为dtrain核dtest   
def split_dfull():
    dtrain = dfull[dfull['target'].notnull()]
    dtest = dfull[dfull['target'].isnull()]
    dtest.drop('target', axis=1, inplace=True)
    return dtrain, dtest

1. 剔除训练集与测试集分布不均匀的特征变量

针对38个特征变量画训练数据和测试数据特征分布直方图

plt.figure(figsize=(20,50),dpi=100)
for i in range(38):
    plt.subplot(10,4,i+1)
    sns.distplot(dtrain.iloc[:,i], color='green')
    sns.distplot(dtest.iloc[:,i], color='red')
    plt.legend(['Train','Test'])
plt.tight_layout()

 从各特征数据分布直方图发现：
V5、V9、V11、V17、V20、V21、V22、V27、V28 特征训练集和测试集分布差异过大
因此为了减小数据分布不均对预测结果的影响，应将上述特征进行剔除

drop_var(['V5','V9','V11','V17','V20','V21','V22','V27','V28'])
dtrain, dtest = split_dfull()

2. 查看各特征与结果变量的相关性，以及是否服从正态分布

利用sns.regplot查看特征与目标变量间的线性关系强弱；
利用sns.distplot查看每个特征分布情况，并进行正态拟合。

plt.figure(figsize=(20,60),dpi=80)
i = 0
for col in dtrain.columns:
    i += 1
    plt.subplot(15,4,i)
    sns.regplot(col, 'target', data = dtrain,
               scatter_kws = {'marker': '.', 's': 5, 'alpha': .6},
               line_kws = {'color': 'k'})
    coef = np.corrcoef(dtrain[col], dtrain['target'])[0][1]
    plt.title(f'coef = {coef}')
    plt.xlabel(col)
    plt.ylabel('target')
    
    i += 1
    plt.subplot(15,4,i)
    sns.distplot(dtrain[col], fit = stats.norm)
    plt.title(f'skew = {stats.skew(dtrain[col])}')
    plt.xlabel(col)
    
plt.tight_layout()

从图中可看出：

• V14、V25、V26、V32、V33、V34、V35 与结果变量相关性很低，予以 删除
• V0、V1、V6、V7、V8、V12、V16、V31 存在明显的数据左偏，可用 对数转换 降低偏态

 3. 删除无关变量

drop_var(['V14','V25','V26','V32','V33','V34','V35'])
dtrain, dtest = split_dfull()

4. 对偏态数据进行正态化转换

# 分布呈明显左偏的特征
piantai = ['V0','V1','V6','V7','V8','V12','V16','V31']

# 创建函数——找到令偏态系数绝对值最小的对数转换的底
def find_min_skew(data):
    subs = list(np.arange(1.01,2,0.01))
    skews = []
    for x in subs:
        skew = abs(stats.skew(np.power(x,data)))
        skews.append(skew)
    min_skew = min(skews)
    i = skews.index(min_skew)
    return subs[i], min_skew

偏态数据对数转换前后对比可视化

plt.figure(figsize=(20,36), dpi=80)

i = 0
for col in piantai:
    
    # 正态化之前的偏态特征与结果变量相关图与正态分布图
    i += 1
    plt.subplot(8,4,i)
    sns.regplot(col, 'target', data = dtrain,
               scatter_kws = {'marker': '.', 's': 5, 'alpha': .6},
               line_kws = {'color': 'k'})
    coef = np.corrcoef(dtrain[col], dtrain['target'])[0][1]
    plt.title(f'coef = {coef}')
    plt.xlabel(col)
    plt.ylabel('target')
    
    i += 1
    plt.subplot(8,4,i)
    sns.distplot(dtrain[col], fit = stats.norm)
    plt.title(f'skew = {stats.skew(dtrain[col])}')
    plt.xlabel(col)
    
    # 找到合适的底数，对dtrain的偏态特征进行对数转换
    sub = find_min_skew(dfull[col])[0]
    dtrain[col] = np.power(sub, dtrain[col])
    
    # 对数转换之后特征与结果变量相关图与正态分布图
    i += 1
    plt.subplot(8,4,i)
    sns.regplot(col, 'target', data = dtrain,
               scatter_kws = {'marker': '.', 's': 5, 'alpha': .6},
               line_kws = {'color': 'k'})
    coef = np.corrcoef(dtrain[col], dtrain['target'])[0][1]
    plt.title(f'coef = {coef}')
    plt.xlabel(col+' transformed')
    plt.ylabel('target')
    
    i += 1
    plt.subplot(8,4,i)
    sns.distplot(dtrain[col], fit = stats.norm)
    plt.title(f'skew = {stats.skew(dtrain[col])}')
    plt.xlabel(col+' transformed')
    
plt.tight_layout()

 偏态系数 明显减小，同时相关系数未有明显影响
证明该方法可以有效实现特征正态化

因此，对训练集和测试集偏态特征同时进行对数转换

for col in piantai:
    sub = find_min_skew(dfull[col])[0]
    dfull[col] = np.power(sub, dfull[col])
dtrain, dtest = split_dfull()

5. 对剩余特征数据进行标准化

采用 z-score标准化 方法 

dfull.iloc[:,:-1] = dfull.iloc[:,:-1].apply(lambda x: (x-x.mean())/x.std())
dtrain, dtest = split_dfull()

 特征工程完成！
查看完成后特征分布直方图

plt.figure(figsize=(20,20),dpi=80)
for i in range(22):
    plt.subplot(6,4,i+1)
    sns.distplot(dtrain.iloc[:,i], color='green')
    sns.distplot(dtest.iloc[:,i], color='red')
    plt.legend(['Train','Test'])
plt.tight_layout()

 3 模型训练

模型训练之前，先用岭回归进行预测，剔除训练集中的异常样本（预测残差在3sigma之外的点）
然后使用以下回归模型进行 模型训练 和 参数调优

• Lasso 回归
• SVR 回归
• XGBoost 回归

将训练数据分割为训练集与验证集 

X = np.array(dtrain.iloc[:,:-1])
y = np.array(dtrain['target'])

X_train, X_valid, y_train, y_valid = train_test_split(X, y, test_size=0.2, random_state=0)

 创建模型评分函数

def score(y, y_pred):
	# 计算均方误差 MSE
    print('MSE = {0}'.format(mean_squared_error(y, y_pred)))
    # 计算模型决定系数 R2
    print('R2 = {0}'.format(r2_score(y, y_pred)))
    
    # 计算预测残差，找异常点
    y = pd.Series(y)
    y_pred = pd.Series(y_pred, index=y.index)
    resid = y - y_pred
    mean_resid = resid.mean()
    std_resid = resid.std()
    z = (resid - mean_resid) / std_resid
    n_outliers = sum(abs(z)>3)
    
    # 图一：真实值vs预计值
    plt.figure(figsize=(18,5), dpi=80)
    plt.subplot(131)
    plt.plot(y, y_pred, '.')
    plt.xlabel('y')
    plt.ylabel('y_pred')
    plt.title('corr = {:.3f}'.format(np.corrcoef(y,y_pred)[0][1]))
    
    # 图二：残差分布散点图
    plt.subplot(132)
    plt.plot(y, y-y_pred, '.')
    plt.xlabel('y')
    plt.ylabel('resid')
    plt.ylim([-3,3])
    plt.title('std_resid = {:.3f}'.format(std_resid))
    
    # 图三：残差z得分直方图
    plt.subplot(133)
    sns.distplot(z, bins=50)
    plt.xlabel('z')
    plt.title('{:.0f} samples with z>3'.format(n_outliers))
    plt.tight_layout()

开始模型训练前，利用岭回归模型预测，剔除异常样本

# 利用RidgeCV函数自动寻找最优参数
ridge = RidgeCV()
ridge.fit(X_train, y_train)
print('best_alpha = {0}'.format(ridge.alpha_))

y_pred = ridge.predict(X_train)
score(y_train, y_pred)

best_alpha = 0.1
MSE = 0.1058302793365256
R2 = 0.888601902633553

 找出异常样本点并剔除

resid = y_train - y_pred
resid = pd.Series(resid, index=range(len(y_train)))
resid_z = (resid-resid.mean()) / resid.std()
outliers = resid_z[abs(resid_z)>3].index
print(f'{len(outliers)} Outliers:')
print(outliers.tolist())

plt.figure(figsize=(14,6),dpi=60)

plt.subplot(121)
plt.plot(y_train, y_pred, '.')
plt.plot(y_train[outliers], y_pred[outliers], 'ro')
plt.title(f'MSE = {mean_squared_error(y_train,y_pred)}')
plt.legend(['Accepted', 'Outliers'])
plt.xlabel('y_train')
plt.ylabel('y_pred')

plt.subplot(122)
sns.distplot(resid_z, bins = 50)
sns.distplot(resid_z.loc[outliers], bins = 50, color = 'r')
plt.legend(['Accepted', 'Outliers'])
plt.xlabel('z')
plt.tight_layout()

27 Outliers:
[89, 692, 768, 797, 829, 1123, 1190, 1198, 1236, 1295, 1324, 1366, 1530, 1544, 1641, 1706, 1718, 1796, 1839, 1983, 2013, 2022, 2025, 2227, 2254, 2279, 2299]

# 异常样本点剔除
X_train = np.array(pd.DataFrame(X_train).drop(outliers,axis=0))
y_train = np.array(pd.Series(y_train).drop(outliers,axis=0))

开始进行模型训练

 模型1: Lasso回归

利用LassoCV自动选择最佳正则化参数

lasso = LassoCV(cv=5)
lasso.fit(X_train, y_train)
print('best_alpha = {0}'.format(lasso.alpha_))

pred_lasso = lasso.predict(X_valid)
score(y_valid, pred_lasso)

best_alpha = 0.0007812269286266984
MSE = 0.11780207761554924
R2 = 0.8866338406111323​​​​​​​

 模型2: 支持向量回归（SVR）

模型参数：

•  kernel [default=‘rbf’]:

        指定要在算法中使用的内核类型。它必是’linear’，‘poly’，‘rbf’，‘sigmoid’，‘precomputed’或者callable之一。
如果没有给出，将使用’rbf’。如果给出了callable，则它用于预先计算内核矩阵。

• C [default=0]:

        错误术语的惩罚参数C.

• gamma [default=‘auto’]:

        ‘rbf’，‘poly’和’sigmoid’的核系数。默认值为’auto’，它使用1/nfeatures，如果gamma='scale’传递，则使用1/nfeatures∗X .std作为gamma的值。

使用sklearn中的网格搜索方法 GridSearchCV 寻找SVR最优模型参数
创建GridSearchCV网格参数搜寻函数，评价标准为最小均方误差，采用K折交叉验证的检验方法

def gsearch(model, param_grid, scoring='neg_mean_squared_error', splits=5, repeats=1, n_jobs=-1):
    # p次k折交叉验证
    rkfold = RepeatedKFold(n_splits=splits, n_repeats=repeats, random_state=0)
    model_gs = GridSearchCV(model, param_grid=param_grid, scoring=scoring, cv=rkfold, verbose=1, n_jobs=-1)
    model_gs.fit(X_train, y_train)
    print('参数最佳取值: {0}'.format(model_gs.best_params_))
    print('最小均方误差: {0}'.format(abs(model_gs.best_score_)))
    return model_gs

使用SVR回归器默认的“rbf”内核，即高斯核
对惩罚参数C与核系数gamma进行网格搜索CV验证

svr = SVR()
cv_params = {'C': np.logspace(0, 3, 4), 'gamma': np.logspace(-4, -1, 4)}
svr = gsearch(svr, cv_params)

参数最佳取值: {'C': 100.0, 'gamma': 0.001}
最小均方误差: 0.08739677966938571

 缩小参数范围进行细调：

svr = SVR()
cv_params = {'C': [1,2,5,10,15,20,30,50,80,100,150,200], 'gamma': [0.0001,0.0005,0.0008,0.001,0.002,0.003,0.005]}
svr = gsearch(svr, cv_params)

参数最佳取值: {'C': 150, 'gamma': 0.0008}
最小均方误差: 0.0872966704265247

 确定最佳参数（C:15, gamma：0.0008）
验证集预测

pred_svr = svr.predict(X_valid)
score(y_valid, pred_svr)

MSE = 0.11160923124447547
R2 = 0.8925934910942511​​​​​​​

模型3: XGB回归（XGBRegressor ）

模型参数：

• n_estimators [default=100]:
  决策树的棵树，即总共迭代的次数
• learning_rate [default=0.3]:
  shrinkage参数，用于更新叶子节点权重时，乘以该系数，避免步长过大。参数值越大，越可能无法收敛。
  把学习率 eta 设置的小一些，小学习率可以使得后面的学习更加仔细。
• min_child_weight [default=1]:
  这个参数默认是 1，是每个叶子里面 h 的和至少是多少，对正负样本不均衡时的 0-1 分类而言，假设 h 在 0.01 附近，min_child_weight 为 1 意味着叶子节点中最少需要包含 100 个样本。这个参数非常影响结果，控制叶子节点中二阶导的和的最小值，该参数值越小，越容易overfitting。
• max_depth [default=6]:
  每颗树的最大深度，树高越深，越容易过拟合。
• max_leaf_nodes:
  最大叶结点数，与max_depth作用有点重合。
• gamma [default=0]:
  后剪枝时，用于控制是否后剪枝的参数。
• max_delta_step [default=0]:
  这个参数在更新步骤中起作用，如果取0表示没有约束，如果取正值则使得更新步骤更加保守。可以防止做太大的更新步子，使更新更加平缓。
• subsample [default=1]:
  样本随机采样，较低的值使得算法更加保守，防止过拟合，但是太小的值也会造成欠拟合。
• colsample_bytree [default=1]:
  列采样，对每棵树的生成用的特征进行列采样.一般设置为： 0.5-1
• reg_lambda [default=1]:
  控制模型复杂度的权重值的L2正则化项参数，参数越大，模型越不容易过拟合。
• reg_alpha [default=0]:
  控制模型复杂程度的权重值的 L1 正则项参数，参数值越大，模型越不容易过拟合。
• scale_pos_weight [default=1]:
  如果取值大于0的话，在类别样本不平衡的情况下有助于快速收敛

# 初始参数值
params = {'learning_rate': 0.1, 'n_estimators': 500, 'max_depth': 5, 'min_child_weight': 1, 'seed': 0,
          'subsample': 0.8, 'colsample_bytree': 0.8, 'gamma': 0, 'reg_alpha': 0, 'reg_lambda': 1}

1. 最佳迭代次数：n_estimators

cv_params = {'n_estimators': [100,200,300,400,500,600,700,800,900,1000,1100,1200]}
xgb = XGBRegressor(**params)
xgb = gsearch(xgb, cv_params)

参数最佳取值: {'n_estimators': 800}
最小均方误差: 0.10106838497021398

# 更新参数
params['n_estimators'] = 800

 2. min_child_weight 以及 max_depth

cv_params = {'max_depth': [3,4,5,6,7,8,9],
            'min_child_weight': [1,2,3,4,5,6,7]}
xgb = XGBRegressor(**params)
xgb = gsearch(xgb, cv_params)

参数最佳取值: {'max_depth': 5, 'min_child_weight': 3}
最小均方误差: 0.09790308113659893

# 更新参数  
params['max_depth'] = 5
params['min_child_weight'] = 3

3. 后剪枝参数 gamma

cv_params = {'gamma': [0,0.01,0.05,0.1,0.2,0.3,0.4,0.5,0.6]}
xgb = XGBRegressor(**params)
xgb = gsearch(xgb, cv_params)

参数最佳取值: {'gamma': 0.1}
最小均方误差: 0.09737878476012207

# 更新参数
params['gamma'] = 0

 4. 样本采样subsample 和 列采样colsample_bytree

cv_params = {'subsample': [0.6,0.7,0.8,0.9],
            'colsample_bytree': [0.6,0.7,0.8,0.9]}
xgb = XGBRegressor(**params)
xgb = gsearch(xgb, cv_params)

参数最佳取值: {'colsample_bytree': 0.6, 'subsample': 0.8}
最小均方误差: 0.09605732520223123

# 更新参数
params['subsample'] = 0.8
params['colsample_bytree'] = 0.8

5. L1正则项参数reg_alpha 和 L2正则项参数reg_lambda

cv_params = {'reg_alpha': [0,0.02,0.05,0.1,1,2,3],
             'reg_lambda': [0,0.02,0.05,0.1,1,2,3]}
xgb = XGBRegressor(**params)
xgb = gsearch(xgb, cv_params)

参数最佳取值: {‘reg_alpha’: 0, ‘reg_lambda’: 1}
最小均方误差: 0.09550672795811381

# 更新参数
params['reg_alpha'] = 0
params['reg_lambda'] = 1

6.最后是learning_rate，一般这时候要调小学习率来测试

cv_params = {'learning_rate': [0.01, 0.02, 0.03, 0.04, 0.05, 0.07, 0.1, 0.2]}
xgb = XGBRegressor(**params)
xgb = gsearch(xgb, cv_params)

参数最佳取值: {‘learning_rate’: 0.04}
最小均方误差: 0.09450026627677707

# 更新参数
params['learning_rate'] = 0.04

参数调优完成，以验证集进行模型误差验证

pred_xgb = xgb.predict(X_valid)
score(y_valid, pred_xgb)

MSE = 0.11754469565827244
R2 = 0.8868815306738465​​​​​​​

4 模型评估

 学习曲线 (learning_curve)

在训练模型时，经常我们是需要判断训练出来的模型是否是过拟合或者欠拟合的，这个时候就需要借助Learning curve来帮助我们进行判断，Learning curve描述的是在不同样本的情况下，训练集和验证集的准确率的变化曲线，通过判断两条曲线的间隔就可以看出模型的大致状态。

models = [lasso, svr, xgb]
model_names = ['Lasso','SVR','XGB']
plt.figure(figsize=(20,5))

for i,m in enumerate(models):
    train_sizes, train_scores, test_scores = learning_curve(m, X, y, cv=5, scoring='neg_mean_squared_error',
                                                            train_sizes=np.linspace(0.1,1.0,5), n_jobs=-1)
    train_scores_mean = -train_scores.mean(axis=1)
    test_scores_mean = -test_scores.mean(axis=1)
    
    plt.subplot(1,3,i+1)
    plt.plot(train_sizes, train_scores_mean, 'o-', label='Train')
    plt.plot(train_sizes, test_scores_mean, '^-', label='Test')
    plt.xlabel('Train_size')
    plt.ylabel('Score')
    plt.ylim([0,0.35])
    plt.title(model_names[i], fontsize=16)
    plt.legend()
    plt.grid()
    
plt.tight_layout()

 

 从学习曲线可看出，Lasso与SVR模型效果还行，XGB回归模型存在轻微过拟合，后续有待进一步参数优化进行改进。

模型加权融合

对于多个模型结果，采取加权融合的办法进行结合，即对各模型预测结果取加权平均，这样可用避免单个模型在预测某一部分数据时产生较大的误差。
为了取得最小的MSE，利用验证集对各种加权方法进行验证，确定三个模型的权重。

def model_mix(pred_1, pred_2, pred_3):
    result = pd.DataFrame(columns=['Lasso','SVR','XGB','Combine'])
    for a in range(5):
        for b in range(1,6):
            for c in range(5):
                y_pred = (a*pred_1 + b*pred_2 + c*pred_3) / (a+b+c)
                mse = mean_squared_error(y_valid, y_pred)
                result = result.append([{'Lasso':a, 'SVR':b, 'XGB':c, 'Combine':mse}], ignore_index=True)
    return result

model_combine = model_mix(pred_lasso, pred_svr, pred_xgb)
model_combine.sort_values(by='Combine', inplace=True)
model_combine.head()

 结果表明：在验证集上，Lasso、SVR、XGBRegressor三种模型分别取权重为1/8、8/5、1/4时得到的预测数据均方误差最小，线下 Min(MSE) = 0.1121

5 模型预测

对3种模型预测结果进行加权融合

X_test = np.array(dtest)
ans_lasso = lasso.predict(X_test)
ans_svr = svr.predict(X_test)
ans_xgb = xgb.predict(X_test)
ans_mix = (ans_lasso + 5 * ans_svr + 2 * ans_xgb) / 8
pd.Series(ans_mix).to_csv('正态+标准化.txt', sep='\t', index=False)
print('Finished!')

Finished!

文章参考   博主：乔木zzz