【Python】沃尔玛销售数据挖掘与可视化分析实战

M5比赛是最近的一次，目标：我们将要预测零售巨头沃尔玛在未来28天提供的销售数据。本文是该任务的一部分：探索性数据分析，后面我们继续还会介绍特征工程、模型训练和预测部分。

我们正在处理42,840个分层时间序列。数据是从美国3个加利福尼亚州(CA)， 德克萨斯州(TX)和威斯康星州(WI)，7个部门的3049种单独产品中获得的。这里的“分层”表示可以在不同级别上汇总数据：商店级别，部门级别，产品类别级别和州级别。销售信息可以追溯到2011年1月至2016年6月。除了销售数量，我们还提供了有关价格，促销和节假日的相应数据。

请注意，我们已经被警告大多数时间序列都包含零值。

kaggle competitions download -c m5-forecasting-accuracy 或者登录 https://www.kaggle.com/c/m5-forecasting-accuracy/data进行下载 或者@公众号：数据STUDIO 后台联系云朵君获取。

使用的工具

import pandas as pd
import numpy as np

import seaborn as sns
from matplotlib import colors
import matplotlib.pyplot as plt
from matplotlib.colors import Normalize

import plotly.express as px
import plotly.graph_objects as go
import plotly.figure_factory as ff
from plotly.subplots import make_subplots
from tqdm import tqdm
import warnings
warnings.filterwarnings('ignore')

特别提出，本次分析使用的是交互式可视化库 plotly。我们之前有详细介绍过该库的使用，需要学习的朋友可以移步:
👉 深入了解 Plotly 高级技术，附实用代码示例
👉 Python Plotly交互可视化详解
👉 当Sklearn遇上Plotly，会擦出怎样的火花？

• Plotly是一个非常著名且强大的开源数据可视化框架，它通过构建基于浏览器显示的web形式 的可交互图表来展示信息。
• Plotly 是基于 HTML 显示的，所以这里推荐使用 Jupyter notebook作为开发工具。
• Plotly_exprress 则是对 Plotly 的高级封装，上手容易，它对 Plotly 的常用绘图函数进行了 封装。缺点是没有 plotly 那样自由度高，类似 Seaborn 和 Matplotlib 的关系。

探索性数据分析（EDA）

所有数据分析和探索并且包含数据清洗，我们熟知的:

• 数据First Look(类型，每列的意义，目标值)
• 数据检查和清洗(主键，空值，异常值)
• 趋势分析(通过作图的手段进行目标值分析，特征值分析)

原始数据分析

读入数据

• calendar

日历数据：包含有关产品销售日期的信息

• date：
• wm_yr_wk：1 + 年 + 第几周
• weekday：周
• wday：日
• month：月
• year：年
• d：关联主键
• event_name_1：事件名称1
• event_type_1：事件类型1
• event_name_2：事件名称2
• event_type_2：事件类型2
• snap_CA：二进制值，是否是优惠日CA
• snap_TX：二进制值，是否是优惠日TX
• snap_WI：二进制值，是否是优惠日WI

分析发现：

• date 不是datetime类型 calendar['date'] = pd.to_datetime(calendar['date']) ，如果downcast之后会被改动过来
• event features 包含很多 'NaN' 可以处理'NaN' ->'no_event'
• snap 是一个奇怪的字段，去了解SNAP的业务含义

美国联邦政府提供称为补充营养援助计划 (SNAP) 的营养援助福利。SNAP 为低收入家庭和个人提供电子福利转账借记卡以购买食品。在许多州，货币福利在每个月的挑选10天进行福利日，然后这10天里每一天，1/10 的人可以通过他们的卡获得福利。

因此该字段非常关键。作为为数不多的外部特征，后续EDA会针对SNAP专门进行分析

• sell_prices

销售表

这里有两个表。

• sales_train_evaluation

• sales_train_validation

量化处理

大致看一下数 据的大小、列名、内容、数据类型。 这里需要用： sell_prices.info() 查看数据情况。

'pandas.core.frame.DataFrame'>
RangeIndex: 6841121 entries, 0 to 6841120
Data columns (total 4 columns):
 #   Column      Dtype  
---  ------      -----  
 0   store_id    object 
 1   item_id     object 
 2   wm_yr_wk    int64  
 3   sell_price  float64
dtypes: float64(1), int64(1), object(2)
memory usage: 208.8+ MB

我们看到这里找个数据占了 208.8+ MB 的内存。通常会在项目中做着做着发现内存消耗太大了，一些表之间的merge操作很慢的时候回过头来做这个量化处理。

