----------------------------数据分析-------------------------------

查看特征的数值类型有哪些，对象类型有哪些：

• 特征一般都是由类别型特征和数值型特征组成，而数值型特征又分为连续型和离散型。
• 类别型特征有时具有非数值关系，有时也具有数值关系。比如‘grade’中的等级A，B，C等，是否只是单纯的分类，还是A优于其他要结合业务判断。
• 数值型特征本是可以直接入模的，但往往风控人员要对其做分箱，转化为WOE编码进而做标准评分卡等操作。从模型效果上来看，特征分箱主要是为了降低变量的复杂性，减少变量噪音对模型的影响，提高自变量和因变量的相关度。从而使模型更加稳定。

numerical_fea = list(data_train.select_dtypes(exclude=['object']).columns)
category_fea = list(filter(lambda x: x not in numerical_fea,list(data_train.columns)))

划分数值型变量中的连续变量和离散型变量：

#过滤数值型类别特征
def get_numerical_serial_fea(data,feas):
    numerical_serial_fea = []
    numerical_noserial_fea = []
    for fea in feas:
        temp = data[fea].nunique()
        #取值数少于10的类别认为数值型类别特征
        if temp <= 10:
            numerical_noserial_fea.append(fea)
            continue
        numerical_serial_fea.append(fea)
    return numerical_serial_fea,numerical_noserial_fea
numerical_serial_fea,numerical_noserial_fea = get_numerical_serial_fea(data_train,numerical_fea)

数值连续型变量可视化分析：

#每个数字特征得分布可视化
f = pd.melt(data_train, value_vars=numerical_serial_fea)
""
col表示值的名称
col_wrap表示一行显示2个子图
sharex、sharey表示是否恭喜x、y轴
""
g = sns.FacetGrid(f, col="variable",  col_wrap=2, sharex=False, sharey=False)
g = g.map(sns.distplot, "value")

g = sns.FacetGrid(tips, col="sex", hue="smoker")
g.map(plt.scatter,"total_bill", "tip", alpha=7)
g.add_legend()

iris = sns.load_dataset('iris')
g = sns.PairGrid(iris)
g.map(plt.scatter)

g = sns.PairGrid(iris,hue="species")
#对角线上显示的图形的类别
g.map_diag(plt.hist)
#非对角线上显示图形的类别
g.map_offdiag(plt.scatter)
#显示图例
g.add_legend()

查看数据集中特征缺失值，唯一值等:

缺失值：

print(f'There are {data_train.isnull().any().sum()} columns in train dataset with missing values.')

上面得到训练集有22列特征有缺失值，进一步查看缺失特征中缺失率大于50%的特征

have_null_fea_dict = (data_train.isnull().sum()/len(data_train)).to_dict()
fea_null_moreThanHalf = {}
for key,value in have_null_fea_dict.items():
    if value > 0.5:
        fea_null_moreThanHalf[key] = value

 查看数据集中特征属性只有一值的特征：

one_value_fea = [col for col in data_train.columns if data_train[col].nunique() <= 1]

可视化缺失特征及缺失率：

# nan可视化
missing = data_train.isnull().sum()/len(data_train)
missing = missing[missing > 0]
missing.sort_values(inplace=True)
missing.plot.bar()

• 纵向了解哪些列存在 “nan”, 并可以把nan的个数打印，主要的目的在于查看某一列nan存在的个数是否真的很大，如果nan存在的过多，说明这一列对label的影响几乎不起作用了，可以考虑删掉。如果缺失值很小一般可以选择填充。
• 另外可以横向比较，如果在数据集中，某些样本数据的大部分列都是缺失的且样本足够的情况下可以考虑删除。

Tips: 比赛大杀器lgb模型可以自动处理缺失值。

重复样本检测：

def duplicate_det(data):
    if data.duplicated(subset=['QUEUE_ID', 'CU', 'DOTTING_TIME']).any() != True:
        print("无重复无样本")
    else:
        print("有重复样本")
        print("重复样本数量:", len(data[data.duplicated(subset=['QUEUE_ID', 'CU', 'DOTTING_TIME'])]))

按天可视化每天的数据：

for key,item in df_train.groupby(by=['month','day']):
    df = pd.DataFrame(item)
    plt.figure()
    plt.title("data {}".format(key))
    plt.plot(df['temp_in'])
    plt.show()

msno可视化缺失值： 

msno.matrix(Train_data.sample(250))

 

msno.bar(Train_data.sample(1000))

 

可视化数值型变量及其log处理：

• 查看某一个数值型变量的分布，查看变量是否符合正态分布，如果不符合正太分布的变量可以log化后再观察下是否符合正态分布。
• 如果想统一处理一批数据变标准化 必须把这些之前已经正态化的数据提出。
• 正态化的原因：一些情况下正态非正态可以让模型更快的收敛，一些模型要求数据正态（eg. GMM、KNN）,保证数据不要过偏态即可，过于偏态可能会影响模型预测结果。

#Ploting Transaction Amount Values Distribution
plt.figure(figsize=(16,12))
plt.suptitle('Transaction Values Distribution', fontsize=22)
plt.subplot(221)
sub_plot_1 = sns.distplot(data_train['loanAmnt'])
sub_plot_1.set_title("loanAmnt Distribuition", fontsize=18)
sub_plot_1.set_xlabel("")
sub_plot_1.set_ylabel("Probability", fontsize=15)

plt.subplot(222)
sub_plot_2 = sns.distplot(np.log(data_train['loanAmnt']))
sub_plot_2.set_title("loanAmnt (Log) Distribuition", fontsize=18)
sub_plot_2.set_xlabel("")
sub_plot_2.set_ylabel("Probability", fontsize=15)

类别型变量分布可视化：

plt.figure(figsize=(8, 8))
sns.barplot(data_train["employmentLength"].value_counts(dropna=False)[:20],
            data_train["employmentLength"].value_counts(dropna=False).keys()[:20])
plt.show()

可视化类别型变量在不同y值上的分布：

train_loan_fr = data_train.loc[data_train['isDefault'] == 1]
train_loan_nofr = data_train.loc[data_train['isDefault'] == 0]

fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(15, 8))
train_loan_fr.groupby('grade')['grade'].count().plot(kind='barh', ax=ax1, title='Count of grade fraud')
train_loan_nofr.groupby('grade')['grade'].count().plot(kind='barh', ax=ax2, title='Count of grade non-fraud')
train_loan_fr.groupby('employmentLength')['employmentLength'].count().plot(kind='barh', ax=ax3, title='Count of employmentLength fraud')
train_loan_nofr.groupby('employmentLength')['employmentLength'].count().plot(kind='barh', ax=ax4, title='Count of employmentLength non-fraud')
plt.show()

可视化连续型变量在不同y值上的分布：

fig, ((ax1, ax2)) = plt.subplots(1, 2, figsize=(15, 6))
data_train.loc[data_train['isDefault'] == 1]['loanAmnt'].apply(np.log) \
    .plot(kind='hist',
          bins=100,
          title='Log Loan Amt - Fraud',
          color='r',
          xlim=(-3, 10),
         ax= ax1)
data_train.loc[data_train['isDefault'] == 0]['loanAmnt'].apply(np.log) \
    .plot(kind='hist',
          bins=100,
          title='Log Loan Amt - Not Fraud',
          color='b',
          xlim=(-3, 10),
         ax=ax2)

 时间格式数据处理及可视化：

#转化成时间格式  issueDateDT特征表示数据日期离数据集中日期最早的日期（2007-06-01）的天数
data_train['issueDate'] = pd.to_datetime(data_train['issueDate'],format='%Y-%m-%d')
startdate = datetime.datetime.strptime('2007-06-01', '%Y-%m-%d')
data_train['issueDateDT'] = data_train['issueDate'].apply(lambda x: x-startdate).dt.days

#转化成时间格式
data_test_a['issueDate'] = pd.to_datetime(data_train['issueDate'],format='%Y-%m-%d')
startdate = datetime.datetime.strptime('2007-06-01', '%Y-%m-%d')
data_test_a['issueDateDT'] = data_test_a['issueDate'].apply(lambda x: x-startdate).dt.days

plt.hist(data_train['issueDateDT'], label='train')
plt.hist(data_test_a['issueDateDT'], label='test')
plt.legend()
plt.title('Distribution of issueDateDT dates')
#train 和 test issueDateDT 日期有重叠 所以使用基于时间的分割进行验证是不明智的

运用确定值填补空缺值：

Train_data['notRepairedDamage'].replace('-', np.nan, inplace=True)

拟合Label可视化曲线：

sns库参考博客

## 1) 总体分布概况（无界约翰逊分布等）
import scipy.stats as st
y = Train_data['price']
plt.figure(1)
plt.title('Johnson SU')
sns.distplot(y, kde=False, fit=st.johnsonsu)
plt.figure(2)
plt.title('Normal')
sns.distplot(y, kde=False, fit=st.norm)
plt.figure(3)
plt.title('Log Normal')
sns.distplot(y, kde=False, fit=st.lognorm)

