前言
在刚刚结束的kaggle比赛M5 Forecasting - Accuracy中,因为是第一次参加,笔者也是花了大量的时间和精力在上面,历时4个月,最终拿到一块银牌(所以我拿到了大学第一个考试挂科。。。求求电磁场老师高抬贵手给点平时分放我一马吧。。。早上看到kaggle成绩异常兴奋,中午考完直接爆炸=。=),当然运气占了很大因素,这次比赛private leadboard的shake up非常大,排名波动几千名都存在的。笔者运气不错,是向上的shake:)。这里总结一下销量预测的基本流程。
从下图可以看到这次比赛排名波动非常大。
正文
读取数据
首先导入我们常用的包,因为是时间序列所以导入datetime包处理时序;因为kaggle的数据集非常大,笔者的内存实在顶不住,需要gc.collect()帮助垃圾回收;numpy、pandas、sklearn自不必说;threading和queue用多线程加速特征工程;笔者常用的融合模型,万金油的lightgbm框架。
from datetime import date, timedelta,datetime
import gc
import pandas as pd
import numpy as np
from sklearn.preprocessing import LabelEncoder
import lightgbm as lgb
from threading import Thread
from queue import Queue
然后我们加上一些宏变量,便于后续的特征处理和日后的代码移植。
date_last = 1913 # 最大天数
max_lags = 1000 # 从文件读取的天数(省内存)
fday = datetime(2016, 4, 25) # 训练集最大日期
start_day = date_last - max_lags # 开始天数
p_horizon = 28 # 预测天数
TARGET = 'sales' # 目标值
导入数据,因为数据特征会分成文件,需要使用pd.merge()进行合并,不要想着就是一个train.csv文件,这是不可能的,因为时序特征会用到rolling(窗口滑动)的方法,需要前面的销量作为后续预测的特征。这里就拿price和sales文件举个例子。
这里解释一下为什么要转化为categories类型,因为例如item_id这种特征属于字符串类型,这是模型无法接受的数据,所以我们需要将其进行编码,比如one-hot,但是种类实在太多了,ont-hot会造成数据稀疏(且占据大量内存),而且树模型可以自行划分特征,所以直接利用categories类型转化为数字编码了,和LabelEncoder其实是一样的,用哪个都行。
为sales增加28个Nan,这是要预测的天数。
# 数据类型
numcols = [f"d_{day}" for day in range(start_day,date_last+1)]
catcols = ['id', 'item_id', 'dept_id','store_id', 'cat_id', 'state_id']
dtype = {numcol:"float32" for numcol in numcols}
dtype.update({col: "category" for col in catcols if col != "id"})
price_dtypes = {'store_id': 'category', 'item_id': 'category', 'wm_yr_wk': 'int16','sell_price':'float32' }
# 将price编码
price=pd.read_csv('数据/sell_prices.csv', dtype = price_dtypes)
for col,col_dtype in price_dtypes.items():
if col_dtype == 'category':
price[col] = price[col].cat.codes.astype('int16')
price[col] -= price[col].min()
# 将train_sale编码
train_sale=pd.read_csv('数据/sales_train_validation.csv',usecols = catcols + numcols, dtype = dtype)
for col in catcols:
if col != "id":
train_sale[col] = train_sale[col].cat.codes.astype('int16')
train_sale[col] -= train_sale[col].min()
# 增加要预测的天数
for day in range(date_last+1, date_last+p_horizon+1):
train_sale[f"d_{day}"] = np.nan
# 将 train_sale和price合并
df = df.merge(price,on=["store_id", "item_id", "wm_yr_wk"],copy=False)
上面不看也问题不大,总而言之,我们得到了以下格式的数据,sales那一列下面的Nan是需要预测的结果,这里用Nan占位。
将时间信息转化为模型可以学习的方式,比如2020/7/1日可以转化为年份特征year:2020,月份特征month:7,星期特征wday:3等。
date_features = {
"wday": "weekday",
"week": "weekofyear",
"month": "month",
"quarter": "quarter",
"year": "year",
"mday": "day",
}
for date_feat_name, date_feat_func in date_features.items():
if date_feat_name in df.columns:
df[date_feat_name] = df[date_feat_name].astype("int16")
else:
df[date_feat_name] = getattr(df["date"].dt, date_feat_func).astype("int16")
特征工程
然后进行我们的重头戏,特征工程,这一步非常重要,直接影响你的模型的预测精度,做机器学习基本一大半时间花在这个上面(当然如果你是神经网络deepNN爱好者,你可能也要花大量的时间在调参上面,所以我非常抗拒神经网络,能不用就不用)。常说特征工程决定你模型的上限,调参是为了接近那个上限。这里给一个万金油写法(将过去的数据作为未来的特征)。
很多人喜欢使用rolling和shift方法,因为那样写起来很简单,而且运行速度快,但是我个人认为不够直观,这里使用自己定义的函数get_timespan直接获取指定区间的数据,虽然麻烦了一点,但是更容易扩展,比如你想到什么更好的特征,可以直接在prepare_dataset上扩展。
def get_timespan(df, time, minus, periods, freq='D'):
return df[pd.date_range(time - timedelta(days=minus), periods=periods, freq=freq)]
# 创建特征
def prepare_dataset(df, time, name_prefix=None):
X = {}
# 这里创建你想要的特征,下面这里是举个例子
# for i in [3, 7, 14, 30, 60, 140]:
# tmp = get_timespan(df, time, i, i)
# X['diff_%s_mean' % i] = tmp.diff(axis=1).mean(axis=1).values
# X['mean_%s_decay' % i] = (tmp * np.power(0.9, np.arange(i)[::-1])).sum(axis=1).values
# X['mean_%s' % i] = tmp.mean(axis=1).values
# X['median_%s' % i] = tmp.median(axis=1).values
# X['min_%s' % i] = tmp.min(axis=1).values
# X['max_%s' % i] = tmp.max(axis=1).values
# X['std_%s' % i] = tmp.std(axis=1).values
X = pd.DataFrame(X, index = df.index) # 转换为dataframe
if name_prefix is not None:
X.columns = ['%s_%s' % (name_prefix, c) for c in X.columns]
if name_prefix == 'sales':
X['date'] = time # 加上时间信息
return X
因为速度感人,这里使用多线程加快处理速度(因为GIL锁的存在,没法实现真正的多线程,加速的效果不尽如意,但是windows下多进程的进程之间互相通信特别麻烦,我懒写)
# 准备数据
print("Preparing dataset...")
num_days = 800 # 需要创建特征的天数
t2016 = fday - timedelta(num_days)
dt = []
thread_list = []
q_day = Queue()
for i in range(num_days):# 加入队列
delta = timedelta(days=i)
q_day.put(t2016 + delta)
def func():
while q_day.empty() == 0:
day = q_day.get()
X_tmp = prepare_dataset(df_sales, day, name_prefix = 'sales') # 创建特征
X_tmp = df.merge(X_tmp, on=['date','id'], copy=False) # 特征与训练集合并
dt.append(X_tmp)
print(day)
gc.collect()
for i in range(4):# 创建线程
t = Thread(target=func)
thread_list.append(t)
for t in thread_list:# 启动线程
time.sleep(3) # 线程之间错开运行
t.start()
for t in thread_list:# 等待线程结束
t.join()
dt = pd.concat(dt, axis=0)
其中函数所需要的数据格式df_sales和创建方法是这样的
# 每种商品每天的销量
df_sales = df.loc[:,['id','date','sales']]
df_sales = df_sales.set_index(['id','date'])['sales'].unstack()
df_sales = df_sales.reset_index()
df_sales = df_sales.set_index('id')
准备训练集和验证集
分出训练集和测试集,因为时序问题的特殊性,应该用过去的数据预测未来,用未来的数据预测过去是不科学的,这也造成了我们做local cv的苦难,不能将数据集一分为四进行验证模型。这里我用了最近的5天数据作为验证集。
# 训练集和验证集
x_train = dt[dt.date < (fday - timedelta(5))]
x_test = dt[(dt.date < fday) & (dt.date > (fday - timedelta(5)))] # 5天的数据
y_train = x_train[TARGET] # 训练集target
y_test = x_test[TARGET] # 测试集target
print('训练集数量:'+str(x_train.shape))
建模预测
终于到了我们激动人心的模型训练环节(谁不想做一个优雅的调参侠呢,特征工程这种累活让别人干最好233)
首先设置模型参数,这里根据自己的情况设定
# lightgbm建模预测
params = {
'boosting_type': 'gbdt',
'objective': 'tweedie',
'tweedie_variance_power': 1.1,
'metric': 'rmse',
'subsample': 0.5,
'subsample_freq': 1,
'learning_rate': 0.03,
'num_leaves': 2**11-1,
'min_data_in_leaf': 2**12-1,
'feature_fraction': 0.5,
'max_bin': 100,
'n_estimators': 1400,
'boost_from_average': False,
'verbose': -1,
# 'device':'gpu',
}
训练模型,这里采用lightgbm的原生接口。
# 转换数据集格式
lgb_train = lgb.Dataset(x_train, label = y_train)
lgb_test = lgb.Dataset(x_test, label = y_test)
# 训练
m_lgb = lgb.train(params, lgb_train, valid_sets = [lgb_test],verbose_eval=100)
m_lgb.save_model("model.lgb")
前面提到我们需要采用滚动窗口预测的方式将过去的数据作为现在的特征,我们只能一天一天预测,将前一天预测的结果作为后一天的特征
# 预测
for days in range(p_horizon):
current_day = fday + timedelta(days=days)
X_tmp = prepare_dataset(df_sales, current_day, name_prefix = 'sales')
# X_tmp2 = prepare_dataset(df_price, current_day, name_prefix = 'price')
# X_tmp = pd.concat([X_tmp, X_tmp2], axis=1)
x_predict = df.merge(X_tmp, on=['date','id'], copy=False)
x_predict = x_predict
y_predict = m_lgb.predict(x_predict,axis=1))
print(current_day) # 打印日期
大功告成!这么一系列操作你应该可以得到一个还算可以的结果,当然争奖牌是远远不够的。后续我们还需要进行模型验证、参数优化等一系列操作,这里就不展开讲了:)
tip
因为数据集一般都非常大,本次比赛的数据全部展开要达到千万级的数据量,占用非常大的内存,在服务器上直接爆内存,只能用家里的电脑借助虚拟内存撑一撑。这里附一个降低精度的代码,可以有效缓解内存不够的窘境。
import pandas as pd
import numpy as np
def reduce_mem_usage(df):
start_mem = df.memory_usage().sum() / 1024**2
for col in df.columns:
col_type = df[col].dtype
if str(col_type)[:3] == 'int' or str(col_type)[:5] == 'float':
c_min = df[col].min()
c_max = df[col].max()
if str(col_type)[:3] == 'int':
if c_min > np.iinfo(np.int8).min and c_max < np.iinfo(np.int8).max:
df[col] = df[col].astype(np.int8)
elif c_min > np.iinfo(np.int16).min and c_max < np.iinfo(np.int16).max:
df[col] = df[col].astype(np.int16)
elif c_min > np.iinfo(np.int32).min and c_max < np.iinfo(np.int32).max:
df[col] = df[col].astype(np.int32)
elif c_min > np.iinfo(np.int64).min and c_max < np.iinfo(np.int64).max:
df[col] = df[col].astype(np.int64)
else:
# pandas没有float16类型, 可能会报错
if c_min > np.finfo(np.float16).min and c_max < np.finfo(np.float16).max:
df[col] = df[col].astype(np.float16)
elif c_min > np.finfo(np.float32).min and c_max < np.finfo(np.float32).max:
df[col] = df[col].astype(np.float32)
else:
df[col] = df[col].astype(np.float64)
elif col_type == 'object':
df[col] = df[col].astype('category')
end_mem = df.memory_usage().sum() / 1024**2
if verbose: print('Mem. usage decreased to {:5.2f} Mb ({:.1f}% reduction)'.format(end_mem, 100 * (start_mem - end_mem) / start_mem))
return df
后言
等我有空整理一下特征工程的代码来讲一下我这次的特征工程是怎么做的,当然这个有空就不知道是啥时候了233