PySpark.ml时间序列特征工程

本文为销量预测第4篇：时间序列特征工程
第1篇：PySpark与DataFrame简介
第2篇：PySpark时间序列数据统计描述，分布特性与内部特性
第3篇：缺失值填充与异常值处理
第5篇：特征选择
第6篇：简单预测模型
第7篇：线性回归与广义线性模型
第8篇：机器学习调参方法
第9篇：销量预测建模中常用的损失函数与模型评估指标
特征工程是将原始数据转化为有用的特征，更好的表示待处理的实际问题，提升数据对预测任务准确性。

未经处理的特征可能有以下问题：

• 量纲不同：量纲不一致，在线性回归等模型中，不同的量纲通过权重系数无法比较其重要性,同时量纲不同也不利于模型训练。可通过标准、归一化来解决;
• 信息冗余：对于某些数值特征，其包含的有效信息转化为区间划分更加合理，例如气温，28度与27.4度，连续的数值并没有那么重要，使用分桶转换表示更合理。同时也可以通过二值化加工"是否为月末"这样的节假日属性;
• 类别特征不能直接纳入模型：某些机器学习算法和模型只接受数值类型的特征输入，此时就需要将类别特征转换为数值特征。通常使用one_hot的方式来转换，从而获得非线性的效果;
• 信息利用率低：不同的机器学习算法和模型对数据中信息的利用是不同，可以利用多项式生成高阶非线性特征。
• 特征过多：在数据集小或者复杂模型上，过多的特征会导致模型过拟合，在特征预处理阶段可以通过降维的方式，提取部分有价值的特征，比如使用主成分分析(PCA)。
  所以下面先展开讲解借助SPARK进行特征预处理的常见方法。

1.特征预处理

1).二值化与分桶

二值化可以将数值型（numerical）的特征经过指定阈值threshold得到布尔型(boolean)数据，将大于阈值的赋值为1，对于数据分布为Bernoulli时的概率估计来说有用。

$$x^{\prime }=\{\begin{array}{l}1,x>\text{ threshold }\\ 0,x\le \text{ threshold }\end{array}$$

from pyspark.ml.feature import Binarizer
#二值化放入的inputCol特征数据类型必须是double类型
binarizer=Binarizer(threshold=10.0,inputCol='feature_v1',outputCol='binarized_feature_v1'
binarizedDataFrame=binarizer.transform(df.select('feature_v1'))

对连续变量进行分桶(Bucketizer)或分位数分桶(QuantileDiscretizer)有以下好处：

1.用粗粒度描述特征，减少过拟合的风险

2.增加稀疏数据的概率，减少计算量

3.减少噪声数据的影响，提升模型的鲁棒性

4.离散后特征便于计算交叉特征，进入非线性，提升表达能力。

分桶代码示例如下：

from pyspark.ml.feature import Bucketizer
#给定边界分桶离散化边界
splits=[-float('inf'),-0.5,0.0,0.5,float('inf')]
bucketizer=Bucketizer(splits=splits,inputCol='feature_v1',outputCol='bucketed_feature_v1')
bucketedData=bucketizer.transform(df.select('feature_v1'))

分位数离散化代码示例如下：

#按分位数分桶离散化——分位数离散化
from pyspark.ml.feature import QuantileDiscretizer

discretizer=QuantileDiscretizer(numBuckets=4,inputCol='feature_v1',outputCol='quantile_feature_v1')   #numBuckets指定分桶数
result=discretizer.fit(df.select('feature_v1')).transform(df.select('feature_v1'))

2).最小最大值标准化(MinMaxScaler)

区间缩放，返回值为缩放到[0, 1]区间的数据,当有新数据加入时，由于max和min的变化，可能需要重新定义；同时MinMaxScaler对异常值敏感。

$$x^{\prime }=\frac{x-min}{max-min}$$

from pyspark.ml.feature import MinMaxScaler
df = spark.createDataFrame([(Vectors.dense([-2.0, 2.3]),),
                      (Vectors.dense([0.0, 0.0]),),
                      (Vectors.dense([0.6, -1.1]),)],
                     ["features"])

min_max_scaler= MinMaxScaler(inputCol='features', outputCol='min_max_norm')
min_max_fit=min_max_scaler.fit(df)
min_max_result=min_max_fit.transform(df)

3).绝对值归一化MaxAbsScaler

