机器学习相关操作分享（四）

# 数据准备阶段
## 导入基本的包
import time
import matplotlib.pyplot as plt
# plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文标签
# plt.rcParams['axes.unicode_minus']=False #用来正常显示负号
# plt.rcParams['font.family'] = ['sans-serif']
import seaborn as sns

import numpy as np
import pandas as pd 
from sklearn.metrics import mean_squared_error,mean_absolute_error

import xgboost as xgb
import lightgbm as lgb
from catboost import CatBoostRegressor

from sklearn.model_selection import KFold,StratifiedKFold
from sklearn.preprocessing import LabelEncoder

pd.set_option('display.max_columns',100)#显示最大列数

import warnings
warnings.filterwarnings("ignore")
## 导入数据
data_path = './data/'
train_data = pd.read_csv(data_path + 'train_dataset.csv')
train_data['type'] = 1
test_data = pd.read_csv(data_path + 'test_dataset.csv')
test_data['type'] = 1
sample_sub = pd.read_csv(data_path + 'submit_example.csv')

print('train：', train_data.shape)
print('test：', test_data.shape)
train_data.head()
# 数据分析
# 基本信息
train_data.info()
# 查看训练集中信用分的统计信息及分布
print(train_data['信用分'].describe())
train_data['信用分'].hist(bins=70)
查阅资料得知，信用评分从统计学上讲，有一个突出的特点：分数特别低的人和分数特别高的人都比较少，大多数人评分中等，大体呈现左偏分布。

本题给的数据也基本符合这个情况，预感两极人群的预测会成为这个题目后期的关键点。
## 单变量分析
features = [f for f in train_data.columns if f not in ['用户编码','信用分']]
for f in features:
    print(f + "的特征分布如下：")
    print(train_data[f].describe())
    if train_data[f].nunique()<20:
        print(train_data[f].value_counts())
    plt.hist(train_data[f], bins=70)
    plt.show()
**根据单变量分析如下**
* 用户年龄：发现290名年龄为0的用户，以及22位100岁以上的用户。推测年龄位0的用户是主办方用0填充了缺失值。
* 用户话费敏感度：在字段说明中，敏感度包含1-5共五级用户。实际用户出现一部分等级为0的用户。推测这些原本也应是缺失值。
* 在一些消费类的账单数据和某些APP的使用次数，数据分布呈现如下长尾形态。

## 多变量分析
主要分析特征与标签的关系
features = [f for f in train_data.columns if f not in ['用户编码','信用分']]
for f in features:
    if train_data[f].nunique()>=20:
        sns.jointplot(x=f,y='信用分',data = train_data)
**根据多变量分析如下**
* 部分特征与信用分存在相关性，用户账单当月总费用（元）、用户当月账户余额（元）等
## 特征之间是否冗余
sns.heatmap(train_data[features].corr(), cmap='Reds')
plt.show()
可以看出，绝大多数特征间线性相关性并不高，最高的为 '用户近6个月平均消费值（元）'和'用户账单当月总费用（元）'，Pearson相关系数达到0.903464，也暂时都保留。
## 查看某个字段不同取值或不同范围，信用分的分布
train_data[train_data['是否大学生客户'] == 1]['信用分'].hist(bins=55)
train_data[train_data['是否大学生客户'] == 0]['信用分'].hist(bins=55)
train_data[train_data['用户话费敏感度'] == 0]['信用分'].hist(bins=55)
train_data[train_data['用户话费敏感度'] == 1]['信用分'].hist(bins=55)
train_data[train_data['用户话费敏感度'] == 2]['信用分'].hist(bins=55)
train_data[train_data['用户话费敏感度'] == 4]['信用分'].hist(bins=55)
train_data.groupby(['用户话费敏感度'])['信用分'].mean()
# 特征工程
## 数据预处理
* 数据合并
* 异常值处理
* 特征分类
data = pd.concat([train_data, test_data], ignore_index=True, sort=True)

data.loc[data['用户年龄'] == 0, '用户年龄'] = None
data.loc[data['用户年龄'] > 100, '用户年龄'] = None
data.loc[data['用户话费敏感度'] == 0, '用户话费敏感度'] = None
data.loc[data['用户近6个月平均消费值（元）'] == 0, '用户近6个月平均消费值（元）'] = None