Label不服从正态分布，所以在进行回归之前，它必须进行转换。虽然对数变换做得很好，但最佳拟合是无界约翰逊分布

相关性分析及其可视化：

price_numeric = Train_data[numeric_features]
correlation = price_numeric.corr()
print(correlation['price'].sort_values(ascending = False),'\n')

f , ax = plt.subplots(figsize = (7, 7))

plt.title('Correlation of Numeric Features with Price',y=1,size=16)

sns.heatmap(correlation,square = True,  vmax=0.8)

参数：annot=True 表示将相关系数显示到表格上，

          fmt="d"指定字体的格式，否则会出现乱码。

查看数值特征的偏度和峰度：

# 2) 查看几个特征得 偏度和峰值
for col in numeric_features:
    print('{:15}'.format(col), 
          'Skewness: {:05.2f}'.format(Train_data[col].skew()) , 
          '   ' ,
          'Kurtosis: {:06.2f}'.format(Train_data[col].kurt())  
         )

Train & Test分布探索：

Train & Test分布主要是为了看数据的分布情况。

下面的方法是最近kaggle上GM搞的一种极其简单的探索训练集&测试集分布的方案,将一个分布的问题变成一个二分类问题,该方法简单有效，如果AUC低于0.6,则我们可以认为训练集和测试集是分布平衡的,反之我们则可以认为训练集和测试集是分布不一致的,这个时候就需要注意AB榜是否会出现翻车的情况等。

train['label'] = 0
test['label']  = 1
trte = pd.concat([train,test],axis=0,ignore_index=True)

import lightgbm as lgb
import time
from sklearn.model_selection import KFold,StratifiedKFold

def validation_prediction_lgb(X,y,feature_names, ratio =1, X_test = None,istest = False):
    n_fold = 5
    folds = StratifiedKFold(n_splits=n_fold, shuffle=True, random_state=42)

    params = {
    'bagging_freq': 5, 
    'boost_from_average':'false',
    'boost': 'gbdt',
    'learning_rate': 0.01,
    'max_depth': 5,
    'metric':'auc',
    'min_data_in_leaf': 50,
    'min_sum_hessian_in_leaf': 10.0,
    'tree_learner': 'serial',
    'objective': 'binary',
    'verbosity': 1}

    importances = pd.DataFrame() 
    if istest:
        prediction = np.zeros(len(X_test))
    models = []
    for fold_n, (train_index, valid_index) in enumerate(folds.split(X,y)):
        print('Fold', fold_n, 'started at', time.ctime())
        X_train, X_valid = X[train_index], X[valid_index]
        y_train, y_valid = y[train_index], y[valid_index]
        weights = [ratio  if val == 1 else 1 for val in y_train]
        
        train_data = lgb.Dataset(X_train, label=y_train,  weight=weights)
        valid_data = lgb.Dataset(X_valid, label=y_valid)

        model = lgb.train(params,train_data,num_boost_round=20000,
                        valid_sets = [train_data, valid_data],verbose_eval=200,early_stopping_rounds = 200)
        
        imp_df = pd.DataFrame() 
        imp_df['feature']  = feature_names
        imp_df['split']    = model.feature_importance()
        imp_df['gain']     = model.feature_importance(importance_type='gain')
        imp_df['fold']     = fold_n + 1
        
        importances = pd.concat([importances, imp_df], axis=0)
        
        models.append(model)
        if istest == True:
            prediction += model.predict(X_test, num_iteration=model.best_iteration)/5
    if istest == True:     
        return models,importances, prediction
    else:
        return models,importances

train_cols = [col for col in trte.columns if col not in ['target','ID_code','label']]
print(len(train_cols))
models,importances = validation_prediction_lgb(trte[train_cols].values, trte['label'].values, train_cols)

小提琴图可视化：

sns.violinplot(x="total_bill",y="day",hue="time",data=tips)

 上面time有2种取值，于是可视化出来的效果不太好，下面指定参数split=True:

sns.violinplot(x="day",y="total_bill",hue="sex",data=tips,split=True)

分类别直方图：

sns.barplot(x="sex",y="survived",hue="class",data=titanic)

查看特征的相对均匀分布：

sns.set()
columns = ['price', 'v_12', 'v_8' , 'v_0', 'power', 'v_5',  'v_2', 'v_6', 'v_1', 'v_14']
sns.pairplot(Train_data[columns],size = 2 ,kind ='scatter',diag_kind='kde')
plt.show()

查看类别特征的unique分布：

for fea in categorical_features:
    print(Train_data[fea].nunique())

将时间戳特征解析为对应时间：

for df in [train_df, test_df]:
    t = pd.to_datetime(df['DOTTING_TIME'], unit='ms')

    # 转成小时
    df['DOTTING_TIME'] = t.dt.hour + t.dt.minute / 60

类别特征取值关于Label箱型图可视化：

# 因为 name和 regionCode的类别太稀疏了，这里我们把不稀疏的几类画一下
categorical_features =['model','brand','bodyType','fuelType','gearbox','notRepairedDamage']
for c in categorical_features:
    Train_data[c] = Train_data[c].astype('category')
    if Train_data[c].isnull().any():
        Train_data[c] = Train_data[c].cat.add_categories(['MISSING'])
        Train_data[c] = Train_data[c].fillna('MISSING')

def boxplot(x, y, **kwargs):
    sns.boxplot(x=x, y=y)
    x=plt.xticks(rotation=90)

f = pd.melt(Train_data, id_vars=['price'], value_vars=categorical_features)
g = sns.FacetGrid(f, col="variable",  col_wrap=2, sharex=False, sharey=False, size=5)
g = g.map(boxplot, "value", "price")

查看训练集测试集中特征属性只有一值的特征：

one_value_fea = [col for col in data_train.columns if data_train[col].nunique() <= 1]
one_value_fea_test = [col for col in data_test_a.columns if data_test_a[col].nunique() <= 1]

---------------------数据预处理与特征工程--------------------

• 异常处理：
  • 通过箱线图（或 3-Sigma）分析删除异常值；
  • BOX-COX 转换（处理有偏分布）；
  • 长尾截断；
• 特征归一化/标准化：
  • 标准化（转换为标准正态分布）；
  • 归一化（抓换到 [0,1] 区间）；
  • 针对幂律分布，可以采用公式： log(1+x1+median)log(1+x1+median)
• 数据分桶：
  • 等频分桶；
  • 等距分桶；
  • Best-KS 分桶（类似利用基尼指数进行二分类）；
  • 卡方分桶；
• 缺失值处理：
  • 不处理（针对类似 XGBoost 等树模型）；
  • 删除（缺失数据太多）；
  • 插值补全，包括均值/中位数/众数/建模预测/多重插补/压缩感知补全/矩阵补全等；
  • 分箱，缺失值一个箱；
• 特征构造：
  • 构造统计量特征，报告计数、求和、比例、标准差等；
  • 时间特征，包括相对时间和绝对时间，节假日，双休日等；
  • 地理信息，包括分箱，分布编码等方法；
  • 非线性变换，包括 log/ 平方/ 根号等；
  • 特征组合，特征交叉；
  • 仁者见仁，智者见智。
• 特征筛选
  • 过滤式（filter）：先对数据进行特征选择，然后在训练学习器，常见的方法有 Relief/方差选择发/相关系数法/卡方检验法/互信息法；
  • 包裹式（wrapper）：直接把最终将要使用的学习器的性能作为特征子集的评价准则，常见方法有 LVM（Las Vegas Wrapper） ；
  • 嵌入式（embedding）：结合过滤式和包裹式，学习器训练过程中自动进行了特征选择，常见的有 lasso 回归；
• 降维
  • PCA/ LDA/ ICA；
  • 特征选择也是一种降维。

采用箱型图原理去除异常值代码：

# 这里我包装了一个异常值处理的代码，可以随便调用。
def outliers_proc(data, col_name, scale=3):
    """
    用于清洗异常值，默认用 box_plot（scale=3）进行清洗
    :param data: 接收 pandas 数据格式
    :param col_name: pandas 列名
    :param scale: 尺度
    :return:
    """

    def box_plot_outliers(data_ser, box_scale):
        """
        利用箱线图去除异常值
        :param data_ser: 接收 pandas.Series 数据格式
        :param box_scale: 箱线图尺度，
        :return:
        """
        iqr = box_scale * (data_ser.quantile(0.75) - data_ser.quantile(0.25))
        val_low = data_ser.quantile(0.25) - iqr
        val_up = data_ser.quantile(0.75) + iqr
        rule_low = (data_ser < val_low)
        rule_up = (data_ser > val_up)
        return (rule_low, rule_up), (val_low, val_up)

    data_n = data.copy()
    data_series = data_n[col_name]
    rule, value = box_plot_outliers(data_series, box_scale=scale)
    index = np.arange(data_series.shape[0])[rule[0] | rule[1]]
    print("Delete number is: {}".format(len(index)))
    data_n = data_n.drop(index)
    data_n.reset_index(drop=True, inplace=True)
    print("Now column number is: {}".format(data_n.shape[0]))
    index_low = np.arange(data_series.shape[0])[rule[0]]
    outliers = data_series.iloc[index_low]
    print("Description of data less than the lower bound is:")
    print(pd.Series(outliers).describe())
    index_up = np.arange(data_series.shape[0])[rule[1]]
    outliers = data_series.iloc[index_up]
    print("Description of data larger than the upper bound is:")
    print(pd.Series(outliers).describe())
    
    fig, ax = plt.subplots(1, 2, figsize=(10, 7))
    sns.boxplot(y=data[col_name], data=data, palette="Set1", ax=ax[0])
    sns.boxplot(y=data_n[col_name], data=data_n, palette="Set1", ax=ax[1])
    return data_n