在原始数据的基础上除以最大值的绝对数，将属性缩放到[-1,1]，不会破坏数据原本稀疏性。

$$x^{\prime }=\frac{x}{\left|x_{\max }\right|}$$

from pyspark.ml.feature import MaxAbsScaler
df = spark.createDataFrame([(Vectors.dense([-2.0, 2.3]),),
                      (Vectors.dense([0.0, 0.0]),),
                      (Vectors.dense([0.6, -1.1]),)],
                     ["features"])

max_abs_scaler= MaxAbsScaler(inputCol='features', outputCol='max_abs_norm')
max_abs_fit=max_abs_scaler.fit(df)
max_abs_result=max_abs_fit.transform(df)

4).特征标准化StandardScaler

标准化的前提是特征服从正态分布，标准化之后数据分布为标准正态分布，标准化消除了数据原本的实际意义。
$$x^{\ast }=\frac{x-\mu }{\sigma }$$

from pyspark.ml.feature import StandardScaler 
scaler = StandardScaler(inputCol="inputs", outputCol="scaled_features")
scaler_fit = scaler.fit(df)
scaled_result = scaler_fit.transform(df)

5).Normalizer （正则化）

Spark中的Normalizer的作用范围是每一行，使每一个行向量的范数变换为一个单位范数，Normalization是对每个样本计算其p-范数，对该样本中每个元素除以该范数，将原始特征Normalizer后,可使得机器学习算法有更好的表现。

$$\begin{gathered} 1范数(L1)：║x║1=│x1│+│x2│+\dots +│xn│ \\ 2范数(L2)：║x{║}^2=\sqrt{(x_1²+x_2²+\dots +x_n²)} \\ \infty 范数：║x║\infty =max(│x_1│，│x_2│，\dots ，│x_n│) \end{gathered}$$

$$pNorm:x^{\ast }=\frac{x}{pNnorm}$$

from pyspark.ml.feature import Normalizer
df = spark.createDataFrame([(Vectors.dense([-2.0, 2.3]),),
                      (Vectors.dense([0.0, 0.0]),),
                      (Vectors.dense([0.6, -1.1]),)],
                     ["features"])
normalizer = Normalizer(inputCol="features", outputCol="normFeatures", p=2.0)
l2NormData = normalizer.transform(df)

6).多项式特征(PolynomialExpansion)

以特征向量(x1,x2)为例，如果degree =2，输出为:

$$\begin{gathered} input=\left(X_1,X_2\right) \\ output=\left(X_1,X_2,X_1^2,X_1{X}_2,X_2^2\right) \end{gathered}$$

提示Spark中的多项式特征没有0次幂项,sklearn.preprocessing.PolynomialFeatures中有参数include_bias,默认为 True 。如果为 True 那么结果中就会有 0 次幂项，即全为1这一列。
$$\left(1,X_1,X_2,X_1^2,X_1{X}_2,X_2^2\right)$$
多项式特征不仅能够能在原特征的基础上形成更高次项，也会生成交互项，获得非线性关系,在带来更强的数据表达能力的同时，也需防止阶数太高可能产生的过拟合问题。关于生成多项式特征，可以在Spark.SQL中手动对多列进行乘积运算。

from pyspark.ml.feature import PolynomialExpansion
from pyspark.ml.linalg import Vectors

df = spark.createDataFrame([(Vectors.dense([-2.0, 2.3]),),
                      (Vectors.dense([0.0, 0.0]),),
                      (Vectors.dense([0.6, -1.1]),)],
                     ["features"])

ploy_df = PolynomialExpansion(degree=3, inputCol="features", outputCol="poly_features")
poly_features = ploy_df.transform(df)

7).独热编码OneHotEncoder