data.rename(columns={'用户编码': 'id', '信用分': 'score'}, inplace=True)

origin_bool_feature = ['当月是否体育场馆消费', '当月是否景点游览', '当月是否看电影', '当月是否到过福州山姆会员店', '当月是否逛过福州仓山万达',
                       '缴费用户当前是否欠费缴费', '是否经常逛商场的人', '是否大学生客户', '是否4G不健康客户', '是否黑名单客户',
                       '用户最近一次缴费距今时长（月）', '用户实名制是否通过核实']

origin_num_feature = ['用户话费敏感度', '用户年龄', '近三个月月均商场出现次数', '当月火车类应用使用次数', '当月飞机类应用使用次数',
                      '当月物流快递类应用使用次数', '用户当月账户余额（元）', '用户网龄（月）', '缴费用户最近一次缴费金额（元）',
                      '当月通话交往圈人数', '当月旅游资讯类应用使用次数', '当月金融理财类应用使用总次数', '当月网购类应用使用次数',
                      '当月视频播放类应用使用次数', '用户账单当月总费用（元）', '用户近6个月平均消费值（元）']

count_feature_list = []
## 基本统计特征
def feature_count(data, features=[]):
    
    # 样本数等于类别数
    if len(set(features)) != len(features):
        print('equal feature !!!!')
        return data
    
    new_feature = 'count'
    
    # 构建特征名
    for i in features:
        new_feature += '_' + i.replace('add_', '')
    try:
        del data[new_feature]
    except:
        pass
    
    # 构造特征
    temp = data.groupby(features).size().reset_index().rename(columns={0: new_feature})
    data = data.merge(temp, 'left', on=features)
    if new_feature not in count_feature_list:
        count_feature_list.append(new_feature)
    return data

fee_feature = ['用户近6个月平均消费值（元）', '用户账单当月总费用（元）', '缴费用户最近一次缴费金额（元）']

for i in fee_feature:
    data = feature_count(data, [i])
data.groupby('用户账单当月总费用（元）').size().reset_index().sort_values(0, ascending=False)[:20]
与实际业务存在很大的联系，如套餐费用，对数值特征进行编码，通过count来反映套餐类别信息
diff_feature = ['fee_del_mean', 'fee_remain_now']
data['five_all'] = data['用户近6个月平均消费值（元）'] * data['用户网龄（月）'].apply(lambda x: min(x, 6)) - data['用户账单当月总费用（元）']
data['fee_del_mean'] = data['用户账单当月总费用（元）'] - data['用户近6个月平均消费值（元）']
# 缴费对于消费的比例
data['fee_remain_now'] = data['缴费用户最近一次缴费金额（元）'] / data['用户账单当月总费用（元）']
# 各类行为总次数
data['次数'] = data[['当月网购类应用使用次数', '当月物流快递类应用使用次数', '当月金融理财类应用使用总次数',
                   '当月视频播放类应用使用次数', '当月飞机类应用使用次数', '当月火车类应用使用次数', '当月旅游资讯类应用使用次数']].sum(axis=1)
for col in ['当月金融理财类应用使用总次数', '当月旅游资讯类应用使用次数']:  # 这两个比较积极向上一点
    data[col + '_百分比'] = data[col] / data['次数']
data['regist_month'] = data['用户网龄（月）'] % 12
num_feature = ['次数', '当月金融理财类应用使用总次数_百分比',
               '当月旅游资讯类应用使用次数_百分比',
               'five_all',
               'regist_month'] + diff_feature + origin_bool_feature + origin_num_feature + count_feature_list
cate_feature = []
for i in num_feature:
    data[i] = data[i].astype(float)