合并训练集和测试集：

"""
通常会将DataFrame类型的训练集和测试集合并起来，方便后面做特征工程
注意：最好记录下前多少行是训练集，方便后面做特征工程
"""
data_df = pd.concat([train_df, test_df], axis=0, ignore_index=True)

线性填充缺失值：

#从scipy库中导入插值需要的方法 interpolate
from scipy import interpolate
#开始使用interpolate进行线性填充
pre_train_df = pre_train_df.interpolate(kind='linear',axis = 0)
test_df = test_df.interpolate(kind='linear',axis = 0)

中位数、众数、均值填补缺失值：

#中位数
df.fillna(df.median(),inplace=True)
#众数
df.fillna(df.mode(),inplace=True)
#平均值
df.fillna(method=a.mean(),inplace=True)

向上向下填充缺失值（可能会造成数据泄露）：

#向上填补缺失值
pre_train_df = pre_train_df.fillna(method="bfill")
#向下填补缺失值
test_df = test_df.fillna(method="ffill")

对str类别特征做LabelEncoder：

cate_feat = ['QUEUE_ID','CU','QUEUE_TYPE']

for i in cate_feat:
    lbl = LabelEncoder() 
    train_df[i] = lbl.fit_transform(train_df[i].astype(str))
    test_df[i] = lbl.fit_transform(test_df[i].astype(str))

一阶差分(常解决时间序列预测问题)：

pre_train_df["outdoorTemp"].shift(3)

 shift()方法简单来说：将该列数据整体向后退3个位置,看下面的输出自己理解，比较简单。 

#同期值
for f2 in tqdm(['outdoorTemp', 'outdoorHum', 'outdoorAtmo', 'indoorHum', 'indoorAtmo']):
    #这样计算同期值好像不太准确
    pre_train_df["ago_1hour_"+f2] = pre_train_df[f2].shift(1*60)
for f2 in tqdm(['outdoorTemp', 'outdoorHum', 'outdoorAtmo', 'indoorHum', 'indoorAtmo']):
    test_df["ago_1hour_"+f2] = test_df[f2].shift(1*2)

#一阶差分
for f3 in tqdm(['outdoorTemp', 'outdoorHum', 'outdoorAtmo', 'indoorHum', 'indoorAtmo']):
    #同期的差距
    pre_train_df["ago_1hour_"+f3+"_trend"] = pre_train_df[f3] - pre_train_df["ago_1hour_"+f3]
for f3 in tqdm(['outdoorTemp', 'outdoorHum', 'outdoorAtmo', 'indoorHum', 'indoorAtmo']):
    test_df["ago_1hour_"+f3+"_trend"] = test_df[f3] - test_df["ago_1hour_"+f3]

n阶差分开窗(常解决时间序列预测问题)：

#统计过去n阶差分开窗-----也称为划窗统计
def get_statis_n_days_num(data,col,n):
    temp = pd.DataFrame()
    for i in range(n):
        #构建移动窗口
        temp = pd.concat([temp,data[col].shift((i+1)*24)],axis = 1)
    #统计窗口中数据的平均值
    data["avg_"+str(n)+"_days_"+col] = temp.mean(axis = 1)
    #中位值
    data["median_"+str(n)+"_days_"+col] = temp.median(axis = 1)
    #最大值
    data["max_"+str(n)+"_days_"+col] = temp.max(axis = 1)
    #最小值
    data["min_"+str(n)+"_days_"+col] = temp.min(axis = 1)
    #方差
    data["std_"+str(n)+"_days_"+col] = temp.std(axis = 1)
    #平均绝对离查
    data["mad_"+str(n)+"_days_"+col] = temp.mad(axis = 1)
    #偏度
    data["skew_"+str(n)+"_days_"+col] = temp.skew(axis = 1)
    #峰度
    data["kurt_"+str(n)+"_days_"+col] = temp.kurt(axis = 1)
    #1分数
    data["q1_"+str(n)+"_days_"+col] = temp.quantile(q = 0.25,axis = 1)
    #3分数
    data["q3_"+str(n)+"_days_"+col] = temp.quantile(q = 0.75,axis = 1)
    #协方差
    data["var_"+str(n)+"_days_"+col] = data["std_"+str(n) + "_days_" + col]/data["avg_"+str(n)+"_days_"+col] #离散系数
    return data

#然后对气压、压强、温度等构造统计特征
for tongji_f in tqdm(['outdoorTemp', 'outdoorHum', 'outdoorAtmo', 'indoorHum', 'indoorAtmo']):
    get_statis_n_days_num(pre_train_df,tongji_f,5)
    get_statis_n_days_num(test_df,tongji_f,5)

基本聚合特征：

group_feats = []
for f in tqdm(['outdoorTemp', 'outdoorHum', 'outdoorAtmo', 'indoorHum', 'indoorAtmo']):
    #对每个小时的数据进行聚合
    data_df['MDH_{}_medi'.format(f)] = data_df.groupby(['month', 'day', 'hour'])[f].transform('median')
    data_df['MDH_{}_mean'.format(f)] = data_df.groupby(['month', 'day', 'hour'])[f].transform('mean')
    data_df['MDH_{}_max'.format(f)] = data_df.groupby(['month', 'day', 'hour'])[f].transform('max')
    data_df['MDH_{}_min'.format(f)] = data_df.groupby(['month', 'day', 'hour'])[f].transform('min')
    data_df['MDH_{}_std'.format(f)] = data_df.groupby(['month', 'day', 'hour'])[f].transform('std')

    group_feats.append('MDH_{}_medi'.format(f))
    group_feats.append('MDH_{}_mean'.format(f))
#统计一天室外最小温度、最大温度、平均温度
data_df["eveyday_max_temp"] = data_df.groupby(["month","day"])["outdoorTemp"].transform("max")
data_df["eveyday_min_temp"] = data_df.groupby(["month","day"])["outdoorTemp"].transform("min")
data_df["eveyday_mean_temp"] = data_df.groupby(["month","day"])["outdoorTemp"].transform("mean")

基本交叉特征：

# 基本交叉特征（包含基础特征、基础特征的统计特征【平均值，中位数】）---常用方法是加减乘除，这里直接是除
for f1 in tqdm(['outdoorTemp', 'outdoorHum', 'outdoorAtmo', 'indoorHum', 'indoorAtmo'] + group_feats):   
    for f2 in ['outdoorTemp', 'outdoorHum', 'outdoorAtmo', 'indoorHum', 'indoorAtmo'] + group_feats:
        if f1 != f2:
            colname = '{}_{}_ratio'.format(f1, f2)
            data_df[colname] = data_df[f1].values / data_df[f2].values

将时间属性转换成连续值：

#将时间具体到哪一天，14,15,16，........,45    
data_df['dt'] = data_df['day'].values + (data_df['month'].values - 3) * 31

提取历史信息特征：

# 历史信息提取
data_df['dt'] = data_df['day'].values + (data_df['month'].values - 3) * 31

for f in ['outdoorTemp', 'outdoorHum', 'outdoorAtmo', 'indoorHum', 'indoorAtmo', 'temperature']:
    tmp_df = pd.DataFrame()
    for t in tqdm(range(15, 45)):
        #当天--->三天前这3天的统计特征
        #agg与前面transform不同的是，groupby的组数是多少，agg的结果的个数就是多少
        tmp = data_df[(data_df['dt'] < t)| (data_df['dt'] > t-3)].groupby(['hour'])[f].agg({'mean'}).reset_index()
        tmp.columns = ['hour', 'hit_{}_mean'.format(f)]
        tmp['dt'] = t
        #tem ---> ['dt' 'hour' 'hit_{}_mean']
        tmp_df = tmp_df.append(tmp)

    data_df = data_df.merge(tmp_df, on=['dt', 'hour'], how='left')

data_df = data_df.fillna(method='bfill')

数值特征离散化并构造统计特征：