我们提供一个 downcast 函数：

sales = downcast(sales_train_validation)
price = downcast(sell_prices)
calendar = downcast(calendar)

处理后再次查看信息，字段 wm_yr_wk 从 int64 降到了 int16 ，字段 sell_price 从 float64 降到了 float32 ，而总体内存降到了 58.8 MB 。

RangeIndex: 6841121 entries, 0 to 6841120
Data columns (total 4 columns):
 #   Column      Dtype   
---  ------      -----   
 0   store_id    category
 1   item_id     category
 2   wm_yr_wk    int16   
 3   sell_price  float32 
dtypes: category(2), float32(1), int16(1)
memory usage: 58.8 MB

总览数据

首先根据各个 id 聚合，并对 item_id 统计，之后用plotly的 treemap 进行数据可视化结果。

group = sales.groupby(['state_id','store_id','cat_id','dept_id'],
                      as_index=False
                     )['item_id'].count().dropna()

从总体看整个数据item的分布。可以看到一共有 30490条数据 ，还可以看下其他细节，如其中加州(CA)数据占比是最大的。

由于代码较长，这里我们不做展示，完整代码获取：请在@公众号：数据STUDIO 后台回复 M5 获取。

抽样调查

为了后面的分析，我们需要先取出所有 id 和所有关键主键名(包含 d_ )。

ids = sorted(list(set(sales_train_validation['id'])))
d_cols = [c for c in sales_train_validation.columns if 'd_' in c]

探索每家店的销售情况

x_1 = sales_train_validation.loc[sales_train_validation['id'] == ids[2]].set_index('id')[d_cols].values[0]
x_2 = sales_train_validation.loc[sales_train_validation['id'] == ids[1]].set_index('id')[d_cols].values[0]
x_3 = sales_train_validation.loc[sales_train_validation['id'] == ids[17]].set_index('id')[d_cols].values[0]

接下来我们可视化查看结果。

一家门店一个单品销量

由于时间序列比较长，这种粗看是看不出更多的细节，因此我们需要更加细分。

因此我们对每家店进行分段查看。

ids = sorted(list(set(sales_train_validation['id'])))
d_cols = [c for c in sales_train_validation.columns if 'd_' in c]
x_1 = sales_train_validation.loc[sales_train_validation['id'] == ids[0]].set_index('id')[d_cols].values[0][:90]
x_2 = sales_train_validation.loc[sales_train_validation['id'] == ids[0]].set_index('id')[d_cols].values[0][350:450]
x_3 = sales_train_validation.loc[sales_train_validation['id'] == ids[0]].set_index('id')[d_cols].values[0][1300:1400]

可视化查看：

可以看出总体分布，并没有太大的波动，仔细分析得意得到如下：

• 销量中有很多0，面对这种销量在现实项目中，先去判断是不是能用业务理解去解释这些0。比如，在过去某便利店的销量预测项目中，发现扑克牌在按照天维度去预测的时候会有大量的0，一周可能只会在某几天卖出2-3副扑克牌，但是如果巧妙的将客户的需求从按照日维度预测未来1-28天，变成扑克牌品类按照周维度去预测，可能在周维度上每周扑克牌销量是稳定恒定值。
• 在没有业务的深度参与的项目中，我们还能用一些所谓的"去燥声"方法，去还原，平滑这些销量，并利用平滑完后的销量去做后续的特征工程的处理。
• 但作为一个实战项目，不会轻易去做去燥声处理。

横向对比不同门店的销量对比

首先取出每家门店不同商品的销售表，并通过关联主键 d_cols 与日历表的日期做一个 merge 。

past_sales = sales_train_validation.set_index('id')[d_cols].T.merge(calendar.set_index('d')['date'],
           left_index=True, right_index=True, validate='1:1').set_index('date')
past_sales.head()

然后对每一家门店做移动平均处理。

data = past_sales[store_items].sum(axis=1).rolling(90).mean()

我们发现，CA_3的销量是最高的，并且相对于第二名，差距还是拉的很大。此外，同在加州的CA_4门店，销量是最低的。

在周期性方面，所有的门店的销量趋势基本一致。因此这里我们可以考虑使用一些方法，将其中的周期性给提取出来。

我们知道，线图可以看出趋势，而箱形图不仅可以看出总体大小，还能看出其分布。

Tips:

这里注意到，我们使用了不同的绘图背景（主题），可以使用如下代码查看可以使用哪些主题模板。

import plotly.io as pio
list(pio.templates)