feature = num_feature + cate_feature
# 训练模型
def get_predict_w(model, data, label='label', feature=[], cate_feature=[], random_state=2018, n_splits=5,
                  model_type='lgb'):
    # 随机数种子
    model.random_state = random_state
    # 交叉验证
    kfold = KFold(n_splits=n_splits, shuffle=True, random_state=random_state)
    
    # 初始化预测结果
    predict_label = 'predict_' + label
    data[predict_label] = 0
    
    # 测试集index提取
    test_index = (data[label].isnull()) | (data[label] == -1)
    # 训练数据提取
    train_data = data[~test_index].reset_index(drop=True)
    # 测试数据提取
    test_data = data[test_index]

    for train_idx, val_idx in kfold.split(train_data):

        train_x = train_data.loc[train_idx][feature]
        train_y = train_data.loc[train_idx][label]

        test_x = train_data.loc[val_idx][feature]
        test_y = train_data.loc[val_idx][label]
        
        if model_type == 'lgb':
            model.fit(train_x, train_y, eval_set=[(train_x, train_y),(test_x, test_y)], early_stopping_rounds=100,
                      eval_metric='mae',
                      categorical_feature=cate_feature,
                      # sample_weight=train_data.loc[train_idx]['sample_weight'],
                      verbose=100)
        elif model_type == 'ctb':
            model.fit(train_x, train_y, eval_set=[(train_x, train_y),(test_x, test_y)], early_stopping_rounds=100,
                      cat_features=cate_feature,
                      # sample_weight=train_data.loc[train_idx]['sample_weight'],
                      verbose=100)
        elif model_type == 'xgb':
            model.fit(train_x, train_y, eval_set=[(train_x, train_y),(test_x, test_y)], early_stopping_rounds=100, 
                      # sample_weight=train_data.loc[train_idx]['sample_weight'],
                      verbose=100)
        
        train_data.loc[val_idx, predict_label] = model.predict(test_x)
        
        # 获取测试集结果
        if len(test_data) != 0:
            test_data[predict_label] = test_data[predict_label] + model.predict(test_data[feature])
    
    # 测试集结果加权平均
    test_data[predict_label] = test_data[predict_label] / n_splits
    
    print(mean_squared_error(train_data[label], train_data[predict_label]) * 5, train_data[predict_label].mean(),
          test_data[predict_label].mean())

    return pd.concat([train_data, test_data], sort=True, ignore_index=True), predict_label
## 单模型
### LightGBM
lgb_model = lgb.LGBMRegressor(
    num_leaves=32, reg_alpha=0., reg_lambda=0.01, objective='mse', metric='mse',
    max_depth=-1, learning_rate=0.1, min_child_samples=50,
    n_estimators=500, subsample=0.7, colsample_bytree=0.45, subsample_freq=5,
)
data.tail()
data, predict_label = get_predict_w(lgb_model, data, label='score',
                                    feature=feature,
                                    random_state=2019, n_splits=5)
data['lgb_mse'] = data[predict_label]
### CatBoost
ctb_params = {
    'n_estimators': 10000,
    'learning_rate': 0.05,
    'random_seed': 4590,
    'reg_lambda': 0.08,
    'subsample': 0.7,
    'bootstrap_type': 'Bernoulli',
    'boosting_type': 'Plain',
    'one_hot_max_size': 10,
    'rsm': 0.5,
    'leaf_estimation_iterations': 5,
    'use_best_model': True,
    'max_depth': 6,
    'verbose': -1,
    'thread_count': 4
}
ctb_model = CatBoostRegressor(**ctb_params)

data, predict_label = get_predict_w(ctb_model, data, label='score',
                                    feature=feature,
                                    random_state=2019, n_splits=5, model_type='ctb')
data['ctb_mse'] = data[predict_label]
### XGBoost
xgb_model = xgb.XGBRegressor(
                        max_depth=6 , learning_rate=0.05, n_estimators=10000, 
                        objective='reg:linear', tree_method = 'hist', subsample=0.8, 
                        colsample_bytree=0.6, min_child_samples=5
                        )
# objective='reg:linear' 线性回归
# 替换inf
data[feature] = data[feature].replace(np.inf, np.nan)
data, predict_label = get_predict_w(xgb_model, data, label='score',
                                    feature=feature,
                                    random_state=2019, n_splits=5, model_type='xgb')

data['xgb_mse'] = data[predict_label]
## 损失函数选择
**MSE均方误差**

因为MSE对error e进行了平方，可以看到，如果e大于1，这个值就会>> |e|。 用了MSE为代价函数的模型因为要最小化这个异常值带来的误差，就会尽量贴近异常值，也就是对outliers（异常值）赋予更大的权重。这样就会影响总体的模型效果。
**MAE平均绝对误差**

相比MSE来说，MAE在数据里有不利于预测结果异常值的情况下撸棒性更好。