# 离散化---为了让训练出来的模型更具鲁棒性
for f in ['outdoorTemp', 'outdoorHum', 'outdoorAtmo', 'indoorHum', 'indoorAtmo']:
    data_df[f + '_20_bin'] = pd.cut(data_df[f], 20, duplicates='drop').apply(lambda x: x.left).astype(int)
    data_df[f + '_50_bin'] = pd.cut(data_df[f], 50, duplicates='drop').apply(lambda x: x.left).astype(int)
    data_df[f + '_100_bin'] = pd.cut(data_df[f], 100, duplicates='drop').apply(lambda x: x.left).astype(int)
    data_df[f + '_200_bin'] = pd.cut(data_df[f], 200, duplicates='drop').apply(lambda x: x.left).astype(int)

统计每个bins中数据的中位数，均值，最大值，最小值：

for f1 in tqdm(
        ['outdoorTemp_20_bin', 'outdoorHum_20_bin', 'outdoorAtmo_20_bin', 'indoorHum_20_bin', 'indoorAtmo_20_bin']):
    for f2 in ['outdoorTemp', 'outdoorHum', 'outdoorAtmo', 'indoorHum', 'indoorAtmo']:
        data_df['{}_{}_medi'.format(f1, f2)] = data_df.groupby([f1])[f2].transform('median')
        data_df['{}_{}_mean'.format(f1, f2)] = data_df.groupby([f1])[f2].transform('mean')
        data_df['{}_{}_max'.format(f1, f2)] = data_df.groupby([f1])[f2].transform('max')
        data_df['{}_{}_min'.format(f1, f2)] = data_df.groupby([f1])[f2].transform('min')

for f1 in tqdm(
        ['outdoorTemp_50_bin', 'outdoorHum_50_bin', 'outdoorAtmo_50_bin', 'indoorHum_50_bin', 'indoorAtmo_50_bin']):
    for f2 in ['outdoorTemp', 'outdoorHum', 'outdoorAtmo', 'indoorHum', 'indoorAtmo']:
        data_df['{}_{}_medi'.format(f1, f2)] = data_df.groupby([f1])[f2].transform('median')
        data_df['{}_{}_mean'.format(f1, f2)] = data_df.groupby([f1])[f2].transform('mean')
        data_df['{}_{}_max'.format(f1, f2)] = data_df.groupby([f1])[f2].transform('max')
        data_df['{}_{}_min'.format(f1, f2)] = data_df.groupby([f1])[f2].transform('min')

for f1 in tqdm(['outdoorTemp_100_bin', 'outdoorHum_100_bin', 'outdoorAtmo_100_bin', 'indoorHum_100_bin',
                'indoorAtmo_100_bin']):
    for f2 in ['outdoorTemp', 'outdoorHum', 'outdoorAtmo', 'indoorHum', 'indoorAtmo']:
        data_df['{}_{}_medi'.format(f1, f2)] = data_df.groupby([f1])[f2].transform('median')
        data_df['{}_{}_mean'.format(f1, f2)] = data_df.groupby([f1])[f2].transform('mean')
        data_df['{}_{}_max'.format(f1, f2)] = data_df.groupby([f1])[f2].transform('max')
        data_df['{}_{}_min'.format(f1, f2)] = data_df.groupby([f1])[f2].transform('min')

for f1 in tqdm(['outdoorTemp_200_bin', 'outdoorHum_200_bin', 'outdoorAtmo_200_bin', 'indoorHum_200_bin',
                'indoorAtmo_200_bin']):
    for f2 in ['outdoorTemp', 'outdoorHum', 'outdoorAtmo', 'indoorHum', 'indoorAtmo']:
        data_df['{}_{}_medi'.format(f1, f2)] = data_df.groupby([f1])[f2].transform('median')
        data_df['{}_{}_mean'.format(f1, f2)] = data_df.groupby([f1])[f2].transform('mean')
        data_df['{}_{}_max'.format(f1, f2)] = data_df.groupby([f1])[f2].transform('max')
        data_df['{}_{}_min'.format(f1, f2)] = data_df.groupby([f1])[f2].transform('min')

 

处理时间特征转换为距某天的时间：

# 使用时间：data['creatDate'] - data['regDate']，反应汽车使用时间，一般来说价格与使用时间成反比
# 不过要注意，数据里有时间出错的格式，所以我们需要 errors='coerce'
data['used_time'] = (pd.to_datetime(data['creatDate'], format='%Y%m%d', errors='coerce') - 
                            pd.to_datetime(data['regDate'], format='%Y%m%d', errors='coerce')).dt.days

apply的灵活应用：

# 从邮编中提取城市信息，相当于加入了先验知识
data['city'] = data['regionCode'].apply(lambda x : str(x)[:-3])
#别忘了设置数据类型，许多模型默认不处理object类型的数据
data['city'] = data['city'].astype("int32")

提取某个类别特征关于Label的统计量（均值、方差.....）：

# 计算某品牌的销售统计量，同学们还可以计算其他特征的统计量
# 这里要以 train 的数据计算统计量
Train_gb = Train_data.groupby("brand")
all_info = {}
for kind, kind_data in Train_gb:
    info = {}
    kind_data = kind_data[kind_data['price'] > 0]
    info['brand_amount'] = len(kind_data)
    info['brand_price_max'] = kind_data.price.max()
    info['brand_price_median'] = kind_data.price.median()
    info['brand_price_min'] = kind_data.price.min()
    info['brand_price_sum'] = kind_data.price.sum()
    info['brand_price_std'] = kind_data.price.std()
    info['brand_price_average'] = round(kind_data.price.sum() / (len(kind_data) + 1), 2)
    all_info[kind] = info
brand_fe = pd.DataFrame(all_info).T.reset_index().rename(columns={"index": "brand"})
data = data.merge(brand_fe, how='left', on='brand')

数据分桶：

# 数据分桶 以 power 为例
# 这时候我们的缺失值也进桶了，
# 为什么要做数据分桶呢，原因有很多，= =
# 1. 离散后稀疏向量内积乘法运算速度更快，计算结果也方便存储，容易扩展；
# 2. 离散后的特征对异常值更具鲁棒性，如 age>30 为 1 否则为 0，对于年龄为 200 的也不会对模型造成很大的干扰；
# 3. LR 属于广义线性模型，表达能力有限，经过离散化后，每个变量有单独的权重，这相当于引入了非线性，能够提升模型的表达能力，加大拟合；
# 4. 离散后特征可以进行特征交叉，提升表达能力，由 M+N 个变量编程 M*N 个变量，进一步引入非线形，提升了表达能力；
# 5. 特征离散后模型更稳定，如用户年龄区间，不会因为用户年龄长了一岁就变化

# 当然还有很多原因，LightGBM 在改进 XGBoost 时就增加了数据分桶，增强了模型的泛化性

bin = [i*10 for i in range(31)]
data['power_bin'] = pd.cut(data['power'], bin, labels=False)
data[['power_bin', 'power']].head()

删除无用的属性列：

# 删除不需要的数据
data = data.drop(['creatDate', 'regDate', 'regionCode'], axis=1)

Log处理：

# 我们对其取 log，在做归一化
from sklearn import preprocessing
min_max_scaler = preprocessing.MinMaxScaler()
data['power'] = np.log(data['power'] + 1) 
data['power'].plot.hist()

归一化：

def max_min(x):
    return (x - np.min(x)) / (np.max(x) - np.min(x))

data['brand_amount'] = ((data['brand_amount'] - np.min(data['brand_amount'])) / 
                        (np.max(data['brand_amount']) - np.min(data['brand_amount'])))
data['brand_price_average'] = ((data['brand_price_average'] - np.min(data['brand_price_average'])) / 
                               (np.max(data['brand_price_average']) - np.min(data['brand_price_average'])))
data['brand_price_max'] = ((data['brand_price_max'] - np.min(data['brand_price_max'])) / 
                           (np.max(data['brand_price_max']) - np.min(data['brand_price_max'])))
data['brand_price_median'] = ((data['brand_price_median'] - np.min(data['brand_price_median'])) /
                              (np.max(data['brand_price_median']) - np.min(data['brand_price_median'])))
data['brand_price_min'] = ((data['brand_price_min'] - np.min(data['brand_price_min'])) / 
                           (np.max(data['brand_price_min']) - np.min(data['brand_price_min'])))
data['brand_price_std'] = ((data['brand_price_std'] - np.min(data['brand_price_std'])) / 
                           (np.max(data['brand_price_std']) - np.min(data['brand_price_std'])))
data['brand_price_sum'] = ((data['brand_price_sum'] - np.min(data['brand_price_sum'])) / 
                           (np.max(data['brand_price_sum']) - np.min(data['brand_price_sum'])))

对类别特征进行one_hat编码：

# 对类别特征进行 OneEncoder
data = pd.get_dummies(data, columns=['model', 'brand', 'bodyType', 'fuelType',
                                     'gearbox', 'notRepairedDamage', 'power_bin'])

相关系数过滤式特征筛选：