上图比较了数据集中每个商店的销售分布。加利福尼亚的商店似乎销售差异最大，这可能表明加利福尼亚的某些地方的增长速度明显快于其他地方，即存在发展差距。另一方面，威斯康星州和德克萨斯州的销售额似乎非常一致，没有太大差异。这表明这些州的发展可能更加一致。此外，加州商店的整体平均销售额似乎也最高。

当然，我们还可以使用直方图来展示结果。（目测这里使用直方图并没有箱形图形象，因此这里并不打算采用直方图）

后续也会介绍小提琴图，一种比箱图更牛的图示类型。

合并表

目前的数据集是每一行sales对应一条时间序列，这样一些外部变量，例如价格，节假日，甚至未来可能遇到的促销、人流量数据是无法直接加在上面的，所以我们需要将这种wide的宽表，变成长表，每一行代表的是某一家店某一个商品某一天的具体信息，包括销量，价格，节假日，促销等等。

下面是一个 pandas melt 函数，从换标合并转换为长表的例子。

首先看下sales表：

首先转换表：

df = pd.melt(sales, 
             id_vars=['id', 'item_id', 'dept_id', 'cat_id', 'store_id', 'state_id'], 
             var_name='d', 
             value_name='sold').dropna()

可看得出，把关联主键从原来的行，变成了现在的列 d ，其值放在列 sold 了。

然后再与日历表合并：

df = pd.merge(df, calendar, on='d', how='left')

最后再与销售价格表合并，最后得到总表：每一行代表的是某一家店某一个商品某一天的具体信息，包括销量，价格，节假日，促销等等。

df = pd.merge(df, sell_prices, 
              on=['store_id','item_id','wm_yr_wk'], how='left') 

Tips:

这里由于表格比较大，在表格转换时，比较耗时，因此我们做了一次操作后，最好把中间结果保持到磁盘。

如果你对parquet使用不熟悉，推荐阅读👉 详解 16 个 Pandas 读与写函数

这里推荐代码：

df.to_parquet('../data/raw_sales_data.parquet')

价格的分布分析

• 价格是按照周变化的, 首先分析门店维度的价格情况。
• 小提琴图(Violin Plot) 用于显示数据分布及其概率密度。这种图表结合了箱形图和密度图的特征，主要用来显示数据的分布形状。

下面按照每个地点、每个门店、每个商品聚合并求平均价格。

group_price_store = df.groupby(['state_id','store_id','item_id'],
                               as_index=False)['sell_price'].mean().dropna()

然后用小提琴图可视化出来。

• 很直观的，能观察到价格的分布在各大state，store似乎比较类似
• 然后可以观察一下最贵的商品和最便宜的商品是哪些

比如: HOBBIES_1_361 在CA的三家店都是最贵的，但是好像没有在 CA_1 的这家点销售，然后查看一下发现其实有销售，放大图可以看到这个SKU在 CA_1 是第二贵的；

HOUSEHOLD_1_060 在TX的几家店都是最贵的， HOBBIES_1_255 在 WI_1 和 WI_2 是最贵的；

总的来说，可以看到 HOBBIES&HOUSEHOLD 大类会贵一些， FOOD 比较便宜。

为了验证 FOOD 这个类别好像价格会便宜一些，可以再画一个按照 Store，cate 分类的图像。

之前有画过by store 销量随着时间变化的曲线，利用plotly.graph_objects可以画出互动性更强的图。

营业额分析

接下来分析一下营业额。营业额等于销量乘以单价。

df['revenue'] = df['sold']*df['sell_price'].astype(np.float32)

然后随机选取一家门店。

# 区分出SNAP/NON_SNAP的
store_sales = df[(df['state_id']=='CA')&(df['store_id']=='CA_1')&(df['date']<='2016-05-22')]
food_sales = store_sales[store_sales['cat_id']=='FOODS']
store_sales = store_sales.groupby(['date','snap_CA'],as_index=False)['sold','revenue'].sum()
snap_sales = store_sales[store_sales['snap_CA']==1]
non_snap_sales = store_sales[store_sales['snap_CA']==0]
food_sales = food_sales.groupby(['date','snap_CA'],as_index=False)['sold','revenue'].sum()
snap_foods = food_sales[food_sales['snap_CA']==1]
non_snap_foods = food_sales[food_sales['snap_CA']==0]
non_snap_foods

我们发现数据会遇到时间不连续的问题，因此需要对不连续的地方进行数据填充。

fill_dates(non_snap_sales)

将食品类单独提出分析。

