最近最近做了个阿里天池的竞赛项目,有关二手车交易价格预测。一般的数据挖掘项目包括探索性数据分析、特征工程、模型选择、模型调优和模型融合。对数据挖掘过程不熟悉的同学可以参考【绝对干货】机器学习模型训练全流程!,比较详细介绍了机器学习模型构造的全流程。
一、探索性数据分析
探索性数据分析的目的,就是先对数据有个整体的把握,熟悉数据,为接下来的特征工程和 模型构造打下基础。探索性数据分析一般会从数据类型、变量的数值分布、缺失值情况和变量之间的相关关系来分析数据。首先,我们导入数据,数据是官方给定的,分为训练集和测试集两个部分。
#导入需要的包
import pandas as pd
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import seaborn as sns
from operator import itemgetter
import datetime
%matplotlib inline
#导入数据
df_train_data=pd.read_csv('used_car_train_20200313.csv',sep=' ')
df_test_data=pd.read_csv('used_car_testB_20200421.csv',sep=' ')
print(df_train_data.shape)
print(df_test_data.shape)
训练集包含150000条数据,31个变量;测试集包含50000条数据,30个变量(缺少交易价格一列)。值得注意的是,探索性数据分析是为接下来的建模做准备的,所以我们在进行探索数据分析时,只能针对训练集,不然会产生数据窥探偏误。测试集的作用是测试我们基于训练集构建的模型的泛化能力,我们的模型要完全独立于测试集进行构造。这就好比我们比赛,你不能既当球员,又当裁判。
1. 训练集数据类型和缺失值分析
#数据备份
df_train=df_train_data.copy()
#查看总体数据
df_train.head(5)
#查看训练集的变量数据类型
df_train.info()
存在缺失值的列有:model、 bodyType、 fuelType、 gearbox。唯一的有个字符串类型的变量为notRepairedDamage。
#查看唯一的字符串类型变量的数值分布
df_train['notRepairedDamage'].value_counts()
训练集各变量的数值分布
#绘制所有变量的柱形图,查看数据
df_train.hist(bins=50,figsize=(20,15))
plt.cla() #清除axes
从上面的柱形图来看,训练数据种分布有点异常的变量有:creatDate、offerType、seller。下面我们具体看一下这些变量数值分布:
#查看训练集中offerType变量
df_train['offerType'].value_counts()
#查看训练集seller变量
df_train['seller'].value_counts()
这些变量的数值分布可以总结为:
- creatDate为汽车开始售卖时间,训练集中的数据大致在2016-03-09至2016-04-03之间(excel打开文件查看)。这个变量可以考虑减去汽车注册日期regDate,得到汽车的使用时长。汽车使用时长可能时影响汽车销售价格的一个重要变量。
- 训练集和测试集中,offerType变量的值全为0。在构建模型时,可以考虑去掉这个变量。
- .训练集中,seller变量的值只有一个值为1,其余的都为0;而测试集中,全为0,可以考虑去掉这个变量。
3. 相关性分析
我们把字符串类型的变量、以及一些无关的变量去掉,来进行相关性分析。
#获得需要的列名
numeric_columns=df_train.select_dtypes(exclude='object').columns
columns=[col for col in numeric_columns if col not in ['SaleID', 'name', 'regDate',
'regionCode','seller', 'offerType', 'creatDate']]
#根据列名提取数据
train_set=df_train[columns]
#计算各列于交易价格的相关性
corr_matrix=train_set.corr()
corr_matrix['price'].sort_values()
二、特征工程
特征工程主要涉及数据的处理、特征的构造和特征帅选。
- 特征构造
我们对照官方给的变量的定义,发现model、 brand、 bodyType、 fuelType等为取几个值的离散变量。当我们采用基于树模型(XGBoost、LightGBM)来解决回归问题时,针对每一棵树,需要做的是确定分裂的变量,以及分裂点。一般采用启发式方法,根据分裂后的收益(误差)是否下降,并选择下降最多的分裂点,来逐步构建一棵树。由于上述的几个变量的值,只是代表类别,并不代表实际的大小,根据值大小的分裂并不能代表实际的意义。所以我们准备把每一个类别切分出来,形成非0即1的二分变量。
我们首先查看,每个变量属性值的数量:
len(df_train['model'].unique()) #结果为249
len(df_train['brand'].unique()) #结果为40
len(df_train['bodyType'].unique()) #结果为9
len(df_train['fuelType'].unique()) #结果为8
这些属性值种,还包含着NAN的类别,由于我在编程的时候,一直识别NAN来作为判断条件,所以我首先将数据种的NAN填充为-1。
#将nan填充为-1
df_train=df_train.fillna(-1)
(1)处理model
#获得model所有属性值,除nan
model_list=[model for model in list(df_train['model'].unique()) if model != -1]
#将所有的属性值切分出来,用1和0表示是否属于该属性,并保留nan值
all_info={}
for n in range(len(model_list)):
print(n)
info=[]
for i in range(len(df_train['model'])):
if df_train['model'][i]==model_list[n]:
info.append(1)
elif df_train['model'][i]==-1:
info.append(np.nan)
else:
info.append(0)
all_info['model'+str(model_list[n])]=info
#将切分出来的属性保存为csv文件
df_info=pd.DataFrame(all_info)
df_info.to_csv('model.csv',encoding='utf_8_sig')
(2)处理brand、 bodyType和fuelType
#获得除nan外的所有属性值列表
brand_list=[brand for brand in list(df_train['brand'].unique()) if brand != -1]
body_list=[bodyType for bodyType in list(df_train['bodyType'].unique()) if bodyType != -1]
fuel_list=[fuelType for fuelType in list(df_train['fuelType'].unique()) if fuelType != -1]
#将brand所有的属性值切分出来,用1和0表示是否属于该属性,并保留nan值
all_info={}
for n in range(len(brand_list)):
print(n)
info=[]
for i in range(len(df_train['brand'])):
if df_train['brand'][i]==brand_list[n]:
info.append(1)
elif df_train['brand'][i]==-1:
info.append(np.nan)
else:
info.append(0)
all_info['brand'+str(brand_list[n])]=info
#将bodyType所有的属性值切分出来,用1和0表示是否属于该属性,并保留nan值
for n in range(len(body_list)):
print(n)
info=[]
for i in range(len(df_train['bodyType'])):
if df_train['bodyType'][i]==body_list[n]:
info.append(1)
elif df_train['bodyType'][i]==-1:
info.append(np.nan)
else:
info.append(0)
all_info['bodyType'+str(body_list[n])]=info
#将fuelType所有的属性值切分出来,用1和0表示是否属于该属性,并保留nan值
for n in range(len(fuel_list)):
print(n)
info=[]
for i in range(len(df_train['fuelType'])):
if df_train['fuelType'][i]==fuel_list[n]:
info.append(1)
elif df_train['fuelType'][i]==-1:
info.append(np.nan)
else:
info.append(0)
all_info['fuelType'+str(fuel_list[n])]=info
#将brand、bodyType和fuelType的结果输出为csv文件
df_info=pd.DataFrame(all_info)
df_info.to_csv('brand_body_fuel.csv',encoding='utf_8_sig')
(3)字符串类型变量处理
我们在探索数据集时,看到notRepairedDamage列为字符串类型,也从它的数值分布,看到它实质上为二分类变量,只是数据中出现了缺失值,使得为字符串类型,这里我们把缺失值设为nan,将数据类型转化为整型。
#处理字符串属性的列
new_col=[]
for n in range(len(df_train['notRepairedDamage'])):
try:
new_col.append(int(float(df_train['notRepairedDamage'][n])))
except:
new_col.append(np.nan)
df_train['notRepairedDamage_new']=new_col
(4)构造使用时间
在探索数据时,看到creatDate为汽车上线日期,regDate为汽车注册日期,我们可以用这两个日期来构造汽车使用时间,者很可能是影响汽车价格的重要变量。
#构建使用时间
df_train['used_time'] = (pd.to_datetime(df_train['creatDate'], format='%Y%m%d', errors='coerce') -
pd.to_datetime(df_train['regDate'], format='%Y%m%d', errors='coerce')).dt.days
- 特征帅选
在特征构造中,我们引入了很多变量,但并不是所有变量对最终构建回归模型都是有用的。如果将所有变量都放入模型中,一方面会耗费更多时间,另一方面也达不到很好的效果。所有这里将进行特征帅选,具体步骤为先通过相关性分析,删除掉与交易价格相关性很弱的变量(因为构建的许多变量和交易价格都是非常弱的相关,所以先进行这一步可以排除掉大量的变量,这可以减少后面根据模型来做特征帅选的时间),接着再根据 lightGBM模型中的变量重要性来去除掉低重要性的变量。
(1)根据相关性进行帅选
#导入构建的变量数据
df_model_data=pd.read_csv('model.csv')
df_other_data=pd.read_csv('brand_body_fuel.csv')
#加入交易价格一列
df_model=df_model_data.copy()
df_other=df_other_data.copy()
df_model['price']=df_train['price']
df_other['price']=df_train['price']
#去掉空值行,以计算相关系数
df_model.dropna(inplace=True)
df_other.dropna(inplace=True)
#计算相关系数
cor=[]
variable=[]
for column in df_model.columns:
cor.append(np.corrcoef(df_model[column],df_model['price'])[0,1])
variable.append(column)
for column in df_other.columns:
cor.append(np.corrcoef(df_other[column],df_other['price'])[0,1])
variable.append(column)
#导出文件
info={}
info['特征']=variable
info['相关性']=cor
info_df=pd.DataFrame(info)
info_df.to_csv('corcoef.csv',encoding='utf_8_sig')
接下来根据相关性来帅选变量,因为这是第一步的特征帅选,后面还有根据模型的帅选,所以我把相关性的阈值设置为较小的0.1。相关性的绝对值小于0.1,将会被删除。未被删除的变量将会形成一个新的DataFrame.
#根据相关系数帅选变量,阈值0.1
info_df=pd.read_csv('corcoef.csv')
df_new_var=pd.DataFrame(df_train['SaleID'])
for n in range(len(info_df['相关性'])):
corcoef=info_df['相关性'][n]
if np.abs(corcoef)>=0.1:
new_variable=info_df['特征'][n]
if new_variable in df_model_data.columns:
df_new_var[new_variable]=df_model_data[new_variable]
elif new_variable in df_other_data.columns:
df_new_var[new_variable]=df_other_data[new_variable]
#将未被删除的新构造的变量加入到数据中
df_train=df_train.merge(df_new_var, how='left',on='SaleID')
(2)根据lightGBM模型进行帅选
我们可以基于树模型的机器学习模型,求取特征重要度,来进行特征帅选。这里选用lightGBM,具体的过程我借鉴了这篇文章里的内容:一个Python特征选择工具,助力实现高效机器学习
#切分特征和标签
train_set=df_train.copy()
x_train=train_set[[col for col in train_set.columns if col not in ['SaleID','name','regDate','notRepairedDamage','regionCode','creatDate',
'price','model','bodyType','fuelType','gearbox','brand']]]
y_train=train_set['price']
接下来根据模型来帅选特征。
features = pd.get_dummies(x_train)
feature_names = list(features.columns)
features = np.array(features)
labels = np.array(y_train).reshape((-1, ))
feature_importance_values = np.zeros(len(feature_names))
task='regression'
early_stopping=True
eval_metric= 'l2'
n_iterations=10
for _ in range(n_iterations):
if task == 'classification':
model = lgb.LGBMClassifier(n_estimators=1000, learning_rate = 0.05, verbose = -1)
else task == 'regression':
model = lgb.LGBMRegressor(n_estimators=1000, learning_rate = 0.05, verbose = -1)
else:
raise ValueError('Task must be either "classification" or "regression"')
# If training using early stopping need a validation set
if early_stopping:
train_features, valid_features, train_labels, valid_labels = train_test_split(features, labels, test_size = 0.15)
# Train the model with early stopping
model.fit(train_features, train_labels, eval_metric = eval_metric,eval_set = [(valid_features, valid_labels)],early_stopping_rounds = 100, verbose = -1)
gc.enable()
del train_features, train_labels, valid_features, valid_labels
gc.collect()
else:
model.fit(features, labels)
# Record the feature importances
feature_importance_values += model.feature_importances_ / n_iterations
feature_importances = pd.DataFrame({'feature': feature_names, 'importance': feature_importance_values})
#按照重要性大小对特征进行排序
feature_importances = feature_importances.sort_values('importance', ascending = False).reset_index(drop = True)
#计算特征的相对重要性,全部特征的相对重要性之和为1
feature_importances['normalized_importance'] = feature_importances['importance'] / feature_importances['importance'].sum()
#计算特征的累计重要性
feature_importances['cumulative_importance'] = np.cumsum(feature_importances['normalized_importance'])
#选取累计重要性大于0.99的特征,这些特征将会被删除掉。
drop_columns=list(feature_importances.query('cumulative_importance>0.99')['feature'])
接下来去掉要删除的列,并查看新数据集的总体情况:
#去掉重要度低的列
x_set=x_train.copy()
x_set.drop(drop_columns,axis=1,inplace=True)
#对数据集总体概览
x_set.info()
看到新数据集中,缺失值主要集中在我们新构造的变量上,一般处理缺失值的方式为:
- 删除有缺失值的行
- 删除有缺失值的列
- 将缺失值设置为某个值
对于user_time、notRepairedDamage_new两列的缺失值比例较大,删除有缺失值的行将会损失掉比较多的数据。由于新数据中各列通过了特征重要性的帅选,直接删除掉列,将会损失掉重要的特征,特别是user_time变量,重要性程度为第一。所以这里,采用将缺失值设置为某个值的方法,将缺失值设置为1。
x_set.fillna(1,inplace=True)
三、模型选择
首先,我们先训练一个线性模型,来对数据进行初探:
from sklearn.linear_model import LinearRegression
lin_reg=LinearRegression()
lin_reg.fit(x_set,y_train)
为了评判模型的拟合优度,这里先定义一计算模型拟合效果的指标(包括MAR、MSE、RMSE)函数:
from sklearn.metrics import mean_squared_error,mean_absolute_error
def model_goodness(model,x,y):
prediction=model.predict(x)
mae=mean_absolute_error(y,prediction)
mse=mean_squared_error(y,prediction)
rmse=np.sqrt(mse)
print('MAE:',mae)
print('MSE:',mse)
print('RMSE:',rmse)
另外,除了模型的拟合效果外,模型的泛化能力也是评判模型的重要指标,甚至是更重要的指标。我们采用交叉验证的方法来验证模型的泛化能力,这里也先定义一个模型泛化能力的指标计算函数:
from sklearn.model_selection import cross_val_score
def display_scores(scores):
print("Scores:", scores)
print("Mean:", scores.mean())
print("Standard deviation:", scores.std())
我们先来看看,构造的线性模型的拟合效果:
训练数据中,交易价格的均值为5923,线性模型的MAE为2639,看起来并不是一个很好的模型。接下来,我尝试一个更优的模型,随机森林。
from sklearn.ensemble import RandomForestRegressor
forest_reg=RandomForestRegressor()
forest_reg.fit(x_set,y_train)
我们看看模型的拟合效果:
有了非常大的改进,但是该模性可能存在过拟合,我们来利用10折交叉验证来验证模型的泛化性能,这里我们采用负的MAE作为得分:
scores=cross_val_score(forest_reg,x_set,y_train,scoring='neg_mean_absolute_error',cv=10)
mae_scores=np.abs(-scores)
display_scores(mae_scores)
10次平均的MAE为655,虽然比构造的模型MAE大,但其值较小,是一个不错的模型。为了选择好的模型,我们进一步构建GBDT、XGBoost和LightGBM模型。
1.GBDT
from sklearn.ensemble import GradientBoostingRegressor
gbrt=GradientBoostingRegressor()
gbrt.fit(x_set,y_train)
scores=cross_val_score(gbrt,x_set,y_train,scoring='neg_mean_absolute_error',cv=10)
mae_scores=np.abs(scores)
display_scores(mae_scores)
GBDT整体上不如随机森林模型,但是GBDT模型的10次交叉验证的平均MAE得分和模型的MAE很相近,模型的泛化能力较强。
2.XGBoost
import lightgbm as lgb
import xgboost as xgb
xgb_reg= xgb.XGBRegressor()
xgb_reg.fit(x_set,y_train)
scores=cross_val_score(xgb_reg,x_set,y_train,scoring='neg_mean_absolute_error',cv=10)
mae_scores=np.abs(scores)
display_scores(mae_scores)
XGBoost模型整体好于GBDT 模型,并且泛化能力也很好。XGBoost是对GBDT的改进,其性能一般较GBDT更优,感兴趣的可以看看:灵魂拷问,你看过Xgboost原文吗?
3.LightGBM
lgb_reg=lgb.LGBMRegressor()
lgb_reg.fit(x_train,y_train)
scores=cross_val_score(lgb_reg,x_set,y_train,scoring='neg_mean_absolute_error',cv=10)
mae_scores=np.abs(scores)
display_scores(mae_scores)
四、模型调优
在模型选择部分,我们对比了不同的模型,发现随机森林、XGBoost和LightGBM三个模型较优,因此我们选择这几个模型来构建预测模型,并进一步进行模型调优,也就是超参数选择。超参数调整的方法优很多,有手动,有自动。这里我们采用自动调参的方法,自动调参现在一般有网格搜索、随机搜索和贝叶斯方法。网格搜索是对所以的参数空间进行搜索,速度一般很慢,但能得到参数空间中的最优参数。随机搜索速度较快,但可能会错过最优参数。贝叶斯方法是利用贝叶斯观点来进行参数调优,具体可以参考Python 环境下的自动化机器学习超参数调优,它可以保证速度的情况下,尽可能搜索到最优超参数组合。由于网格搜索实在是太慢了,这里我采用随机搜索方法和贝叶斯方法两种来进行超参数搜索。
1.利用随机搜索对随机森林模型进行调优
我们利用sklearn.model_selection模块中的RandomizedSearchCV来进行随机搜索,搜索的超参数包括bootstrap,最大特征数max_features,树的最大深度max_depth,n_estimators。
from sklearn.model_selection import RandomizedSearchCV
#2.设置参数空间
from hyperopt import hp
space_forest = {
'bootstrap':[True,False],
'max_features':list(range(0,25,1)),
'max_depth': list(range(0, 100, 1)),
'n_estimators': list(range(30, 150, 1))
}
#随机搜索,利用5折交叉验证得分来作为模型优劣的判断标准
forest_reg=RandomForestRegressor()
random_search=RandomizedSearchCV(forest_reg, space_forest,cv=5,scoring='neg_mean_squared_error')
#得到最优参数
random_search.best_params_
2.利用贝叶斯方法对LightBoost进行调优
python中的hypreopt包可以进行贝叶斯方法的调优,这篇文章里Python 环境下的自动化机器学习超参数调优,有详细的介绍。
#2.定义参数空间
from hyperopt import hp
space = {
'num_leaves': hp.quniform('num_leaves', 30, 150, 1),
'learning_rate': hp.loguniform('learning_rate', np.log(0.01), np.log(0.2)),
'subsample_for_bin': hp.quniform('subsample_for_bin', 20000, 300000, 20000),
'max_depth': hp.quniform('max_depth', 0, 100, 1),
'n_estimators': hp.quniform('n_estimators', 30, 150, 1)
}
#定义优化函数,即为5折交叉验证的得分
from sklearn.model_selection import cross_val_score
def objective(params, n_folds=5):
num_leaf=int(params['num_leaves'])
estimator=int(params['n_estimators'])
rate=params['learning_rate']
sub_for_bin=int(params['subsample_for_bin'])
max_dep=int(params['max_depth'])
lgb_reg=lgb.LGBMRegressor(num_leaves=num_leaf,n_estimators = estimator,learning_rate=rate,subsample_for_bin=sub_for_bin,max_depth=max_dep)
lgb_reg.fit(x_set,y_train)
scores=cross_val_score(lgb_reg,x_set,y_train,scoring='neg_mean_absolute_error',cv=5)
mae_scores=np.abs(scores)
loss=mae_scores.mean()
return loss
#寻找到使优化函数最小超参数组合,利用hyperopt中的fmin来求最小化
from hyperopt import Trials,fmin,tpe
best = fmin(fn = objective, space = space, algo = tpe.suggest, max_evals = 500)
最小化的交叉验证得分为586,好于未调优的模型。
五、模型融合
模型融合是对调优后的模型的结果进行融合,以提高预测效果。模型融合的方式有:
- 简单加权融合
- stacking/blending
- boosting/bagging(在xgboost,Adaboost,GBDT中已经用到)
这里我们采用简单加权融合。
#构造调优后的模型
forest_reg=RandomForestRegressor(n_estimators=73, max_features= 11, max_depth=33,bootstrap=False)
xgb_reg=xgb.XGBRegressor(colsample_bytree=0.8400547784951151,gamma=0.24562245606367,learning_rate=0.07356238705262594,max_depth=13,n_estimators=135)
lgb_reg=lgb.LGBMRegressor(learning_rate=0.1281370851088587,max_depth=24,n_estimators=148,num_leaves=146,subsample_for_bin=200000)
#将训练集切分出验证集
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error, mean_absolute_error
x_tr,x_val,y_tr,y_val = train_test_split(x_set,y_train,test_size=0.3)
#模型训练
forest_reg.fit(x_tr,y_tr)
xgb_reg.fit(x_tr,y_tr)
lgb_reg.fit(x_tr,y_tr)
#用验证集进行预测
forest_pre=forest_reg.predict(x_val)
xgb_pre=xgb_reg.predict(x_val)
lgb_pre=lgb_reg.predict(x_val)
#输出预测结果的MAE
MAE_forest=mean_absolute_error(forest_pre,y_val)
MAE_xgb=mean_absolute_error(xgb_pre,y_val)
MAE_lgb=mean_absolute_error(lgb_pre,y_val)
print('MAE for RandomForest:',MAE_forest)
print('MAE for xgb:',MAE_xgb)
print('MAE for lgb:',MAE_lgb)
#采用加权融合模型
val_Weighted = (1-2*MAE_lgb/(MAE_xgb+MAE_lgb+MAE_forest))*lgb_pre+(1-2*MAE_xgb/(MAE_xgb+MAE_lgb+MAE_forest))*xgb_pre+(1-2*MAE_forest/(MAE_xgb+MAE_lgb+MAE_forest))*forest_pre
我们看到模型融合的结果好于,单个模型的效果。接下来,我们可以将上面过程中确定的最优模型和模型的融合方式应用到测试集,完成整个的项目。