可以这么想？哪个常数能够最小化我们的MSE？ 答案是中值。因为在有异常值的时候，中值的代表性要好于均值。所以MAE的撸棒性要高于MSE。
**可以通过不同损失函数的尝试，构建结果的差异性，进行最终的融合**
**LightGBM使用MAE平均绝对误差**
lgb_model = lgb.LGBMRegressor(
    num_leaves=32, reg_alpha=0., reg_lambda=0.01, objective='mae',
    max_depth=-1, learning_rate=0.1, min_child_samples=50,
    n_estimators=500, subsample=0.7, colsample_bytree=0.45, subsample_freq=5,
)

data, predict_label = get_predict_w(lgb_model, data, label='score',
                                    feature=feature,
                                    random_state=2019, n_splits=5)

data['lgb_mae'] = data[predict_label]
print(data['lgb_mse'].describe())
data['lgb_mse'].hist(bins=70)
print(data['lgb_mae'].describe())
data['lgb_mae'].hist(bins=70)
## 加权融合
* 构建有差异的结果
* 模型差异
* 样本差异
* 特征差异
* 损失函数差异
* 训练目标差异，对于树模型而言，更容易学习到稳定的结果，如果目标的值方差很大，可以选择进行log变换
all_score = ['lgb_mse', 'ctb_mse', 'xgb_mse']
data['t_label'] = data['lgb_mse'] * 0.5 + data['ctb_mse'] * 0.3 + data['xgb_mse'] * 0.2
## Stacking
**Stacking是一种表示学习(representation learning)**

stacking集成学习框架的对于基分类器的两个要求: 差异化要大、准确性要高

**Stacking的输出层选择简单的模型，如逻辑回归等**

为了降低过拟合的问题，第二层分类器应该是较为简单的分类器，广义线性如逻辑回归是一个不错的选择。在特征提取的过程中，我们已经使用了复杂的非线性变换，因此在输出层不需要复杂的分类器。
### 第二层仅为学习到的特征
def stack_model(oof_1, oof_2, oof_3, predictions_1, predictions_2, predictions_3, y, eval_type='regression'):
   
    train_stack = np.vstack([oof_1, oof_2, oof_3]).transpose()
    test_stack = np.vstack([predictions_1, predictions_2, predictions_3]).transpose()
    
    from sklearn.model_selection import RepeatedKFold
    folds = RepeatedKFold(n_splits=5, n_repeats=1, random_state=2018)
    oof = np.zeros(train_stack.shape[0])
    predictions = np.zeros(test_stack.shape[0])

    for fold_, (trn_idx, val_idx) in enumerate(folds.split(train_stack, y)):
        print("fold n°{}".format(fold_+1))
        trn_data, trn_y = train_stack[trn_idx], y[trn_idx]
        val_data, val_y = train_stack[val_idx], y[val_idx]
        print("-" * 10 + "Stacking " + str(fold_) + "-" * 10)
        clf = BayesianRidge()
        clf.fit(trn_data, trn_y)

        oof[val_idx] = clf.predict(val_data)
        predictions += clf.predict(test_stack) / (n_splits * n_repeats)
    if eval_type == 'regression':
        print('mean: ',np.sqrt(mean_squared_error(y, oof)))
    if eval_type == 'binary':
        print('mean: ',log_loss(y, oof))
    
    return oof, predictions
# oof_stack , predictions_stack  = stack_model(oof_lgb[0] , oof_xgb[0] , oof_cat[0] , predictions_lgb[0] , predictions_xgb[0] , predictions_cat[0] , target)
### 学习到的特征+原始特征

**Stacking是否需要多层？第一层的分类器是否越多越好？**

stacking的表示学习不是来自于多层堆叠的效果，而是来自于不同学习器对于不同特征的学习能力，并有效的结合起来。一般来看，2层对于stacking足够了。多层的stacking会面临更加复杂的过拟合问题，且收益有限。

第一层分类器的数量对于特征学习应该有所帮助，经验角度看越多的基分类器越好。即使有所重复和高依赖性，我们依然可以通过特征选择来处理，问题不大。
**stacking与深度学习不同之处**

* stacking需要宽度，深度学习不需要

* 深度学习需要深度，而stacking不需要

**但stacking和深度学习都共同需要面临**