这里定义了一个函数 plot_metric(df,state,store,metric) ，针对目标列进行画图分析，画出随着时间推移，某一个state，store的目标列变化图像。其中专门把 SNAP 时间段和 NON_SNAP 时间段分开来看。其中函数参数 metric 目标列名字，画图分析的目标列。

• 根据销量

fig = plot_metric(df,'CA','CA_1','sold')
fig.show()

• 根据销售额绘图，只需要跟换参数即可：

fig = plot_metric(df,'CA','CA_1','revenue')
fig.show()

特征工程

一个完整的特征工程可以包括但不限于:

1. 特征构建

• 销量预测：历史销量lag，历史销量的统计值(MA)，窗口统计值(Window max,min)
• 各类维度表信息的合并(商品，门店，价格，促销，节假日)

2. 特征选择

• 去除数据泄漏的特征
• 去除常数特征
• 去除相关性高的特征

3. 类别型特征的编码

对于时间序列问题，尤其是销量预测类型的表格类数据，特征工程大都有以下几类：

• 对于EDA发现的结果或者没有处理的部分进行处理（比如Nan值，在EDA过程中我们可能发现有Nan，在做特征第一步要把他们进行处理）
• lag features，历史数据，例如前一天的数据
• 针对lag features 的统计值，例如MA，一阶差分甚至二阶差分
• 各维度上的统计值，例如各item维度上的售卖量之和；然后对这些值还能做一些比值类的特征，比如城市A占总全国销售量的百分比。
• 将date列提取出时间相关的特征dayofweek, month, year, isMonthEnd, isMonthStart

对于特征工程里创建的特征，重点还是要思考它对于最后预测值的意义，一些比值，或者统计值是否能直接或者间接帮助到预测。

推荐特征工程相关文章👇
👉 一文带你用sklearn做特征工程
👉 特征工程：时间特征构造以及时间序列特征构造
👉 特征工程函数代码大全
👉 如何做好多元时间序列特征工程？

特征构建

首先处理缺失值。将 envent 变量中缺失值用 no_envent 填充。

cat=['event_name_1','event_type_1','event_name_2','event_type_2']
for i in cat:
    print(i)
    calendar[i] = calendar[i].cat.add_categories(["no_event"]).fillna('no_event')

ps：这里注意，由于前面已经设置了 cat 为分类变了，因此在用一个新分类值填充缺失值时，需要先新增一个类别值。否则该段代码将报错：

pandas ValueError: Cannot setitem on a Categorical with a new category, set the categories first

提取时间相关特征

通常做法是 date 字段转换为时间戳，再分别提取：

df['Year'] = df['Date'].dt.year
df['Year'] = df['Date'].dt.year
df['Month'] = df['Date'].dt.month
df['Day'] = df['Date'].dt.day
df['DayOfWeek'] = df['Date'].dt.dayofweek 
df['WeekOfYear'] = df['Date'].dt.weekofyear

本文中，还在这里也是提取了很多时间相关的特征：

• 增加 ' is_weekend ' 特征，代表是否是周末
• 增加 ' month_day '，代表这是一个月的第几天
• 增加 ' month_week_number '，代表这是这个月的第几周
• 增加 ' events_per_day '，代表一天有几个events

处理完成后，将所有的表合并：

sales=pd.melt(sales_train_evaluation,id_vars=['id','item_id','dept_id','cat_id','store_id','state_id'],var_name='d',value_name='demand')
sales=pd.merge(sales,calendar,on='d',how='left')
sales=pd.merge(sales,sell_prices,on=['item_id','store_id','wm_yr_wk'],how='left')

sales.info()

Int64Index: 58327370 entries, 0 to 58327369
Data columns (total 26 columns):
 #   Column             Dtype   
---  ------             -----   
 0   id                 category
 1   item_id            category
 2   dept_id            category
 3   cat_id             category
 4   store_id           category
 5   state_id           category
 6   d                  object  
 7   demand             int16   
 8   date               object  
 9   wm_yr_wk           int16   
 10  weekday            category
 11  wday               int8    
 12  month              int8    
 13  year               int16   
 14  event_name_1       category
 15  event_type_1       category
 16  event_name_2       category
 17  event_type_2       category
 18  snap_CA            int8    
 19  snap_TX            int8    
 20  snap_WI            int8    
 21  is_weekend         int8    
 22  month_day          int8    
 23  month_week_number  int8    
 24  events_per_day     int8    
 25  sell_price         float32 
dtypes: category(11), float32(1), int16(3), int8(9), object(2)
memory usage: 3.0+ GB

创建历史数据lag features

首先通过pandas的 shift 来我完成。