就是把数据变成(1,0,0,…,0),(0,1,0,0,…,0)，该特征属性有多少类别就有多少维。Spark中在处理OneHot之前一般先要转换成字符串索引(StringIndexer)，将字符串列编码为标签索引列,再做OneHot处理，示例如下:

df = df.withColumn('dayofweek', dayofweek('dt'))
df = df.withColumn("dayofweek", df["dayofweek"].cast(StringType()))
dayofweek_ind = StringIndexer(inputCol='dayofweek', outputCol='dayofweek_index')
dayofweek_ind_model = dayofweek_ind.fit(df)
dayofweek_ind_ = dayofweek_ind_model.transform(df)
onehotencoder = OneHotEncoder(inputCol='dayofweek_index', outputCol='dayofweek_Vec')
df = onehotencoder.transform(dayofweek_ind_)

8).降维 PCA（主成分分析 ）

维数灾难是机器学习中常见的现象，随着特征数增加，需要处理的数据相对于特征形成的空间而言比较稀疏，由有限训练数据拟合的模型可以很好的适用于训练数据，但对于未知的测试数据，很大几率距离模型空间较远，训练的模型不能处理这些新的未知数据点，从而形成“过拟合”的现象。在特征预处理阶段，可以通过降维的方式减轻维度灾难，常用的方法有主成分分析(PCA)。比如销售量，销售额，进店客流等属于高度相关的特征，针对数据集较小或者模型复杂度高时，如需使用全部特征，且为避免过拟合，此时就可以选择降维手段。

PCA主要包含以下几个步骤：
　　1、标准化样本矩阵中的原始数据；
　　2、获取标准化矩阵的协方差矩阵；　　
　　3、计算协方差矩阵的特征值和特征向量；
　　4、依照特征值的大小，挑选主要的特征向量；
　　5、生成指定维度的新特征。

#从hive中读取最新的特征列
def read_importance_feature():
    """
    :return: list of importance of feature
    """
    importance_feature = spark.sql("""select feature from temp.selection_result where cum_sum<0.99 and update_date 
    in (select max(update_date) as update_date from app.selection_result)""").select("feature").collect()
    importance_list = [row.feature for row in importance_feature]
    return importance_list

inputCols=read_importance_feature()
#读取数据
df=spark.sql("""select * from temp.dataset_feature'""")
df = df.na.fill(0)

#先把特征转换为向量

feature_vector = VectorAssembler(inputCols=inputCols, outputCol="original_features")
output = feature_vector.transform(df)
features_label = output.select("shop_number", "item_number", "dt", "original_features", "label")

#放入向量
pca = PCA(k=7,inputCol="original_features",outputCol="features")
model = pca.fit(features_label)
pca_result = model.transform(features_label).select("shop_number", "item_number", "dt","features", "label")

Spark中还有其他的特征预处理方式，如关于文本的StopWordsRemover、分词Tokenizer,正则匹配取词RegexTokenizer,TF-IDF词编码等,因与销量预测任务相关度降低，此处也就略去不表，感兴趣的读者可查询其他相关材料。

完成以上特征预处理以后，下面讲解在销量预测中最常用的特征工程。

2.日期特征

日期特征是时序中较为重要的一类特征，可以基于此计算得到序列关于日期的季节性规律。

1).日期拆解类

把带有日期的数据拆解到不同的日期粒度，比如‘2021-01-02’，可以得到年，月，日，季度等基础特征，同时Spark.SQL支持以下方式：

特征名称	Spark.SQL
年份	year
季度	quarter
月份	month
日	day
分	minute
一年中的第n周	weekofyear
星期几	dayofweek
月中第几天	dayofmonth

2).日期判断类

在基础的日期信息上，还可以进一步加工，比如，对月中第几天，可以加工为是否月初和是否月末等信息：

以下列举常用的日期特征衍生：

• 是否月初
• 是否月末
• 是否季节初
• 是否季节末
• 是否年初
• 是否年尾
• 是否周末
• 是否为节假日
• 是否营业时间

星期的one_hot编码以及手动生成二值化特征"是否月末"等特征衍生方式方式可以参考如下代码：

