Kaggle便利店竞赛，使用xgboost模型销量预测，附完整Python代码

使用商店、促销和竞争对手数据预测销售

Rossmann在7个欧洲国家经营着3000多家药店。目前，Rossmann商店经理的任务是提前6周预测他们的日销售额。商店的销售受到许多因素的影响，包括促销、竞争、学校和国家假日、季节性和地域性。由于数以千计的管理者根据自己的特殊情况预测销售，结果的准确性可能会有很大的差异。在他们的第一次Kaggle竞争中，Rossmann要求预测德国1115家商店的6周日销售额。可靠的销售预测使商店经理能够制定有效的员工时间表，提高生产力和积极性。

训练集、测试集文件

train.csv-包括销售在内的历史数据

test.csv-不包括销售的历史数据

sample_submission.csv-格式正确的示例提交文件

store.csv-有关存储的补充信息

数据字段描述

Store - 每个商店的唯一ID

Sales - 销售额

Customers - 销售客户数

Open - 商店是否营业 0=关闭，1=开业

StateHoliday - 假日。通常所有商店都在国定假日关门. a = 公共假日, b = 复活节假日, c = 圣诞节, 0 = 无

SchoolHoliday - 学校假期

StoreType - 店铺类型: a, b, c, d

Assortment - 产品组合级别: a = 基本, b = 附加, c = 扩展

CompetitionDistance - 距离最近的竞争对手距离（米）

CompetitionOpenSince[Month/Year] - 竞争对手开业年月 year and month of the time the nearest competitor was opened

Promo - 指店铺当日是否在进行促销

Promo2 - 指店铺是否在进行连续促销 0 = 未参与, 1 = 正在参与

Promo2Since[Year/Week] - 商店开始参与Promo2的年和日历周

一.导入数据

<code> 
import
 numpy 
as
 np
import
 pandas 
as
 pd
import
 matplotlib.pyplot 
as
 plt
import
 seaborn 
as
 sns
import
 xgboost 
as
 xgb
from
 time  
import
 time
import
 pickle

 
store = pd.read_csv(
'store.csv'
)
train = pd.read_csv(
'train.csv'
)
test = pd.read_csv(
'test.csv'
)/<code>

<code> 
 
print
(store.head())
print
(train.head())
print
(test.head())
 
 
print
(train.isnull().sum())
 
print
(test.isnull().sum())
 
 
print
(
test
[
test
[
'Open'
].isnull()])
 
test.fillna(1, inplace=True)
 
print
(store.isnull().sum())
 
 
 
store.fillna(0, inplace=True)
 
strain = train[train[
'Sales'
] > 0]
strain.loc[strain[
'Store'
] == 1, [
'Date'
, 
'Sales'
]].plot(x=
'Date'
, y=
'Sales'
, title=
'Store1'
, figsize=(16, 4))
plt.show()/<code>

从图中可以看出店铺的销售额是有周期性变化的，一年中11,12月份销量相对较高，可能是季节因素或者促销等原因

此外从2014年6-9月份的销量来看，6,7月份的销售趋势与8,9月份类似，而我们需要预测的6周在2015年8,9月份，因此我们可以把2015年6,7月份最近6周的1115家店的数据留出作为测试数据，用于模型的优化和验证

<code> 
 
train
 = train[train[
'Sales'
] > 
0
]
 
train = pd.merge(train, store, 
on
=
'Store'
, how=
'left'
)
test = pd.merge(test, store, 
on
=
'Store'
, how=
'left'
)
print(train.
info
())/<code>

二.特征工程

<code>
for
 
data
  
in
 [train, test]:
    # 将时间特征进行拆分和转化
    
data
[
'year'
] = 
data
[
'Date'
].apply(lambda x: x.split(
'-'
)[
0
])
    
data
[
'year'
] = 
data
[
'year'
].astype(int)
    
data
[
'month'
] = 
data
[
'Date'
].apply(lambda x: x.split(
'-'
)[
1
])
    
data
[
'month'
] = 
data
[
'month'
].astype(int)
    
data
[
'day'
] = 
data
[
'Date' 
].apply(lambda x: x.split(
'-'
)[
2
])
    
data
[
'day'
] = 
data
[
'day'
].astype(int)
    # 将
