ODM和OLAP实现时序预测(三)
Part 3 时序预测之多步预测
这是系列文章中的最后一部分。这部分囊括了如何使用数据挖掘方法做多步
(
开循环
)
预测。正如
part1
中所讲,多步预测可以比单步预测做多次预测,不像单步预测那样,递归多步自然预测不能单独通过
SQL PREDICTION
函数计算。然而,我们可以使用带
MODEL
语句的
PREDICTION
函数来实现多步预测。下面,我将讲解如何使用这种方法。
方法
多步预测需要的方法与在
part1
和
part2
中所描述的方法相同。唯一的区别是我们如何记录模型。
预测
例如只包括头七天销售数据的销售时序
(ts_data)
,根据
part1
和
part2
中描述的方法创建回归数据挖掘模型
(ts_model)
,我们要使用
ts_model
模型预测接下来的七天
(8-14
天
)
的销售数据。正如我们在
part2
中所说明的那样,可以使用
PREDICTION
函数来进行单步预测,对于这个例子,下面的查询能够产生单步预测结果。
SELECT day, sales,
PREDICTION(ts_model USING LAG(sales,1) OVER (ORDER BY day) as l1, LAG(sales,2) OVER (ORDER BY day) as l2, LAG(sales,3) OVER (ORDER BY day) as l3, LAG(sales,4) OVER (ORDER BY day) as l4, LAG(sales,5) OVER (ORDER BY day) as l5, LAG(sales,6) OVER (ORDER BY day) as l6, LAG(sales,7) OVER (ORDER BY day) as l7 ) pred FROM ts_data a; |
因为模型要求早先七天的销售数据作为输入,最初我们只能获得第八天的销售数据,要预测
9-14
天的销售数据,一旦当天的销售数据可用就需要在每天结束运行这个查询来预测下一天的销售数据,尽管这也提供了些有用的信息,但并没有给我们留下更多的时间来做计划。理想情况,我们需要一次就生成接下来
7
天的销售预测。这就需要多步预测来帮忙。在多步预测过程中当有序列值时需要实际的序列值,没有真实的序列值时使用前一个预测值。
多步预测-未转换数据情况(No Transformation Case)
可以使用
SQL PREDICTION
函数和
MODEL
语句生成多步预测。例子如下:
SELECT d day, s sales, p pred
FROM ts_data a MODEL DIMENSION BY (day d) MEASURES (sales s, CAST(NULL AS NUMBER) p) RULES( p[FOR d FROM 1 TO 14 INCREMENT 1] = PREDICTION(ts_model USING NVL(s[CV()-1], p[CV()-1]) as l1, NVL(s[CV()-2], p[CV()-2]) as l2, NVL(s[CV()-3], p[CV()-3]) as l3, NVL(s[CV()-4], p[CV()-4]) as l4, NVL(s[CV()-5], p[CV()-5]) as l5, NVL(s[CV()-6], p[CV()-6]) as l6, NVL(s[CV()-7], p[CV()-7]) as l7) ) ORDER BY d |
MODEL
语句定义了一个多维数据表样式的数据视图,维的
index
是数据表中的行,度量是列。
RULE
语句可以为指定行范围的每个列定义计算规则。如果指定的行范围超过实际范围,会插入新行。对于
MODEL
语句的介绍可以参见
The SQL Model Query of Oracle 10g
白皮书和
"Announcing the New Model"
文档。
要使用
SQL MODEL
语句做多步预测,我们需要定义预测度量
(measure)
。度量使用
PREDICTION
函数计算并且历史时期范围和预测时期的范围一样。
我们分析一下上面的查询。查询描述了“数据表
(spreadsheet)
”视图数据。
DIMENSION
语句是指使用
day
列作为行的
index
。
MEASURES
语句定义了两个度量,原始的
sales
和
pred
列,
pred
列初始化为空,用来保存预测值。
RULES
语句说明如何计算
pred
列,
d(day
维度
)
从
1
到
14
,递增
1
,其中范围
8-14
并没有给出原始数据而且处于预测期。
Pred
列的规则适用
PREDICTION
函数计算预测值。
NVL
函数返回指定天的不为空的数据,如果为空则返回预测列。语法
s[CV()-1]
描述了早先天的
sales(s)
的销售值,如果
d=7
则
s[cv()-1]
返回的是
d=6
的销售数据。这种记号法提供了一种转换方式指定作为时序预测模型需要的迟滞
(lagged)
值。
通常,我们想计算一个时期的多步预测,当试图评估模型质量时,下面是一个例子。大多数情况下,我们使用下面的方法修改上面的查询:
SELECT d day, s sales, p pred FROM ts_data a MODEL DIMENSION BY (day d) MEASURES (sales s, CAST(NULL AS NUMBER) ap, CAST(NULL AS NUMBER) p) RULES( ap [FOR d FROM 1 TO 7 INCREMENT 1] = s[CV()], p[FOR d FROM 1 TO 14 INCREMENT 1] = PREDICTION(ts_model USING NVL(ap[CV()-1], p[CV()-1]) as l1, NVL(ap[CV()-2], p[CV()-2]) as l2, NVL(ap[CV()-3], p[CV()-3]) as l3, NVL(ap[CV()-4], p[CV()-4]) as l4, NVL(ap[CV()-5], p[CV()-5]) as l5, NVL(ap[CV()-6], p[CV()-6]) as l6, NVL(ap[CV()-7], p[CV()-7]) as l7) ) ORDER BY d |
上面的这个查询中,
MEASURES
语句定义了三个维度,原始的销售数据
sales
,辅助列
ap
,和列
p
初始为空保存测试数据。一旦预测被计算了,我们就可以通过比较预测值和实际值来评价模型。典型的评价包括:
Root Mean Square Error(RMSE)
,
Mean Absolute Error(MAE)
,
Mean Absolute Percentage Error(MAPE)
,以及最大的错误。
多步预测-发生数据转换情况(Transformation Case)
如果在模型构建前数据被转换了
(
见
part2
部分的例子
)
则转换需要包含在查询计算预测语句中。比如
part2
部分中使用了
LOG
转换,差分和常态
(z-scroe)
等方式处理数据,构建了
ts_model
回归模型。我们需要修改上面的查询:
SELECT d day, s sales, p pred FROM ts_data a MODEL DIMENSION BY (day d) MEASURES (sales s, CAST(NULL AS NUMBER) ts, CAST(NULL AS NUMBER) tp, CAST(NULL AS NUMBER) np, CAST(NULL AS NUMBER) dp, CAST(NULL AS NUMBER) p) RULES( ts[FOR d FROM 1 TO 14 INCREMENT 1] = (LOG(10,s[CV()]) - LOG(10,s[CV()-1]) + 0.010578734)/0.193587249 , tp[FOR d FROM 1 TO 14 INCREMENT 1] = PREDICTION(ts_model USING NVL(ts[CV()-1], tp[CV()-1]) as l1, NVL(ts[CV()-2], tp[CV()-2]) as l2, NVL(ts[CV()-3], tp[CV()-3]) as l3, NVL(ts[CV()-4], tp[CV()-4]) as l4, NVL(ts[CV()-5], tp[CV()-5]) as l5, NVL(ts[CV()-6], tp[CV()-6]) as l6, NVL(ts[CV()-7], tp[CV()-7]) as l7), np[FOR d FROM 1 TO 14 INCREMENT 1] = 0.193587249 * tp[CV()] - 0.010578734 , dp[FOR d FROM 1 TO 14 INCREMENT 1] = np[CV()] + NVL(LOG(10,s[CV()-1]),dp[CV()-1]) , p[FOR d FROM 1 TO 14 INCREMENT 1] = POWER(10,dp[CV()]) ) ORDER BY d |
DIMENSION
语句使用
day
列作为行的
index
。
MEASURE
语句定义了下面的度量:
l
Sales
是原始的销售数据
l
Ts
表示所有转换过后的数据,
LOG
转换、差分、
z-score
常态化
l
Tp
包含用于训练和预测时期的预测值,预测以转换数据相同的度量方式提供
l
Np
处理预测值的常态转换
l
Dp
处理差分赚
l
P
是与原始值相同方式的预测值,未转换的,销售序列,使用
POWER
函数来抵消
LOG(10)
的转换。
对于在原始序列做的每种转换我们都应该引入相反的转换步骤“
reverse transform
”,将预测值转换成与原始值相同的度量方式。上面的查询让我们见识到了
MODEL
的能力,
MODEL
语句能够以一种简洁的方式定义多重计算。
假设我们有预测期的实际值,我们可以引入一个新的度量
ap
,并且在所有的规则中使用
ap
替换
s
。修改后的查询如下:
SELECT d day, s sales, p pred FROM ts_data a MODEL DIMENSION BY (day d) MEASURES (sales s, CAST(NULL AS NUMBER) ap, CAST(NULL AS NUMBER) ts, CAST(NULL AS NUMBER) tp, CAST(NULL AS NUMBER) np, CAST(NULL AS NUMBER) dp, CAST(NULL AS NUMBER) p) RULES( ap [FOR d FROM 1 TO 7 INCREMENT 1] = s[CV()], ts[FOR d FROM 1 TO 14 INCREMENT 1] = (LOG(10,ap[CV()]) - LOG(10,ap[CV()-1]) + 0.010578734)/0.193587249, tp[FOR d FROM 1 TO 14 INCREMENT 1] = PREDICTION(ts_model USING NVL(ts[CV()-1], tp[CV()-1]) as l1, NVL(ts[CV()-2], tp[CV()-2]) as l2, NVL(ts[CV()-3], tp[CV()-3]) as l3, NVL(ts[CV()-4], tp[CV()-4]) as l4, NVL(ts[CV()-5], tp[CV()-5]) as l5, NVL(ts[CV()-6], tp[CV()-6]) as l6, NVL(ts[CV()-7], tp[CV()-7]) as l7), np[FOR d FROM 1 TO 14 INCREMENT 1] = 0.193587249 * tp[CV()] - 0.010578734, dp[FOR d FROM 1 TO 14 INCREMENT 1] = np[CV()] + NVL(LOG(10,ap[CV()-1]),dp[CV()-1]), p[FOR d FROM 1 TO 14 INCREMENT 1] = POWER(10,dp[CV()]) ) ORDER BY d |
示例1:航班乘客预测
我们使用
part2
中用到的
airline_SVM
模型来进行多步预测,与单步预测结果做比较。下面的查询生成所有的训练数据和
12
个月的预测数据
(
月份编号从
132-144)
。
SELECT m month, p passengers, pred FROM airline a MODEL DIMENSION BY (month m) MEASURES (a.passengers p, CAST(NULL AS NUMBER) ap, CAST(NULL AS NUMBER) tp, CAST(NULL AS NUMBER) tpred, CAST(NULL AS NUMBER) npred, CAST(NULL AS NUMBER) dpred, CAST(NULL AS NUMBER) pred) RULES( ap [FOR m FROM 1 TO 131 INCREMENT 1] = p[CV()], tp[FOR m FROM 1 TO 131 INCREMENT 1] = (LOG(10,ap[CV()]) - LOG(10,ap[CV()-1]) - 0.003919158)/0.046271162, tpred[FOR m FROM 1 TO 144 INCREMENT 1] = PREDICTION(airline_SVM USING NVL(tp[CV()-1],tpred[CV()-1]) as l1, NVL(tp[CV()-2],tpred[CV()-2]) as l2, NVL(tp[CV()-3],tpred[CV()-3]) as l3, NVL(tp[CV()-4],tpred[CV()-4]) as l4, NVL(tp[CV()-5],tpred[CV()-5]) as l5, NVL(tp[CV()-6],tpred[CV()-6]) as l6, NVL(tp[CV()-7],tpred[CV()-7]) as l7, NVL(tp[CV()-8],tpred[CV()-8]) as l8, NVL(tp[CV()-9],tpred[CV()-9]) as l9, NVL(tp[CV()-10],tpred[CV()-10]) as l10, NVL(tp[CV()-11],tpred[CV()-11]) as l11, NVL(tp[CV()-12],tpred[CV()-12]) as l12), npred[FOR m FROM 1 TO 144 INCREMENT 1] = 0.003919158 + 0.046271162 * tpred[CV()], dpred[FOR m FROM 1 TO 144 INCREMENT 1] = npred[CV()] + NVL(LOG(10,p[CV()-1]),dpred[CV()-1]), pred[FOR m FROM 1 TO 144 INCREMENT 1] = POWER(10,dpred[CV()]) ) ORDER BY m; |
因为我们有预测期的真实数据,而且我们也在构建模型前将数据转换了。
RULES
段描述如何计算每个度量:
l
Ap
代表训练期的原始数值,从
132-144
的预测期置为
null
l
Tp
转换数据,转换方法:
LOG
,差分,
Z-SCORE
常态。
l
Tpred
计算训练的预测值并且使用
PREDICTION
函数预测。规则对
tp
序列可用时使用真实转换数据,否则使用预测数据。提供的预测数据与转换后的数据的方式一样
l
Npred
反向转换预测值的常态化数据
l
Dpred
反向转换差分。如果真实值可用则使用实际值,否则使用
dpred
前值。因为在构建模型时,差分是在
LOG
转换后执行的,当使用时机值时,我们需要
LOG
转换序列。
l
Pred
是与未转换的原值相同视角的预测数据,使用
POWER
函数来反向转换
LOG
的结果。
下图是一个各种方法的比较结果
Model
|
RMSE Training
|
MAE Training
|
RMSE Test
|
MAE Test
|
SVM (one-step)
|
3.7
|
3.2
|
18.9
|
13.4
|
AR (one-step)
|
10.1
|
7.7
|
19.3
|
15.2
|
SVM (multi-step)
|
3.7
|
3.2
|
11.5
|
9.3
|
示例2:电力负载预测
在这个例子中,我们要解决的问题是预测每天最大的负载,基于早先的电力负载值和额外的数据。训练数据由
1997-1998
年间的半小时间隔的负载。对这个例子有一个比赛,结果在
here
获胜者是
Lin
,用的是支持向量机回归算法。这个例子中,我会按照
Lin
的方法。使用下面的输入来构建模型:
l
以前的
7
个最大负载值
l
每周的天数
l
0-1
的标识来指出是否是假日
首先,我从不同文件
(1997,1998,1/1999)
中将数据合并在一个文件中,然后将数据加载到表
load_all
中,这个表结构为
day_id number,day Date,day_week varchar2(1),max_load number,holiday number
,
day_id
是一个从
1-761
的顺序的范围;
day_week
是对于给定
day
的
week
的天数,
day_week
很容易通过使用
to_char(day,’d’)
来计算出来。
因为序列没有任何明显的趋势或差异,我们只需要常态化
(normalize)max_load
列来准备构建模型的数据。我常态化序列到
[0,1]
的范围内。
CREATE VIEW load_norm AS SELECT day_id, day, max_load, holiday, day_week, (max_load-464)/(876-464) tl FROM load_all |
464
和
876
分别是在
97-1-1
到
99-1-31
范围内的
Min(max_load)
和
MAX(max_load)
。
在构建模型前,我们需要选择哪个准备序列的迟滞
(lag)
值作为未来值的预测器。正如在
part2
中所讨论的,这等同于选择时间窗口并且在窗口中包含所有的
lag
值。我们使用
part2
中的
xcorr
过程来计算自相关。带有头
20
个
lag
的
Tl
序列可以计算:
BEGIN xcorr(p_in_table=> 'load_norm', p_out_table => 'my_corr', p_seq_col => 'day_id', p_base_col=> 'tl', p_lag_col => 'tl', p_max_lag => 20); END; / |
最大的自相关期限是发生在
lag7
,这看上去好像是有原因的,我们希望电力负载能表现出很明显的周期性,基于这种分析,我们可能包括在窗口大小为
7
的预测器。
常态化序列并且选择
lagged
变量后,我们准备创模型来预测未来值。首先我们需要创建带
7
个
lagged
变量的视图:
CREATE VIEW load_lag as SELECT day_id, day, max_load, tl, holiday, day_week, LAG(tl, 1)OVER (ORDER BY day) L1, LAG(tl, 2)OVER (ORDER BY day) L2, LAG(tl, 3)OVER (ORDER BY day) L3, LAG(tl, 4)OVER (ORDER BY day) L4, LAG(tl, 5)OVER (ORDER BY day) L5, LAG(tl, 6)OVER (ORDER BY day) L6, LAG(tl, 7)OVER (ORDER BY day) L7 FROM load_norm |
接下来,我们创建训练集,使用
load_lag
视图中行的子集。因为我们想测试模型的预测能力,因此我们需要在历史数据上进行训练并且保存预测结果用于测试。抽取
731-761
用于测试数据集。我们也需要过滤头
7
行将一些
lagged
变量值置为
NULL
。不像
Lin
,我我觉得包括
8
月和
9
月的训练集更好,下面的视图创建训练集。
CREATE VIEW load_train AS SELECT day_id, day, max_load, tl, holiday, day_week, L1, L2, L3, L4, L5, L6, L7 FROM load_lag WHERE day_id < 731 and day_id > 7; |
最后,使用
ODMR
或数据挖掘
API
构建模型,对于这个例子,我们使用
PL/SQL API
和默认的
SVM
回归算法构建模型。
BEGIN DBMS_DATA_MINING.CREATE_MODEL( model_name => 'load_SVM', mining_function => dbms_data_mining.regression, data_table_name => 'load_train', case_id_column_name => 'day_id', target_column_name=> 'tl'); END; |
上面的语句创建了
SVM
回归模型,叫做
load_SVM
,使用
load_train
视图的数据作为训练数据。
要应用预测,我们只需要应用模型到新数据上。例如,应用
load_SVM
模型到
load_lag
视图上为训练数据
(1-730)
和测试数据
(731-761)
生成多步预测。
SELECT m day_id, p max_load, pred FROM load_all a MODEL DIMENSION BY (day_id m) MEASURES (a.max_load p, a.day_week dw, a.holiday h, CAST(NULL AS NUMBER) ap, CAST(NULL AS NUMBER) tp, CAST(NULL AS NUMBER) tpred, CAST(NULL AS NUMBER) pred) RULES( ap [FOR m FROM 1 TO 730 INCREMENT 1] = p[CV()], h[FOR m FROM 731 TO 761 INCREMENT 1] = 0, tp[FOR m FROM 1 TO 730 INCREMENT 1] = (ap[CV()] - 464)/(876-464), tpred[FOR m FROM 1 TO 761 INCREMENT 1] = PREDICTION(load_SVM USING NVL(tp[CV()-1],tpred[CV()-1]) as l1, NVL(tp[CV()-2],tpred[CV()-2]) as l2, NVL(tp[CV()-3],tpred[CV()-3]) as l3, NVL(tp[CV()-4],tpred[CV()-4]) as l4, NVL(tp[CV()-5],tpred[CV()-5]) as l5, NVL(tp[CV()-6],tpred[CV()-6]) as l6, NVL(tp[CV()-7],tpred[CV()-7]) as l7, dw[CV()] as day_week, h[CV()] as holiday ), pred[FOR m FROM 1 TO 761 INCREMENT 1] = 464 + (876-464) * tpred[CV()] ) ORDER BY m; |
因为我们只应用了常态化准备数据,因此我们只需要在
RULE
段做一步反向转换即可。这个模型与以前模型的主要差别在于额外的度量捕捉输入:
day_week
和
holiday
。与
Lin
一样,我在查询的第二个规则中设置预测期的
holiday flag(h)
为
0
。
Model
|
MAPE Training
|
MAX Load Error Training
|
MAPE Test
|
MAX Load Error Test
|
SVM (one-step)
|
2.2
|
77.4
|
1.7
|
39.1
|
SVM (multi-step)
|
2.2
|
77.0
|
1.9
|
50.4
|
通过上面的表,我们可以看到,多步方法执行的更好,
MAPE
值一样,
MAX Load forecast erro
比单步大。
总结
这系列文章帮助我们如何使用
Oracle Data Mining
来创建强大的时序预测模型。
附记
在
part2
的航班例子中,常态化
shift
和衡量参数使用完整数据进行计算。更好地方法是使用训练数据计算这些参数。这些可以选择的常态化方案也影响航班例子的结果。我在
part2
中修改了常态化,以便只适用训练数据计算出来的参数。所有的例子和数据在这里可以获得。
here