* 黑箱与解释问题
* 严重的过拟合问题
## 样本权重
**样本权重参数： sample_weight**
* 第一类：样本不平衡问题

样本不平衡，导致样本不是总体样本的无偏估计，从而可能导致我们的模型预测能力下降。遇到这种情况，我们可以通过调节样本权重来尝试解决这个问题

* 第二类：误差较大的样本

给予误差交大的样本，很难学习的样本更大的权重
### 选择误差大的样本

### 调整样本权重
data['temp_label'] = data['lgb_mse']
data.loc[data.id.isin(ab_id), 'temp_label'] = None
data['sample_weight'] = data['temp_label'] + 200
data['sample_weight'] = data['sample_weight'] / data['sample_weight'].mean()

##top up amount, 充值金额是整数，和小数，应该对应不同的充值途径？
def produce_offline_feature(train_data):
    train_data['不同充值途径']=0
    train_data['不同充值途径'][(train_data['缴费用户最近一次缴费金额（元）']%10==0)&train_data['缴费用户最近一次缴费金额（元）']!=0]=1
    return train_data

train_data=produce_offline_feature(train_data)
test_data=produce_offline_feature(test_data)
##看importance，当月话费 和最近半年平均话费都很高，算一下当月/半年 -->稳定性
def produce_fee_rate(train_data):
    train_data['当前费用稳定性']=train_data['用户账单当月总费用（元）']/(train_data['用户近6个月平均消费值（元）']+1)
    
    ##当月话费/当月账户余额
    train_data['用户余额比例']=train_data['用户账单当月总费用（元）']/(train_data['用户当月账户余额（元）']+1)
    return train_data

train_data=produce_offline_feature(train_data)
test_data=produce_offline_feature(test_data)
#获取特征
def get_features(data):
    data.loc[data['用户年龄']==0,'用户年龄'] = None
    data.loc[data['用户话费敏感度'] == 0, '用户话费敏感度'] = None
    data.loc[data['用户账单当月总费用（元）'] == 0, '用户账单当月总费用（元）'] = None
    data.loc[data['用户近6个月平均消费值（元）'] == 0, '用户近6个月平均消费值（元）'] = None
    data['缴费金额是否能覆盖当月账单']    = data['缴费用户最近一次缴费金额（元）'] - data['用户账单当月总费用（元）']
    data['最近一次缴费是否超过平均消费额'] = data['缴费用户最近一次缴费金额（元）'] - data['用户近6个月平均消费值（元）']
    data['当月账单是否超过平均消费额'] = data['用户账单当月总费用（元）'] - data['用户近6个月平均消费值（元）']
    
    #映射年龄
    def map_age(x):
        if x<=18:
            return 1
        elif x<=30:
            return 2
        elif x<=35:
            return 3
        elif x<=45:
            return 4
        else:
            return 5
    data['是否大学生_黑名单']=data['是否大学生客户']+data['是否黑名单客户']
    data['是否去过高档商场']=data['当月是否到过福州山姆会员店']+data['当月是否逛过福州仓山万达']
    data['是否去过高档商场']=data['是否去过高档商场'].map(lambda x:1 if x>=1 else 0)
    data['是否_商场_电影']=data['是否去过高档商场']*data['当月是否看电影']
    data['是否_商场_体育馆']=data['是否去过高档商场']*data['当月是否体育场馆消费']
    data['是否_商场_旅游']=data['是否去过高档商场']*data['当月是否景点游览']
    data['是否_电影_体育馆']=data['当月是否看电影']*data['当月是否体育场馆消费']
    data['是否_电影_旅游']=data['当月是否看电影']*data['当月是否景点游览']
    data['是否_旅游_体育馆']=data['当月是否景点游览']*data['当月是否体育场馆消费']
    
    data['是否_商场_旅游_体育馆']=data['是否去过高档商场']*data['当月是否景点游览']*data['当月是否体育场馆消费']
    data['是否_商场_电影_体育馆']=data['是否去过高档商场']*data['当月是否看电影']*data['当月是否体育场馆消费']
    data['是否_商场_电影_旅游']=data['是否去过高档商场']*data['当月是否看电影']*data['当月是否景点游览']
    data['是否_体育馆_电影_旅游']=data['当月是否体育场馆消费']*data['当月是否看电影']*data['当月是否景点游览']
    
    data['是否_商场_体育馆_电影_旅游']=data['是否去过高档商场']*data['当月是否体育场馆消费']*\
                                   data['当月是否看电影']*data['当月是否景点游览']
    
    discretize_features=['当月物流快递类应用使用次数','当月飞机类应用使用次数',\
                         '当月火车类应用使用次数','当月旅游资讯类应用使用次数']
    data['交通类应用使用次数比']=(data['当月飞机类应用使用次数'] + 1) / (data['当月火车类应用使用次数'] + 1)
    
    data['6个月平均占比总费用']=data['用户近6个月平均消费值（元）']/data['用户账单当月总费用（元）']+1
    
    def map_discretize(x):
        if x==0:
            return 0
        elif x<=5:
            return 1
        elif x<=15:
            return 2
        elif x<=50:
            return 3
        elif x<=100:
            return 4
        else:
            return 5
        
    for col in discretize_features[:]:
        data[col]=data[col].map(lambda x:map_discretize(x))
    
    return data