'PromoInterval'
特征转化为
'IsPromoMonth'
特征，表示某天某店铺是否处于促销月，
1
表示是，
0
表示否
    # 提示下：这里尽量不要用循环，用这种广播的形式，会快很多。循环可能会让你等的想哭
    month2str = {
1
: 
'Jan'
, 
2
: 
'Feb'
, 
3
: 
'Mar'
, 
4
: 
'Apr'
, 
5
: 
'May'
, 
6
:  
'Jun'
, 
7
: 
'Jul'
, 
8
: 
'Aug'
, 
9
: 
'Sep'
, 
10
: 
'Oct'
,
                 
11
: 
'Nov'
, 
12
: 
'Dec'
}
    
data
[
'monthstr'
] = 
data
[
'month'
].map(month2str)
    
data
[
'IsPromoMonth'
] = 
data
.apply(
        lambda x: 
0
 
if
 x[
'PromoInterval'
] == 
0
 
else
  
1
 
if
 x[
'monthstr'
] 
in
 x[
'PromoInterval'
] 
else
 
0
, axis=
1
)
    # 将存在其它字符表示分类的特征转化为数字
    mappings = {
'0'
: 
0
, 
'a'
: 
1
, 
'b'
: 
2
, 
'c'
: 
3
, 
'd'
: 
4
}
    
data
[
'StoreType'
].replace(mappings, inplace=True)
    
data
[
'Assortment'
].replace(mappings, inplace=True)
    
data
[
'StateHoliday'
].replace(mappings, inplace=True)

# 删掉训练和测试数据集中不需要的特征
df_train = train.drop([ 
'Date'
, 
'Customers'
, 
'Open'
, 
'PromoInterval'
, 
'monthstr'
], axis=
1
)
df_test = test.drop([
'Id'
, 
'Date'
, 
'Open'
, 
'PromoInterval'
, 
'monthstr'
], axis=
1
)
# 如上所述，保留训练集中最近六周的数据用于后续模型的测试
Xtrain = df_train[
6
 * 
7
 * 
1115
:]
Xtest = df_train[:
6
 * 
7
 * 
1115
]

# 大家从表上可以看下相关性
plt.subplots(figsize=(
24
, 
20
))
sns.heatmap(df_train.corr(), cmap=
'RdYlGn'
, annot=True, vmin=-
0.1
, vmax=
0.1
, center=
0 
)
plt.show()/<code>

<code> 
 
ytrain
 = np.log1p(Xtrain[
'Sales'
])
ytest
 = np.log1p(Xtest[
'Sales'
])

Xtrain
 = Xtrain.drop([
'Sales'
], axis=
1
)
 
Xtest
 = Xtest.drop([
'Sales'
], axis=
1
)/<code>

三.模型构建

<code> 
def
 
rmspe
(y, yhat)
:

    
return
 np.sqrt(np.mean((yhat / y - 
1
) ** 
2
))
def
 
rmspe_xg
(yhat, y)
:

    y = np.expm1(y.get_label())
    yhat = np.expm1(yhat)
    
return
 
'rmspe'
, rmspe(y, yhat)

 
 
params = {
'objective'
: 
'reg:linear'
,
          
'booster'
: 
'gbtree'
,
          
'eta'
: 
0.03
,
           
'max_depth'
: 
10
,
          
'subsample'
: 
0.9
,
          
'colsample_bytree'
: 
0.7
,
          
'silent'
: 
1
,
          
'seed'
: 
10
}
num_boost_round = 
6000
dtrain = xgb.DMatrix(Xtrain, ytrain)
dvalid = xgb.DMatrix(Xtest, ytest)
watchlist = [(dtrain, 
'train'
), (dvalid, 
'eval'
)]/<code>

三.模型训练

<code> 
print(
'Train a XGBoost model'
)
start = time()
gbm = xgb.train(params, dtrain, num_boost_round, evals=watchlist,
                early_stopping_rounds=
100
, feval=rmspe_xg, verbose_eval=
True
)
pickle.dump(gbm, open(
"pima.pickle.dat"
, 
"wb"
))
end = time()
print(
'Train time is {:.2f} s.'
.format(end - start))
'''
Train time is 923.86 s.
训练花费15分钟。。
'''
/<code>

四.结果优化

<code>gbm = pickle.load(open(
"pima.pickle.dat"
, 
"rb"
))
 