# 相关性分析
print(data['power'].corr(data['price'], method='spearman'))
print(data['kilometer'].corr(data['price'], method='spearman'))
print(data['brand_amount'].corr(data['price'], method='spearman'))
print(data['brand_price_average'].corr(data['price'], method='spearman'))
print(data['brand_price_max'].corr(data['price'], method='spearman'))
print(data['brand_price_median'].corr(data['price'], method='spearman'))

构建组合特征：

## 组合特征
data = pd.concat([train, test], axis=0, ignore_index=True)
for col1 in [['QUEUE_ID'],['CU'],['QUEUE_TYPE'],['QUEUE_ID','CU'],['QUEUE_ID','QUEUE_TYPE'],['CU','QUEUE_TYPE']]:
    name = '_'.join(col1)
    for col2 in ['CPU_USAGE','MEM_USAGE','DISK_USAGE']:
        tmp = data.groupby(col1)[col2].agg({'mean','median','std'}).reset_index()
        tmp.columns = col1 + [name+'_'+col2+'_'+f for f in tmp.columns if f not in col1]
        train_df = train_df.merge(tmp, on=col1, how='left')
        test_df = test_df.merge(tmp, on=col1, how='left')

构建衰减求和特征：

## 衰减求和特征
#power(x, y) 函数，计算 x 的 y 次方
#np.arange(5)----array([0, 1, 2, 3, 4])
pw09 = np.power(0.9, np.arange(5))
pw08 = np.power(0.8, np.arange(5))
pw07 = np.power(0.7, np.arange(5))
pw06 = np.power(0.6, np.arange(5))
pw05 = np.power(0.5, np.arange(5))
for col in ['CPU_USAGE','MEM_USAGE','DISK_USAGE']:
    for i, pw in enumerate([pw09,pw08,pw07,pw06,pw05]):
        colname = col+'_pw_'+str(9-i)
        train_df[colname] = train_df[col]*pw[0] + \
                            train_df[col+'_shift1']*pw[1] + \
                            train_df[col+'_shift2']*pw[2] + \
                            train_df[col+'_shift3']*pw[3] + \
                            train_df[col+'_shift4']*pw[4]
        
        test_df[colname] = test_df[col]*pw[0] + \
                           test_df[col+'_shift1']*pw[1] + \
                           test_df[col+'_shift2']*pw[2] + \
                           test_df[col+'_shift3']*pw[3] + \
                           test_df[col+'_shift4']*pw[4]

构建历史平移特征：

方差过滤式特征筛选：

这是通过特征本身的方差来筛选特征的类。比如一个特征本身的方差很小，就表示样本在这个特征上基本没有差 异，可能特征中的大多数值都一样，甚至整个特征的取值都相同，那这个特征对于样本区分没有什么作用。所以无 论接下来的特征工程要做什么，都要优先消除方差为0的特征。VarianceThreshold有重要参数threshold，表示方 差的阈值，表示舍弃所有方差小于threshold的特征，不填默认为0，即删除所有的记录都相同的特征。

from sklearn.feature_selection import VarianceThreshold
#实例化，不填参数默认方差为0
selector = VarianceThreshold() 
#获取删除不合格特征之后的新特征矩阵
X_var0 = selector.fit_transform(X) 
#也可以直接写成 X = VairanceThreshold().fit_transform(X)
X_var0.shape

当特征是二分类时，特征的取值就是伯努利随机变量，这些变量的方差可以计算为：

其中X是特征矩阵，p是二分类特征中的一类在这个特征中所占的概率。

#若特征是伯努利随机变量，假设p=0.8，即二分类特征中某种分类占到80%以上的时候删除特征
X_bvar = VarianceThreshold(.8 * (1 - .8)).fit_transform(X)
X_bvar.shape

卡方过滤特征筛选(离散型Label):

卡方过滤是专门针对离散型标签（即分类问题）的相关性过滤。卡方检验类feature_selection.chi2计算每个非负特征和标签之间的卡方统计量，并依照卡方统计量由高到低为特征排名。再结合feature_selection.SelectKBest 这个可以输入”评分标准“来选出前K个分数最高的特征的类，我们可以借此除去最可能独立于标签，与我们分类目 的无关的特征。 另外，如果卡方检验检测到某个特征中所有的值都相同，会提示我们使用方差先进行方差过滤。并且，刚才我们已 经验证过，当我们使用方差过滤筛选掉一半的特征后，模型的表现时提升的。因此在这里，我们使用threshold=中 位数时完成的方差过滤的数据来做卡方检验（如果方差过滤后模型的表现反而降低了，那我们就不会使用方差过滤 后的数据，而是使用原数据）：

from sklearn.ensemble import RandomForestClassifier as RFC
from sklearn.model_selection import cross_val_score
from sklearn.feature_selection import SelectKBest
from sklearn.feature_selection import chi2
#假设在这里我一直我需要300个特征
X_fschi = SelectKBest(chi2, k=300).fit_transform(X_fsvar, y)
X_fschi.shape

那如何设置一个最佳的K值呢？在现实数据中，数据量很大，模型很复杂的时候，我们也许不能先去跑一遍模型看 看效果，而是希望最开始就能够选择一个最优的超参数k。那第一个方法，就是我们之前提过的学习曲线：

#======【TIME WARNING: 5 mins】======#
%matplotlib inline
import matplotlib.pyplot as plt
score = []
for i in range(390,200,-10):
    X_fschi = SelectKBest(chi2, k=i).fit_transform(X_fsvar, y)
    once = cross_val_score(RFC(n_estimators=10,random_state=0),X_fschi,y,cv=5).mean()
    score.append(once)
plt.plot(range(350,200,-10),score)
plt.show()

卡方检验的本质是推测两组数据之间的差异，其检验的原假设是”两组数据是相互独立的”。卡方检验返回卡方值和 P值两个统计量，其中卡方值很难界定有效的范围，而p值，我们一般使用0.01或0.05作为显著性水平，即p值判断 的边界，具体我们可以这样来看:

 

P值	<=0.05或0.01	>0.05或0.01
数据差异	差异不是自然形成的	这些差异是很自然的样本误差
相关性	两组数据是相关的	两组数据是相互独立的
原假设	拒绝原假设，接受备择假设	接受原假设

 

从特征工程的角度，我们希望选取卡方值很大，p值小于0.05的特征，即和标签是相关联的特征。而调用 SelectKBest之前，我们可以直接从chi2实例化后的模型中获得各个特征所对应的卡方值和P值。

chivalue, pvalues_chi = chi2(X_fsvar,y)
chivalue
pvalues_chi
#k取多少？我们想要消除所有p值大于设定值，比如0.05或0.01的特征：
k = chivalue.shape[0] - (pvalues_chi > 0.05).sum()
#X_fschi = SelectKBest(chi2, k=填写具体的k).fit_transform(X_fsvar, y)
#cross_val_score(RFC(n_estimators=10,random_state=0),X_fschi,y,cv=5).mean()

F检验特征筛选:

F检验，又称ANOVA，方差齐性检验，是用来捕捉每个特征与标签之间的线性关系的过滤方法。它即可以做回归也可以做分类，因此包含feature_selection.f_classif（F检验分类）和feature_selection.f_regression（F检验回归）两个类。其中F检验分类用于标签是离散型变量的数据，而F检验回归用于标签是连续型变量的数据。 和卡方检验一样，这两个类需要和类SelectKBest连用，并且我们也可以直接通过输出的统计量来判断我们到底要 设置一个什么样的K。需要注意的是，F检验在数据服从正态分布时效果会非常稳定，因此如果使用F检验过滤，我 们会先将数据转换成服从正态分布的方式。 F检验的本质是寻找两组数据之间的线性关系，其原假设是”数据不存在显著的线性关系“。它返回F值和p值两个统 计量。和卡方过滤一样，我们希望选取p值小于0.05或0.01的特征，这些特征与标签时显著线性相关的，而p值大于 0.05或0.01的特征则被我们认为是和标签没有显著线性关系的特征，应该被删除。以F检验的分类为例，我们继续 在数字数据集上来进行特征选择：

from sklearn.feature_selection import f_classif
F, pvalues_f = f_classif(X_fsvar,y)
F
pvalues_f
k = F.shape[0] - (pvalues_f > 0.05).sum()
#X_fsF = SelectKBest(f_classif, k=填写具体的k).fit_transform(X_fsvar, y)
#cross_val_score(RFC(n_estimators=10,random_state=0),X_fsF,y,cv=5).mean()

互信息法特征筛选:

互信息法是用来捕捉每个特征与标签之间的任意关系（包括线性和非线性关系）的过滤方法。和F检验相似，它既可以做回归也可以做分类，并且包含两个类feature_selection.mutual_info_classif（互信息分类）和 feature_selection.mutual_info_regression（互信息回归）。这两个类的用法和参数都和F检验一模一样，不过 互信息法比F检验更加强大，F检验只能够找出线性关系，而互信息法可以找出任意关系。 互信息法不返回p值或F值类似的统计量，它返回“每个特征与目标之间的互信息量的估计”，这个估计量在[0,1]之间 取值，为0则表示两个变量独立，为1则表示两个变量完全相关。以互信息分类为例的代码如下：