from pyspark.sql.functions import *
from pyspark.sql.types import *
from pyspark.ml.feature import OneHotEncoder
from pyspark.ml.feature import StringIndexer

df=df.withColumn('year',year('dt'))
df=df.withColumn('quarter',quarter('dt'))
df=df.withColumn('month',month('dt'))
df=df.withColumn('day',dayofmonth('dt'))
df=df.withColumn('dayofweek',dayofweek('dt'))
df=df.withColumn('weekofyear',weekofyear('dt'))

#是否月末编码,cast
df = df.withColumn('day', df["day"].cast(StringType()))

df = df.withColumn('month_end',when(df['day'] <=25,0).otherwise(1))

#星期编码--将星期转化为0-1变量
dayofweek_ind = StringIndexer(inputCol='dayofweek', outputCol='dayofweek_index')
dayofweek_ind_model = dayofweek_ind.fit(df)

dayofweek_ind_ = dayofweek_ind_model.transform(df)
onehotencoder = OneHotEncoder(inputCol='dayofweek_index', outputCol='dayofweek_Vec')
df = onehotencoder.transform(dayofweek_ind_)

需特意阐述一点，星期几虽然是整数类型，可以直接纳入机器学习模型中做训练，但是[1,2,3,4,5,6,7]的取值中，如，4/1=4，并不能说星期四是周一的4倍，不能说明周四在销售数据上从时间上看比周一大，也就是此时的星期几数据虽然在数据类型上可以是整数，但是其意义不具备连续型数据的可比较与可加性，所以需要作为类别变量做特殊处理，所以特征加工还是应该遵循常识和逻辑，不得无脑把加工的特征直接丢进模型中，否则会训练出错误模型，导致上线预测效果不稳定或者非常差。

3).节假日处理

关于日期的处理，还有一类比较特殊在时序领域需特意关注的是节假日信息，此处单独拿出来讲解，该部分内容参照Prophet库中对节假日的处理方式，即手工维护一个节假日表，包含节假日名称，日期，前后受节假日影响的天数，以下以儿童节和”618“为例。

import pandas as pd
children_day = pd.DataFrame({
        'holiday': 'children_day',
        'ds': pd.to_datetime(['2019-06-01', '2020-06-01']),
        'lower_window': -1,
        'upper_window': 0,})
shopping_618 = pd.DataFrame({
        'holiday': 'shopping_618',
        'ds': pd.to_datetime(['2019-06-18', '2020-06-18']),
        'lower_window': 0,
        'upper_window': 1,})

holidays_df = pd.concat((children_day,shopping_618))
holidays_set = holidays_df[['ds','holiday','lower_window','upper_window']].reset_index()

以上通过spark.sql内置的函数对日期进行拆解，同时使用pyspark中的ml.feature模块处理，one_hot和特征的类型转换，也因此展示了spark.sql的灵活和spark中机器学习模型对于数据特征处理的强大，后面也会介绍另一个特征加工利器Spark.UDF函数，用以生成更加复杂的特征。

3.统计特征

1).滞后（lag）特征