lags = [1,2,3,4,5]
for lag in lags:
    example_df['sales_lag_'+str(lag)] = example_df.sort_values(['date']).groupby(['store','sku'])['sales'].shift(lag)

计算lag feature的统计值

固定历史时间lag feature的平均值

比如 lag1-lag3 的平均值，代表过去三天的销量平均值，可以进一步引申到移动平均值 Moving Average 值。

常见的能想到两种方法：

1. 目前我们已经有lag1-lag3，再做平均值，如果要做移动平均值，对于每个lag值需要有一个不同的权重的话，将需要设置每个权重的list。
2. 用pandas自带的rolling函数，设定好一个窗口，然后计算滚动平均值。

# 方法1
ma_len = 3




    

ma_weights = [1 / ma_len for _ in range(ma_len)]
example_df['sales_MA_lag_1_3']= example_df[['sales_lag_1', 'sales_lag_2', 'sales_lag_3']].mul(ma_weights).sum(1)

# 方法2
example_df['rolling_mean_1_3'] = example_df.groupby(['store','sku'])['sales'].transform(lambda x: x.rolling(window=3, min_periods=1).mean().shift(1))

通过对比 sales_MA_lag_1_3 和 rolling_mean_1_3 列已经发现了两者有略微的不同。方法1中，存在缺失值时，将缺失值当作0处理，在物理意义上缺乏合理性。而在方法2中，通过控制参数 min_periods 设置最小区间观察值，即非缺失值个数。

累计历史时间lag feature的平均值

有固定窗口大小的 rolling window ，就会有非固定窗口大小的 expanding window 。

• rolling注重与随着时间变化，距离当天最近的7天的销量和（平均），忽略过去时间久的一些值
• expanding注重发生到当前位置的历史值都会在这个特征里面体现。

example_df['expanding_mean_1_3'] = example_df.groupby(['store','sku'])['sales'].transform(lambda x: x.expanding(min_periods=3).mean().shift(1))

当然，同理你还能做很多特征：

• 窗口的min max sum std
• 趋势类特征 e.g. lag2- lag1 lag3-lag2
• 比值类

特征选择

特征选择的通常有如下作用：

1. 简化模型，使模型更易于理解:去除不相关的特征会降低学习任务的难度。并且可解释性能对模型效果的稳定性有更多的把握
2. 改善性能:节省存储和计算开销
3. 改善通用性、降低过拟合风险:减轻维数灾难，特征的增多会大大增加模型 的搜索空间，大多数模型所需要的训练样本随着特征数量的增加而显著增加。特征的增加虽然能更好地拟合训练数据，但也可能增加方差

如何挑选特征：连续型特征和离散型特征

1. 手动选择不想要的特征，这类特征将导致数据泄漏，或者datetime类型无法被模型读入。
2. 挑选出常数的特征，这类特征在树模型中没有任何作用，且会影响速度。
3. 根据列中元素的类型，粗筛选出连续型特征和离散型特征。
4. 对数值型特征再进行一步计算，筛选出相关性极高的特征。

👉 特征选择技术总结
👉 最新特征筛选方法--Deep Lasso
👉 特征选择技术总结
👉 详解 5 大常用的特征选择方法！
👉 特征选择与特征提取最全总结
👉 特征选择与提取最全总结之过滤法

特征选择常用方法：

# 去除leak的特征，目标值，datetime值等不需要的特征
excluded_cols = ['Date', 'Sales','log_Sales', 'Customers','PromoInterval','monthStr'] 
init_cols = whole_df.columns
features = columns_minus(init_cols, excluded_cols)
# 去除常数项特征，如open全都是1
constant_cols = [col for col in features if whole_df[col].nunique() == 1]
features = columns_minus(features, constant_cols)
# 剩下特征区分numeric & categorical
num_features = whole_df[features].select_dtypes(include=[np.number]).columns.tolist() cate_features = columns_minus(features, num_features)
# num_features可以去计算correlation，默认是pearson correlation
corr_matrix = whole_df[num_features].corr().abs()
# 选择左上矩阵，目的为了观察哪些是high correlation的
upper = corr_matrix.where(np.triu(np.ones(corr_matrix.shape), k=1).astype(np.bool))
# 去除相关性高于某个threshold的特征
high_corr_cols = [col for col in num_features if any(upper[col] >= 0.9)]

例如我们本次筛选出相关性大于0.9的特征：

因此可以选择其中一列进行删除处理即可。

划分数据集

关于时间序列交叉验证，可以参考：

train=part_sales.loc[part_sales['day'