from sklearn.feature_selection import mutual_info_classif as MIC
result = MIC(X_fsvar,y)
k = result.shape[0] - sum(result <= 0)
#X_fsmic = SelectKBest(MIC, k=填写具体的k).fit_transform(X_fsvar, y)
#cross_val_score(RFC(n_estimators=10,random_state=0),X_fsmic,y,cv=5).mean()

包裹式特征筛选：

!pip install mlxtend  # 不要点，下载速度很慢

# k_feature 太大会很难跑，没服务器，所以提前 interrupt 了
from mlxtend.feature_selection import SequentialFeatureSelector as SFS
from sklearn.linear_model import LinearRegression
sfs = SFS(LinearRegression(),
           k_features=10,
           forward=True,
           floating=False,
           scoring = 'r2',
           cv = 0)
x = data.drop(['price'], axis=1)
x = x.fillna(0)
y = data['price']
sfs.fit(x, y)
""
将筛选的特征打印出来
""
sfs.k_feature_names_ 

# 画出来，可以看到边际效益
from mlxtend.plotting import plot_sequential_feature_selection as plot_sfs
import matplotlib.pyplot as plt
fig1 = plot_sfs(sfs.get_metric_dict(), kind='std_dev')
plt.grid()
plt.show()

嵌入式特征选择：

 

总结：

• 特征工程是比赛中最至关重要的的一块，特别的传统的比赛，大家的模型可能都差不多，调参带来的效果增幅是非常有限的，但特征工程的好坏往往会决定了最终的排名和成绩。
• 特征工程的主要目的还是在于将数据转换为能更好地表示潜在问题的特征，从而提高机器学习的性能。比如，异常值处理是为了去除噪声，填补缺失值可以加入先验知识等。
• 特征构造也属于特征工程的一部分，其目的是为了增强数据的表达。
• 有些比赛的特征是匿名特征，这导致我们并不清楚特征相互直接的关联性，这时我们就只有单纯基于特征进行处理，比如装箱，groupby，agg 等这样一些操作进行一些特征统计，此外还可以对特征进行进一步的 log，exp 等变换，或者对多个特征进行四则运算（如上面我们算出的使用时长），多项式组合等然后进行筛选。由于特性的匿名性其实限制了很多对于特征的处理，当然有些时候用 NN 去提取一些特征也会达到意想不到的良好效果。
• 对于知道特征含义（非匿名）的特征工程，特别是在工业类型比赛中，会基于信号处理，频域提取，丰度，偏度等构建更为有实际意义的特征，这就是结合背景的特征构建，在推荐系统中也是这样的，各种类型点击率统计，各时段统计，加用户属性的统计等等，这样一种特征构建往往要深入分析背后的业务逻辑或者说物理原理，从而才能更好的找到 magic。
• 当然特征工程其实是和模型结合在一起的，这就是为什么要为 LR NN 做分桶和特征归一化的原因，而对于特征的处理效果和特征重要性等往往要通过模型来验证。
• 总的来说，特征工程是一个入门简单，但想精通非常难的一件事。

 

----------------------------建模调参-------------------------------

lgb使用模板：

""
trn_x、val_x、train_x：为我们切分好DataFrame类型的数据集，不含有Lable那一列
trn_y、val_y、train_y：我们对label进行分片得到的series
""
train_matrix = lgb.Dataset(trn_x, label=trn_y)
valid_matrix = lgb.Dataset(val_x, label=val_y)
data_matrix = lgb.Dataset(train_x, label=train_y)

params = {
        #弱分类器使用gbdt
        'boosting_type': 'gbdt',
        #损失函数使用均方误差
        'objective': 'mse',
        #为了防止过拟合，该节点的信息增益如果小于该值就不在往下分裂
        'min_child_weight': 5,
        #为了防止过拟合，设定的没棵树叶子节点数目
        'num_leaves': 2 ** 8,
        ""
        feature_fraction 小于 1.0, LightGBM 将会在每次迭代中随机选择部分特征. 例如, 如果设置 
        为 0.8, 将会在每棵树训练之前选择 80% 的特征

        ""         
        'feature_fraction': 0.5,
        #但是它将在不进行重采样的情况下随机选择部分数据
        'bagging_fraction': 0.5,
        #每几次迭代执行bagging
        'bagging_freq': 1,
        #学习率
        'learning_rate': 0.001,
        #种子，这个在比赛的时候可能会改
        'seed': 2020
     }

#kfold = KFold(n_splits=5, random_state=42)
#for fold_id, (trn_idx, val_idx) in enumerate(kfold.split(df_train, df_train[target])):
"""
使用fit进行训练，不用封装DMatrix(train_matrix = lgb.Dataset(trn_x, label=trn_y))
"""
#训练5000次，1000次早停
model = lgb.train(params, train_matrix, 50000, valid_sets=[train_matrix, valid_matrix], 
        verbose_eval=500,
        early_stopping_rounds=1000)

"""
训练过程打印重要性属性
"""
print(pd.DataFrame({
        'column': feature_names,
        'importance': lgb_model.feature_importances_,
        }).sort_values(by='importance',ascending = False))

model2 = lgb.train(params, data_matrix, model.best_iteration)
#模型训练完成，开始预测label
val_pred = model.predict(val_x, num_iteration=model2.best_iteration).reshape(-1, 1)

test_pred = model.predict(test_x, num_iteration=model2.best_iteration).reshape(-1, 1)

xgb使用模板：

#使用Xgboost自带的读取格式DMatrix(),并忽略np.nan数据
train_matrix = xgb.DMatrix(trn_x, label=trn_y, missing=np.nan)
valid_matrix = xgb.DMatrix(val_x, label=val_y, missing=np.nan)
test_matrix = xgb.DMatrix(test_x, missing=np.nan)
#mae 均方根误差
params = {'booster': 'gbtree',
          'eval_metric': 'mae',
          'min_child_weight': 5,
          'max_depth': 8,
          'subsample': 0.5,
          'colsample_bytree': 0.5,
          'eta': 0.001,
          'seed': 2020,
          'nthread': 36,
          'silent': True,
        }

watchlist = [(train_matrix, 'train'), (valid_matrix, 'eval')]

model = xgb.train(params, train_matrix, num_boost_round=50000, evals=watchlist, 
              verbose_eval=500,
              early_stopping_rounds=1000)

val_pred = model.predict(valid_matrix, ntree_limit=model.best_ntree_limit).reshape(-1, 1)

test_pred = model.predict(test_matrix, ntree_limit=model.best_ntree_limit).reshape(-1, 1)

cat使用模板：

params = {'learning_rate': 0.001, 
          'depth': 5, 
          'l2_leaf_reg': 10, 
          'bootstrap_type': 'Bernoulli',
          'od_type': 'Iter', 
          'od_wait': 50,
          'random_seed': 11, 
          'allow_writing_files': False}

model = cat(iterations=20000, **params)

model.fit(trn_x, trn_y, eval_set=(val_x, val_y),cat_features=[],use_best_model=True,verbose=500)

val_pred = model.predict(val_x)
test_pred = model.predict(test_x)

SGDRegressor（随机梯度下降回归）使用模板：

params = {
            'loss': 'squared_loss',
            'penalty': 'l2',
            'alpha': 0.00001,
            'random_state': 2020,
        }

model = SGDRegressor(**params)
model.fit(trn_x, trn_y)
val_pred = model.predict(val_x)
test_pred = model.predict(test_x)

ridge(岭回归)使用模板：

params = {
         'alpha': 1.0,
         'random_state': 2020,
       }

model = Ridge(**params)
model.fit(trn_x, trn_y)
val_pred = model.predict(val_x)
test_pred = model.predict(test_x)

LSTM-CNN使用模板：

下列代码取自华为举办的“大数据时代的Serverless工作负载预测”的baseline.

# coding:utf-8
import os
import numpy as np
import pandas as pd
from tensorflow.keras import Input, Model
from tensorflow.keras.initializers import he_normal
from tensorflow.keras.layers import LSTM, Dense, Conv1D, LeakyReLU, Concatenate
from tensorflow.keras.optimizers import Adam
from sklearn.preprocessing import MinMaxScaler

train = pd.read_csv('E:/Datas/train.csv')
train = train[train.STATUS == 'available']
train = train[train.PLATFORM == 'x86_64']
train = train[train.RESOURCE_TYPE == 'vm']
train = train.reset_index(drop=True)

"""
将数据按照 QUEUE_ID 和 时间进行排序
"""
train = train.sort_values(by=['QUEUE_ID', 'DOTTING_TIME']).reset_index(drop=True)
test = pd.read_csv('E:/Datas/evaluation_public.csv')
test = test.sort_values(by=['ID', 'DOTTING_TIME']).reset_index(drop=True)