with lag_sale as 
(
select store_id,sku_id,sale_date,sale_qty,
lag(sale_qty,1) over(partition by store_id,sku_id order by sale_date) as lag1qty,
lag(sale_qty,2) over(partition by store_id,sku_id order by sale_date) as lag2qty,
lag(sale_qty,3) over(partition by store_id,sku_id order by sale_date) as lag3qty,
lag(sale_qty,4) over(partition by store_id,sku_id order by sale_date) as lag4qty,
lag(sale_qty,5) over(partition by store_id,sku_id order by sale_date) as lag5qty,
lag(sale_qty,6) over(partition by store_id,sku_id order by sale_date) as lag6qty,
lag(sale_qty,7) over(partition by store_id,sku_id order by sale_date) as lag7qty,
lag(sale_qty,14) over(partition by store_id,sku_id order by sale_date) as lag14qty,
lag(sale_qty,21) over(partition by store_id,sku_id order by sale_date) as lag21qty,
lag(sale_qty,28) over(partition by store_id,sku_id order by sale_date) as lag28qty,
lag(sale_qty,35) over(partition by store_id,sku_id order by sale_date) as lag35qty,
from  dataset_fix_with_future
)
select
a.store_id,
a.sku_id,
a.sale_date,
a.sale_qty,
nvl(b.lag1qty,0) lag1qty,
nvl(b.lag2qty,0) lag2qty,
nvl(b.lag3qty,0) lag3qty,
nvl(b.lag4qty,0) lag4qty,
nvl(b.lag5qty,0) lag5qty,
nvl(b.lag6qty,0) lag6qty,
nvl(b.lag7qty,0) lag7qty,
nvl(b.lag14qty,0) lag14qty,
nvl(b.lag21qty,0) lag21qty,
nvl(b.lag28qty,0) lag28qty,
nvl(b.lag35qty,0) lag35qty,
nvl(b.lag7qty/b.lag14qty,1) as qty_slope,
nvl(b.lag7qty-b.lag14qty,0) as qty_diff
from temp.dataset_future a
left join lag_sale b
on a.store_id=b.store_id and a.sku_id=b.sku_id and a.sale_date=b.sale_date

以上代码生成的特征有：

• 使用窗口函数lag对生成滞后特征；

• 其中qty_slope为最近两个周期的比例；

• 同时把二者相减生成增长特征qty_diff；

• nvl函数对null值进行填补为0。

2). 滑窗统计特征

滑窗统计特征是机器学习算法处理时序问题最经典的处理方式之一，通常情况下都是最重要的特征类。窗口大小不宜过大或者过小，通常去到序列中半个或者一个周期为佳，比如对包含多个年份的数据时间序列中，滑窗以3个时间点（月份），对于处理天这个粒度上的序列数据，如果存在以星期为周期的序列上，则取7作为窗口大小，如果窗口太小，则对于序列的波动太敏感，针对这样类似于这样的"超参数",可以结合业务背景和时间序列理论和作图分析进行人为设定，如果对待分析建模的数据没有相关背景支撑，则借助机器学习对超参数的确定方式，设置若干个可能的取值，使用模型训练效果最好的参数取值，同时，在处理序列较长或者存在多种周期季节模式的序列时，也可以使用多种不同大小的窗口函数，比如，针对存在180天的序列，除了使用7天的滑窗，也可以同时取30天的窗口。如下图7.7。

with lag_windows_df as (
SELECT 
store_id,
sku_id,
sale_date,
sale_qty,
avg(lag1qty) over(partions BY store_id,sku_id order by sale_date rows between 6 preceding and current row) as lag1_7_avg,
max(lag1qty) over(partions BY store_id,sku_id order by sale_date rows between 6 preceding and current row) as lag1_7_max,
min(lag1qty) over(partions BY store_id,sku_id order by sale_date rows between 6 preceding and current row) as lag1_7_min,
stddev_samp(lag1qty) over(partions BY store_id,sku_id order by sale_date rows between 6 preceding and current row) lag1_7_std,
skewness(lag1qty) over(partions BY store_id,sku_id order by sale_date rows between 6 preceding and current row) as lag1_7_skew,
kurtosis(lag1qty) over(partions BY store_id,sku_id order by sale_date rows between 6 preceding and current row) as lag1_7_kurt
from temp.dataset_future)
select
store_id,
sku_id,
sale_date,
sale_qty,
lag1_7_avg,
lag1_7_max,
lag1_7_min,
lag1_7_std,
lag1_7_skew,
lag1_7_kurt,
nvl(lag1_7_std/lag1_7_avg,1) as cv_1_7
from lag_windows_df

使用over partition by窗口函数，统计窗口期内的AVG,STD,MAX等指标。
特征工程是一个需长期持久化完善的建模任务之一，其重要性怎么强调都不过分，也是日常工作花费时间最多的地方，需要结合业务发挥创造性。以上所讲解的方法和处理方式只是其中一部分，限于使用SPARK这一工具与篇幅，同时考虑内容的适普性，只书写了以上内容。