train_data=get_features(train_data)
test_data=get_features(test_data)
def base_process(data):
    transform_value_feature=['用户年龄','用户网龄（月）','当月通话交往圈人数','近三个月月均商场出现次数','当月网购类应用使用次数',\
                             '当月物流快递类应用使用次数','当月金融理财类应用使用总次数','当月视频播放类应用使用次数',\
                             '当月飞机类应用使用次数','当月火车类应用使用次数','当月旅游资讯类应用使用次数']
    user_fea=['缴费用户最近一次缴费金额（元）','用户近6个月平均消费值（元）','用户账单当月总费用（元）','用户当月账户余额（元）']
    log_features=['当月网购类应用使用次数','当月金融理财类应用使用总次数','当月物流快递类应用使用次数','当月视频播放类应用使用次数']
    
    #处理离散点
    for col in transform_value_feature+user_fea+log_features:
        #取出最高99.9%值
        ulimit=np.percentile(train_data[col].values,99.9)
        #取出最低0.1%值
        llimit=np.percentile(train_data[col].values,0.1)
        train_data.loc[train_data[col]>ulimit,col]=ulimit
        train_data.loc[train_data[col]<llimit,col]=llimit
        
    for col in user_fea+log_features:
        data[col]=data[col].map(lambda x:np.log1p(x))
    
    return data

# train_data=base_process(train_data)
# test_data=base_process(test_data)
lgb_model = lgb.LGBMRegressor(
    num_leaves=32, reg_alpha=0., reg_lambda=0.01, objective='mse', metric='mae',
    max_depth=-1, learning_rate=0.01, min_child_samples=50,
    n_estimators=15000, subsample=0.7, colsample_bytree=0.45, subsample_freq=5,
)
train_label = train_data['信用分']
train = train_data.drop(['用户编码','信用分'], axis=1)
test = test_data.drop(['用户编码'], axis=1)
folds = KFold(n_splits=5, shuffle=False, random_state=2019)
oof_lgb = np.zeros(train.shape[0])
predictions_lgb = np.zeros(test.shape[0])
for fold_, (trn_idx, val_idx) in enumerate(folds.split(train, train_label)):
    print("fold n°{}".format(fold_+1))
    trn_x, trn_y = train.loc[trn_idx].values, train_label.loc[trn_idx].values
    val_x, val_y = train.loc[val_idx].values, train_label.loc[val_idx].values
    
    lgb_model.fit(trn_x, trn_y, 
              eval_set=[(trn_x, trn_y),(val_x, val_y)], 
              verbose=500, early_stopping_rounds=300)
    
    oof_lgb[val_idx] = lgb_model.predict(train.loc[val_idx].values)
    predictions_lgb += lgb_model.predict(test) / folds.n_splits

print("CV mse score: {:<8.5f}".format(mean_squared_error(train_label , oof_lgb)))
print("CV mae score: {:<8.5f}".format(mean_absolute_error(train_label, oof_lgb)))
sample_sub = sample_sub[['id']]
sample_sub['score'] = predictions_lgb
sample_sub['score'] = sample_sub['score'].apply(lambda x: int(np.round(x)))
sample_sub.to_csv('output/sub.csv', index=False)