"""
构造一些特征转化和简单的差分特征
"""
train['used_cpu'] = train['CU'] * train['CPU_USAGE'] / 100
train['used_mem'] = train['CU'] * 4 * train['MEM_USAGE'] / 100
test['used_cpu'] = test['CU'] * test['CPU_USAGE'] / 100
test['used_mem'] = test['CU'] * 4 * test['MEM_USAGE'] / 100
train['to_run_jobs'] = train['LAUNCHING_JOB_NUMS'] - train['RUNNING_JOB_NUMS']
test['to_run_jobs'] = test['LAUNCHING_JOB_NUMS'] - test['RUNNING_JOB_NUMS']
# 差分
train['used_cpu_diff1'] = train.groupby(['QUEUE_ID'])['used_cpu'].diff(1).fillna(0)
train['used_mem_diff1'] = train.groupby(['QUEUE_ID'])['used_mem'].diff(1).fillna(0)
train['used_disk_diff1'] = train.groupby(['QUEUE_ID'])['DISK_USAGE'].diff(1).fillna(0)
train['to_run_jobs_diff1'] = train.groupby(['QUEUE_ID'])['to_run_jobs'].diff(1).fillna(0)
train['launching_diff1'] = train.groupby(['QUEUE_ID'])['LAUNCHING_JOB_NUMS'].diff(1).fillna(0)
test['used_cpu_diff1'] = test.groupby(['QUEUE_ID'])['used_cpu'].diff(1).fillna(0)
test['used_mem_diff1'] = test.groupby(['QUEUE_ID'])['used_mem'].diff(1).fillna(0)
test['used_disk_diff1'] = test.groupby(['QUEUE_ID'])['DISK_USAGE'].diff(1).fillna(0)
test['to_run_jobs_diff1'] = test.groupby(['QUEUE_ID'])['to_run_jobs'].diff(1).fillna(0)
test['launching_diff1'] = test.groupby(['QUEUE_ID'])['LAUNCHING_JOB_NUMS'].diff(1).fillna(0)

"""
处理一下数据的量纲问题方便训练，这里没有采用sklearn的一些标准化方法，因为后续逆变换不方便，可以自行尝试。
"""
train['CPU_USAGE'] = train['CPU_USAGE'] / 100
train['MEM_USAGE'] = train['MEM_USAGE'] / 100
train['LAUNCHING_JOB_NUMS'] = train['LAUNCHING_JOB_NUMS'] / 150
train['DISK_USAGE'] = train['DISK_USAGE'] / 30
train['used_cpu'] = train['used_cpu'] / 1000
train['used_mem'] = train['used_mem'] / 4000
train['to_run_jobs'] = train['to_run_jobs'] / 150
train['used_cpu_diff1'] = train['used_cpu_diff1'] / 1000
train['used_mem_diff1'] = train['used_mem_diff1'] / 4000
train['used_disk_diff1'] = train['used_disk_diff1'] / 30
train['launching_diff1'] = train['launching_diff1'] / 150

test['CPU_USAGE'] = test['CPU_USAGE'] / 100
test['MEM_USAGE'] = test['MEM_USAGE'] / 100
test['LAUNCHING_JOB_NUMS'] = test['LAUNCHING_JOB_NUMS'] / 150
test['DISK_USAGE'] = test['DISK_USAGE'] / 30
test['used_cpu'] = test['used_cpu'] / 1000
test['used_mem'] = test['used_mem'] / 4000
test['to_run_jobs'] = test['to_run_jobs'] / 150
test['used_cpu_diff1'] = test['used_cpu_diff1'] / 1000
test['used_mem_diff1'] = test['used_mem_diff1'] / 4000
test['used_disk_diff1'] = test['used_disk_diff1'] / 30
test['launching_diff1'] = test['launching_diff1'] / 150

"""
定义打算使用的特征名称
"""
feature_columns = ['used_cpu','used_mem','to_run_jobs','used_cpu_diff1','used_mem_diff1', 
                   'launching_diff1','used_disk_diff1','CPU_USAGE','DISK_USAGE','MEM_USAGE','LAUNCHING_JOB_NUMS']

"""
下面提供一个函数用于将原始的时间顺序的数据转化为监督学习的数据，即网络的输入输出形式。
"""
def series_to_superviesed(x_timeseries, y_timeseries, n_memory_step, n_forcast_step, split = None):
    '''
        x_timeseries: input time series data, numpy array, (time_step, features)
        y_timeseries: target time series data,  numpy array, (time_step, features)
        n_memory_step: number of memory step in supervised learning, int = 5
        n_forcast_step: number of forcase step in supervised learning, int = 5
        split: portion of data to be used as train set, float, e.g. 0.8
    '''
    assert len(x_timeseries.shape) == 2, 'x_timeseries must be shape of (time_step, features)'
    assert len(y_timeseries.shape) == 2, 'y_timeseries must be shape of (time_step, features)' 
    
    #输入序列的长度为input_step
    input_step, input_feature = x_timeseries.shape
    #输出序列的长度为output_step
    output_step, output_feature = y_timeseries.shape
    #输入序列和输出序列的长度必须保持一致
    assert input_step == output_step, 'number of time_step of x_timeseries and y_timeseries are not consistent!'

    #5个一组开始构建数据,构建时，索引一直到倒数后10个
    n_RNN_sample=input_step-n_forcast_step-n_memory_step+1
    
    RNN_x=np.zeros((n_RNN_sample,n_memory_step, input_feature))
    RNN_y=np.zeros((n_RNN_sample,n_forcast_step, output_feature))

    for n in range(n_RNN_sample):
        RNN_x[n,:,:]=x_timeseries[n:n+n_memory_step,:]
        RNN_y[n,:,:]=y_timeseries[n+n_memory_step:n+n_memory_step+n_forcast_step,:]
        
    if split != None:
        assert (split <=0.9) & (split >= 0.1), 'split not in reasonable range'
        return RNN_x[:int(split*len(RNN_x))], RNN_y[:int(split*len(RNN_x))],\
               RNN_x[int(split*len(RNN_x))+1:], RNN_y[int(split*len(RNN_x))+1:]
    else:
        return RNN_x, RNN_y, None, None

"""
划分训练数据和测试数据（本地）
"""
X_train = list()
Y_train = list()
X_test = list()
Y_test = list()
for id_ in tqdm(train.QUEUE_ID.unique()):
  tmp = train[train.QUEUE_ID == id_].reset_index(drop=True)
  x_train,y_train,x_test,y_test = series_to_superviesed(np.array(tmp[feature_columns]),np.array(tmp[['CPU_USAGE',
                                                                                                     'LAUNCHING_JOB_NUMS'
                                                                                                    ]]),5,5,split=0.9)
  X_train.append(x_train)
  Y_train.append(y_train)
  X_test.append(x_test)
  Y_test.append(y_test)

X_train = np.vstack(np.array(X_train))  
Y_train = np.vstack(np.array(Y_train))  
X_test = np.vstack(np.array(X_test))  
Y_test = np.vstack(np.array(Y_test))

from tensorflow.keras.layers import BatchNormalization,Dropout,Dense,Activation,Input,TimeDistributed,Conv1D,MaxPooling1D,Flatten,LSTM,Bidirectional
from tensorflow.keras.models import Sequential
from tensorflow.keras.models import Model

"""
CNN-LSTM网络搭建
这里只训练了一个模型，因此 Y_train 只取其中一个维度
"""
n_steps = 5
n_seq = 1
n_features = 11

model = Sequential()
"""
TimeDistributed可以理解
(32, 10, 16)--（32,10,8）
对每个特征向量都执行dense(8)的操作
"""
#Conv1D可以理解成对格式化数据的行进行kernel_size大小的卷积，进行特征提取的过程,特征维度变为了filters，即64
model.add(TimeDistributed(Conv1D(filters=64, kernel_size=3, activation='relu'), input_shape=(None, n_steps, n_features)))
model.add(TimeDistributed(MaxPooling1D(pool_size=3)))
model.add(TimeDistributed(Flatten()))
model.add(Bidirectional(LSTM(units=1,return_sequences=True)))
model.add(TimeDistributed(Dense(5)))
model.compile(optimizer='rmsprop', loss='mse')
model.summary()

X_train = X_train.reshape((X_train.shape[0], n_seq, n_steps, n_features))
#只预测了CPU_USAGE
Y_train = Y_train[:,:,0].reshape(-1,1,5)
model.fit(X_train, Y_train, epochs=1)

"""
可视化lstm关于真实事件序列预测的拟合程度
"""
X_test = X_test.reshape((X_test.shape[0], n_seq, n_steps, n_features))
y_pred = model.predict(X_test).reshape(-1,5,1)
Y_test = Y_test[:,:,0].reshape(-1,5,1)
import matplotlib.pyplot as plt
for n in range(5):
    plt.subplot(5,1,n+1)
    #预测的下5个时刻的数据可视化
    plt.plot(Y_test[700:1100,n,:],'b', label = 'True')
    plt.plot(y_pred[700:1100,n,:],'r', label = 'Predict')
    plt.legend()