print(
'validating'
)
Xtest.sort_index(inplace=
True
)
ytest.sort_index(inplace=
True
)
yhat = gbm.predict(xgb.DMatrix(Xtest))
error = rmspe(np.expm1(ytest), np.expm1(yhat))
print(
'RMSPE: {:.6f}'
.format(error))
'''
validating
RMSPE: 0.128683
'''

 
res = pd.DataFrame(data=ytest)
res[
'Predicition'
] = yhat
res = pd.merge(Xtest, res, left_index=
True
, right_index=
True
)
res[
'Ratio'
] = res[
'Predicition'
] / res[
'Sales'
]
res[
'Error'
] = abs(res[
'Ratio'
] - 
1
)
res[
'Weight'
] = res[
'Sales'
] / res[
'Predicition'
]
res.head()

 
col_1 = [
'Sales'
, 
'Predicition'
]
col_2 = [
'Ratio'
]
L = np.random.randint(low=
1
, high=
1115 
, size=
3
)
print(
'Mean Ratio of predition and real sales data is {}:store all'
.format(res[
'Ratio'
].mean()))
for
 i 
in
 L:
    s1 = pd.DataFrame(res[res[
'Store'
] == i], columns=col_1)
    s2 = pd.DataFrame(res[res[
'Store'
] == i], columns=col_2)
    s1.plot(title=
'Comparation of predition and real sales data:store {}'
.format(i), figsize=(
12
, 
4
))
    s2.plot(title=
'Ratio of predition and real sales data: store {}'
.format(i), figsize=(
12
, 
4
))
    print(
'Mean Ratio of predition and real sales data is {}:store {}'
.format(s2[
'Ratio'
].mean(), i))


 
res.sort_values([
'Error'
], ascending=
False
, inplace=
True
)
print(res[:
10
])
/<code>

五.模型优化

<code> 
 
print(
'weight correction'
)
W=[(
0.990
+(i/
1000
)) 
for
 i 
in
 range(
20
)]
S=[]
for
 w 
in
 W:
    error=rmspe(np.expm1(ytest),np.expm1(yhat*w))
    print(
'RMSPE for {:.3f}:{:.6f}'
.format(w,error))
    S.append(error)
Score=pd.Series(S,index=W)
Score.plot()
BS=Score[Score.values==Score.values.min()]
print(
'Best weight for Score:{}'
.format(BS))
'''
weight correction
RMSPE for 0.990:0.131899
RMSPE for 0.991:0.129076
RMSPE for 0.992:0.126723
……
Best weight for Score:0.996    0.122779
dtype: float64
'''
 
 
plt.show()/<code>

<code> 
L=range(
1115
)
W_ho=[]
W_test=[]
for
 i 
in
 L:
    s1=pd.DataFrame(res[res[
'Store'
]==i+
1
],columns=col_1)
    s2=pd.DataFrame(df_test[df_test[
'Store'
]==i+
1
])
    W1=[(
0.990
+(i/
1000
)) 
for
 i 
in
 range(
20
)]
    S=[]
    
for
 w 
in
 W1:
        error=rmspe(np.expm1(s1[
'Sales'
]),np.expm1(s1[
'Predicition'
]*w))
        S.append(error)
    Score=pd.Series(S,index=W1)
    BS=Score[Score.values==Score.values.min()]
    a=np.array(BS.index.values)
    b_ho=a.repeat(len(s1))
    b_test=a.repeat(len(s2))
    W_ho.extend(b_ho.tolist())
    W_test.extend(b_test.tolist())
 
Xtest=Xtest.sort_values(by=
'Store'
)
Xtest[
'W_ho'
]=W_ho
Xtest=Xtest.sort_index()
W_ho=list(Xtest[
'W_ho'
].values)
Xtest.drop([
'W_ho'
],axis=
1
,inplace=
True
)

df_test=df_test.sort_values(by=
'Store'
)
df_test[
'W_test' 
]=W_test
df_test=df_test.sort_index()
W_test=list(df_test[
'W_test'
].values)
df_test.drop([
'W_test'
],axis=
1
,inplace=
True
)

 
yhat_new=yhat*W_ho
error=rmspe(np.expm1(ytest),np.expm1(yhat_new))
print(
'RMSPE for weight corretion {:.6f}'
.format(error))
'''
RMSPE for weight corretion 0.116168
相对于整体校正的0.122779的得分又有不小的提高
'''
/<code>