python中构建灰色预测模型 灰色预测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的年和日历周

一.导入数据

# 导入需要的库import numpy as npimport pandas as pdimport matplotlib.pyplot as pltimport seaborn as snsimport xgboost as xgbfrom time import timeimport pickle# 导入数据集store = pd.read_csv('store.csv')train = pd.read_csv('train.csv')test = pd.read_csv('test.csv')

# 1.检查数据# 可以看前几行观察下数据的基本情况print(store.head())print(train.head())print(test.head())# 对缺失数据检查# 检查训练集数据有无缺失print(train.isnull().sum())# 检查测试集数据有无缺失print(test.isnull().sum())# 发现open有11行缺失值# 看看缺哪些，通过查看train里622号店的营业情况发现，622号店周一到周六都是营业的print(test[test['Open'].isnull()])# 所以我们认为缺失的部分是应该正常营业的，用1填充test.fillna(1, inplace=True)# 接下来看看store的缺失值,store的缺失值比较多print(store.isnull().sum())# 下面对缺失数据进行填充# 店铺竞争数据缺失，而且缺失的都是对应的。原因不明，而且数量也比较多，如果用中值或均值来填充，有失偏颇。暂且填0，解释意义就是刚开业# 店铺促销信息的缺失是因为没有参加促销活动，所以我们以0填充store.fillna(0, inplace=True)# 分析店铺销量随时间的变化strain = train[train['Sales'] > 0]strain.loc[strain['Store'] == 1, ['Date', 'Sales']].plot(x='Date', y='Sales', , figsize=(16, 4))plt.show()

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

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

# 合并数据# 我们只需要销售额大于0的数据train = train[train['Sales'] > 0]# 把store基本信息合并到训练和测试数据集上train = pd.merge(train, store, on='Store', how='left')test = pd.merge(test, store, on='Store', how='left')print(train.info())

二.特征工程

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()

# 提取后续模型训练的数据集# 拆分特征与标签，并将标签取对数处理ytrain = np.log1p(Xtrain['Sales'])ytest = np.log1p(Xtest['Sales'])Xtrain = Xtrain.drop(['Sales'], axis=1)Xtest = Xtest.drop(['Sales'], axis=1)

三.模型构建

# 定义评价函数，可以传入后面模型中替代模型本身的损失函数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 = 6000dtrain = xgb.DMatrix(Xtrain, ytrain)dvalid = xgb.DMatrix(Xtest, ytest)watchlist = [(dtrain, 'train'), (dvalid, 'eval')]

三.模型训练

# 模型训练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分钟。。'''

四.结果优化

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))'''validatingRMSPE: 0.128683'''# 构建保留数据集预测结果res = pd.DataFrame(data=ytest)res['Predicition'] = yhatres = 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(.format(i), figsize=(12, 4))    s2.plot(.format(i), figsize=(12, 4))    print('Mean Ratio of predition and real sales data is {}:store {}'.format(s2['Ratio'].mean(), i))# 分析偏差最大的10个预测结果res.sort_values(['Error'], ascending=False, inplace=True)print(res[:10])#从分析结果来看，初始模型已经可以比较好的预测保留数据集的销售趋势，但相对真实值，模型的预测值整体要偏高一些。从对偏差数据分析来看，偏差最大的3个数据也是明显偏高。因此，我们可以以保留数据集为标准对模型进行偏差校正。

五.模型优化

#模型优化#偏差整体校正优化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 correctionRMSPE for 0.990:0.131899RMSPE for 0.991:0.129076RMSPE for 0.992:0.126723……Best weight for Score:0.996    0.122779dtype: float64'''#当校正系数为0.996时，保留数据集的RMSPE得分最低：0.122779,相对于初始模型0.128683得分有很大的提升。#因为每个店铺都有自己的特点，而我们设计的模型对不同的店铺偏差并不完全相同，所以我们需要根据不同的店铺进行一个细致的校正。plt.show()

#细致校正：以不同的店铺分组进行细致校正，每个店铺分别计算可以取得最佳RMSPE得分的校正系数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_hoXtest=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_testdf_test=df_test.sort_index()W_test=list(df_test['W_test'].values)df_test.drop(['W_test'],axis=1,inplace=True)#计算校正后整体数据的RMSPE得分yhat_new=yhat*W_hoerror=rmspe(np.expm1(ytest),np.expm1(yhat_new))print('RMSPE for weight corretion {:.6f}'.format(error))'''RMSPE for weight corretion 0.116168相对于整体校正的0.122779的得分又有不小的提高'''