plt.show()

CNN使用模板：

 本模板截取自Kaggle比赛“Web Traffic Time Series Forecasting比赛（时序问题）”，其中数据集如下所示：

其中评价指标如下所示：

我们对时间序列数据进行了多个角度的观察，发现我们的每个网页的流量分布都不太一样，有些网页有着比较明显的周期性，有些网页则有着非常明显的突变，还有些网页则有着明显的上升或者下降的趋势。本文的特征工程全部都是中位数，在开源的很多高分kernel中，我们发现中位数能取得非常不错的效果，同时kaggle 的 GM CPMP也通过可视化的形式阐述了SMAPE指标在常见的分布，例如均匀分布等，只需要中位数便可以取得非常不错的效果，所以此处我们便只使用中位数作为我们的特征，最终我们使用的特征包含：

1. 历史上每一周的中位数；

2. 历史上每隔四周的中位数；

3. 历史上每隔七周的中位数；

4. 历史上每隔一年的的中位数；

5. 整个时间序列的中位数；

我们发现很多对中位数进行简单的菲波那切数列的操作就可以在排行榜上取得非常不错的结果，所以我们完全有理由相信NN学出来的基于中位数的组合可以取得更好的结果。如果大家读Discussion的话，就会发现其实很多金牌的选手也是这么想并且是这么做的。好了，废话不多说，直接上NN方案，该方案在线上可以取得不错的效果，除了效果不错外，它还具有一定的可解释性。

我们的NN框架如下，下入中的展示的特征和实际提取的特征有些许区别，此处我们以下图为例阐述该NN框架的优点。我们将滑动7天的数据作为一个模块，将滑动14天的数据作为一个模块，将滑动21天的数据作为一个模块，同时以最后一天的数据作为一个模块；为了更好的提取长时间序列之间的关系，此处我们以Dialted CNN作为中间的处理框架来提取长序列之间的潜在关系，并最后分别输出一个值r1,r2,r3,r4（原先就是一个中位数，所以r4还是一个中位数），最终我们加一个简单的Dense层，并用SMAPE Loss进行优化。

"""
由于框架中并未封装SMAPE的损失函数，这里我们需要自定义该损失函数

"""
def smape_loss(y_true, y_pred):

    diff = K.abs((y_true - y_pred) / K.clip(K.abs(y_true) + K.abs(y_pred), K.epsilon(), None))

    return 200 * K.mean(diff, axis=-1)


"""
开始构建模型，其中这个膨胀率dilation_rate不是太理解
"""
median_7_data_raw         = Input(shape=(Median7_LEN,))

median_28_data_raw        = Input(shape=(Median28_LEN,))

median_49_data_raw        = Input(shape=(Median49_LEN,))

median_365_data_raw       = Input(shape=(Median365_LEN,))

median_median_data_raw    = Input(shape=(MedianMedian_LEN,))

median_7_data      = Reshape((Median7_LEN, 1))(median_7_data_raw)

median_28_data     = Reshape((Median28_LEN, 1))(median_28_data_raw)

median_49_data     = Reshape((Median49_LEN, 1))(median_49_data_raw)

median_365_data    = Reshape((Median365_LEN, 1))(median_365_data_raw)

median_median_data = Reshape((MedianMedian_LEN, 1))(median_median_data_raw)

conv_result1 = Activation('relu')(concatenate([Conv1D(64, 3, padding='same')(median_7_data), Conv1D(64, 3, dilation_rate=3, padding='same')(median_7_data),]))

conv_result1 = Activation('relu')(concatenate([Conv1D(64, 3, padding='same')(conv_result1), Conv1D(64, 3, dilation_rate=3, padding='same')(conv_result1),]))

conv_result1 = Activation('relu')(concatenate([Conv1D(64, 3, padding='same')(conv_result1),Conv1D(64, 3, dilation_rate=3, padding='same')(conv_result1),]))

conv_result1 = (Reshape((conv_result1.shape[1] * conv_result1.shape[2],1,))(conv_result1))

conv_result2 = Activation('relu')(concatenate([Conv1D(64, 3, padding='same')(median_28_data),Conv1D(32, 3, dilation_rate=3, padding='same')(median_28_data),]))

conv_result2 = Activation('relu')(concatenate([Conv1D(64, 3, padding='same')(conv_result2),Conv1D(32, 3, dilation_rate=3, padding='same')(conv_result2),]))

conv_result2 = Activation('relu')(concatenate([Conv1D(64, 3, padding='same')(conv_result2),Conv1D(32, 3, dilation_rate=3, padding='same')(conv_result2),]))

conv_result2 = (Reshape((conv_result2.shape[1] * conv_result2.shape[2],1,))(conv_result2))

conv_result3 = Activation('relu')(Conv1D(32, 3, padding='same')(median_49_data)) 

conv_result3 = Activation('relu')(Conv1D(32, 3, padding='same')(conv_result3)) 

conv_result3 = (Reshape((conv_result3.shape[1] * conv_result3.shape[2],1,))(conv_result3))

res        = concatenate([conv_result1,conv_result2, conv_result3, median_365_data, median_median_data],axis=1)

res        = (Reshape((res.shape[1] * res.shape[2],))(res)) 

result     = Dropout(0.2)(Dense(1024)(res))

result     = Dropout(0.2)(Dense(256)(result))

result     = Dense(62)(result)

model = Model(inputs=[median_7_data_raw, median_28_data_raw,median_49_data_raw,median_365_data_raw,median_median_data_raw], outputs= result)

model.compile(optimizer=Adam(lr=1e-3), loss=smape_loss)

训练模型前corr()进行特征筛选：

m_type = 'lgb'
for label in ['CPU_USAGE_1','LAUNCHING_JOB_NUMS_1','CPU_USAGE_2','LAUNCHING_JOB_NUMS_2','CPU_USAGE_3','LAUNCHING_JOB_NUMS_3',
              'CPU_USAGE_4','LAUNCHING_JOB_NUMS_4','CPU_USAGE_5','LAUNCHING_JOB_NUMS_5']:
    print('############## {} ##############'.format(label))
    
    # 特征选择
    tmp = train_df[features+[label]].corr()
    tmp[label] = tmp[label].apply(lambda x:abs(x))
    feats = tmp.sort_values(label,ascending=False).reset_index()['index'].tolist()[1:150]
    
    x_train = train_df[feats]
    x_test = test_df[feats]
    
    y_train = train_df[label]
    
    lgb_train, lgb_test = cv_model(lgb, x_train, y_train, x_test, m_type)
    
    test_df[label] = lgb_test

多折交叉验证训练模型模板：

def cv_model(clf, train_x, train_y, test_x, clf_name):
    folds = 5
    seed = 2019
    kf = KFold(n_splits=folds, shuffle=True, random_state=seed)

    oof = np.zeros(train_x.shape[0])
    pred = np.zeros(test_x.shape[0])

    cv_scores = []

    for i, (train_index, valid_index) in enumerate(kf.split(train_x, train_y)):
        trn_x, trn_y, val_x, val_y = train_x.iloc[train_index], train_y[train_index], train_x.iloc[valid_index], train_y[valid_index]

        if clf_name == "lgb":
            train_matrix = clf.Dataset(trn_x, label=trn_y)
            valid_matrix = clf.Dataset(val_x, label=val_y)

            params = {
                'boosting_type': 'gbdt',
                'objective': 'regression',
                'metric': 'mse',
                'num_leaves': 2 ** 7,
                'feature_fraction': 0.8,
                'bagging_fraction': 0.8,
                'learning_rate': 0.1,
                'nthread': 28,
                'silent': True,
                'verbose': -1,
            }

            model = clf.train(params, train_matrix, 50000, 
                              valid_sets=[train_matrix, valid_matrix], 
                              categorical_feature=[],
                              verbose_eval=1000,early_stopping_rounds=200)
            val_pred = model.predict(val_x, num_iteration=model.best_iteration)
            test_pred = model.predict(test_x, num_iteration=model.best_iteration)
        oof[valid_index] = val_pred
        pred += test_pred / kf.n_splits
        cv_scores.append(mean_absolute_error(val_y, val_pred))
    
    return oof, pred

----------------------------模型融合-------------------------------

将Dataframe以csv格式存储到本地：

data.to_csv('data_for_tree.csv', index=0)

--------------数据挖掘比赛总结的一些tricks----------------

• 代码做好备份，若使用jupyter notebook来做比赛，将代码名称命名为项目名+日期，每个代码记录下自己做的尝试与取得效果。
• 用excell以日期的形式将调参的历史给记录下来，以便后面回滚。
• 将预测的结果和真实值曲线可视化出来，以便对预测值做一些修正。
• 对分布不一致的数据进行分开建模。