摘要

我们考虑在电网负荷预测中的应用,其中广义加性模型是合适的,但数据集的大小可能会使其在现有方法中难以使用。因此,在模型中的平滑项用惩罚回归样条表示的情况下,我们为大数据集开发了实用的广义加性模型拟合方法。这些方法使用迭代更新方案来获得模型矩阵的因子,而在任何时候只需要计算模型矩阵的子块。我们表明,可以以合理的方式进行有效的平滑参数估计。当新数据可用时,电网负荷预测问题需要更新模型拟合,以及处理电网负荷中剩余自相关的一些方法。为这些问题提供了方法,并介绍了并行实现。这些方法允许使用适当的计算机硬件来估计大型数据集的广义可加模型,而电网负荷预测问题说明了降秩样条平滑方法在处理复杂建模问题中的实用性。

1.简介

包含数万到数百万个响应观测值的回归问题现在很常见。有时,这样大的数据集需要全新的建模方法,但有时现有的模型类是合适的,只要它们在计算上可行。本文考虑使用适度的计算机硬件使广义可加模型(GAM)估计对大数据集可行的问题,在这种情况下,平滑参数必须作为模型拟合的一部分进行估计。

我们的动力来自一个特别的应用程序。图。1显示了自2002年9月1日起,法国国家电网每半小时的千兆瓦负荷。法国能源公司Electricitéde France(EDF)在使用GAM进行短期负荷预测方面取得了相当大的成功,该预测基于一些协变量,尤其是24小时前的负荷。然而,使用现有的GAM拟合方法,一次拟合整个数据集在计算上是不可行的,取而代之,每天每半小时必须拟合48个单独的模型(Pierrot和Goude,2011). 显然,最好使用单一模型。

 (a) 自2002年9月1日以来,法国国家电网的负荷为每天半小时(单位:千兆瓦),以及(b)数据最后100天的半小时千兆瓦负荷
图1。

(a) 自2002年9月1日以来,法国国家电网的负荷为每天半小时(单位:千兆瓦),以及(b)数据最后100天的半小时千兆瓦负荷

尽管对现有的GAM估计方法具有挑战性,但按照目前的标准,我们的激励示例相对温和。例如,许多统计建模工作都致力于阐明空气污染与呼吸系统死亡率之间的关系。一种方法是使用GAM将死亡率分解为在时间上平滑变化的背景成分和污染物影响(平滑或线性),同时将观察到的死亡计数视为泊松分布。彭和韦尔蒂(2004)收集了美国108个城市约5000天的每日污染死亡率数据,分为三个年龄组。木材分析(2006)第5.3节(仅针对芝加哥)提出了非常强烈的臭氧-温度相互作用效应,如果可重复,这将非常重要。然而,这种影响对几天的数据很敏感,因此应该在彭和韦尔蒂数据集中剩下的美国城市进行测试。理想情况下,我们希望同时将GAM拟合到数据集中的所有120万个观测值。如在线支持材料所示,这种拟合远远超出了现有GAM估计方法的范围,但在使用以下开发的方法的中等计算硬件上是可行的。

GAM和相关模型(例如Wood)的当前装配方法(2011))当应用于包含数万个观测值的数据集时,它们是相当有效和健壮的,但在这一点上,它们往往会变得过于内存密集,因此包含数十万或数百万数据的更大数据集是遥不可及的。困难在于模型的模型矩阵可能变得太大:如果n个第页分别是模型矩阵的行数和列数,以及M(M)是平滑参数的数量,则GAM拟合方法的内存要求通常为O(运行)(2),可能会变得太大而无法处理。在这里,我们将展示如何使用更新模型矩阵分解的简单策略来避免在GAM上下文中形成整个模型矩阵。最重要的是,我们展示了如何在这种情况下调整平滑参数估计方法。

我们考虑的一般模型类可以写成

{E类()}=A类θ+j个L(左)j个(f)j个
(1)

哪里是其中之一n个从指数族分布(或至少在标度参数内已知的均值-方差关系)中观察单变量响应变量,是已知的平滑单调链接函数,A类是一个n个-行模型矩阵,θ是未知参数的向量,L(左)ij公司是已知的线性函数(f)j个一个或多个变量的未知平滑函数,平滑度未知。与每个(f)j个是偏离平稳性的度量J型j个((f)).

最常见的方程式示例(1)是GAM(Hastie和Tibshirani,1986, 1990)当L(左)ij公司是评价函数,但其他例子有变系数模型、泛函广义线性模型(如马克思和艾勒(1999))和结构化加性回归模型(例如Fahrmeir. (2004))。有多种方法可用于估计模型类的特定成员。我们将重点关注以下情况:(f)j个用中等秩惩罚回归样条表示(例如Parker和Rice(1985)艾勒和马克思(1996))。在这种情况下,可以为整个模型类获得非常有效的计算方法(例如Wood(2000)),通过惩罚迭代加权最小二乘法进行拟合,通过广义交叉验证(GCV)、限制最大似然(REML)或类似方法进行平滑度选择(参见Wood(2011)).

2.高斯身份链接情况

首先考虑以下情况:具有方差的独立正态分布ϕ、和是标识函数。这个(f)j个每个都通过使用线性基展开来表示(例如B类-样条或薄板回归样条基),以及J型j个被选为基系数的二次型。在这种情况下,预期响应的模型可以重写为

E类()=X(X)β
(2)

哪里n个×第页模型矩阵X(X)包含A类和评估的基函数,以及β包含θ以及所有基础系数。我们假设第页<n个,因为这里介绍的方法只在这种情况下才有实际意义,我们估计β通过最小化

负极X(X)β2+j个λj个βT型S公司j个β
(3)

哪里S公司j个是已知系数的矩阵,如下所示J型j个((f)j个)=βT型S公司j个β(S公司j个第页×第页,但其非零块通常小于此值),并且λj个是控制拟合-平滑权衡的平滑参数(f)j个(这里的符号有点草率,因为可能有几个平滑参数与一个参数相关(f)j个.)给定λ,表达式(3)可以很容易地最小化,以给出系数估计值β^λ.

λ更尴尬。一种基于尝试最小化预测误差的方法是GCV,它寻求最小化

V(V)(λ)=n个负极X(X)β^λ2{n个负极信托收据(如果λ)}2.

关于平滑参数,其中tr(如果λ)是模型的有效自由度,如果λ=(X(X)T型X(X)+S公司λ)−1X(X)T型X(X)S公司λj个λj个S公司j个其他可能性为REML,V(V)第页,或马洛斯的C类第页这在附录A中有介绍。牛顿法通常用于优化V(V)关于日志(λ),使用β^λ通过直接最小化每次试验的表达式(3)获得λ:木材(2004,2011)提供计算稳定的数值方法的详细信息。

现在假设模型矩阵是第一个二维码分解为正交列n个×第页因素和上三角形第页×第页因素R(右)以便X(X)=二维码.如果我们还形成(f)=T型和‖第页2=‖2−‖(f)2则表达式(3)变为

(f)负极R(右)β2+第页2+j个λj个βT型S公司j个β
(4)

相当常规的计算表明

V(V)(λ)=n个(f)负极R(右)β^λ2+第页2{n个负极信托收据(如果λ)}2

如果λ=(R(右)T型R(右)+S公司λ)−1R(右)T型R(右).

关键是一旦我们R(右),(f)和‖第页2然后我们就有了安装所需的一切,并且X(X)不再扮演其他角色。因此,如果我们不需要形成X(X)总的来说,我们可以在不产生高计算机内存成本的情况下估计模型。附录A表明,使用Mallows的C类第页和REML,并讨论了一些潜在的替代方法。

事实上R(右),(f)和‖第页2可以以只需要X(X)使用基于迭代更新二维码-分解,或使用Choleski分解方法不太稳定。附录B提供了完整的细节,并表明该方法自然会为大型可加模型带来高效的在线更新方法。

2.1. 相关错误

第节中讨论的短期负荷预测问题4还需要我们对残差自相关和简单的AR进行建模(第页)相关结构很容易处理。首先将高斯身份链接模型(2)修改为=X(X)β+e(电子),其中协方差矩阵e(电子)ϕΣΣ是自动累进AR(第页)相关矩阵。然后是Choleski因子C类属于Σ−1是带状的并且ɛ=总工程师独立且分布一致N个(0,ϕ). 因此,如果~=Cy公司X(X)~=CX公司,我们有

~=X(X)~β+ε,
(5)

其形式为(2),因此可以使用前面章节的方法进行估算β唯一的修改是,如果使用REML进行估算ρ则log-REML分数必须通过以下方式进行修正C类,但鉴于此C类是三角形的,所需的对数行列式很容易获得。从计算上来说,简单的一维搜索可以用于ρ,每个ρ-需要重新安装模型的值。请注意C类事实上,它可以在不成形的情况下获得Σ−1是什么构成AR(第页)计算上可行的模型:X(X)~~包括计算成本低廉的相邻行的加权差分X(X)而不是昂贵的全矩阵乘法。

3.广义加性模型拟合

在广义情况下,未知函数及其惩罚的表示与简单高斯恒等式链路情况完全相同。所有的变化是,该模型变成了一个过度参数化的广义线性模型,

{E类()}=X(X)β,

通过惩罚似然最大化来估计,代替惩罚最小二乘法(参见示例Green和Silverman(1994))。用于最大化惩罚可能性的算法被迭代地重新加权最小二乘法(PIRLS)惩罚,其进行如下,其中V(V)是这样的函数:var()=¼V(μ)和μ=E类().

第一次初始化μ^=+ξη^=(μ^)哪里ξ添加少量(通常为0)以确保(μ^)存在。然后重复以下步骤以达到收敛。

  • 步骤1:表单z(z)=(μ^)(负极μ^)+η^w个=V(V)(μ^)负极1/2(μ^)负极1.

  • 第2步:将w个在对角矩阵中W公司,最小化表达式(3)的加权版本,
    Wz公司负极WX公司β2+j个λj个βT型S公司j个β,
    关于β以获得β^和更新η^=X(X)β^、和μ^=负极1(η^).

对于中等大小的数据集,迭代PIRLS算法以每次试验收敛是最可靠的λ,并进行估算λ通过使用GCV的通用版本,C类第页或拉普拉斯近似REML(见Wood(2008,2011))。在大数据集情况下,这种方法的缺点是需要几倍于X(X)有效计算平滑度选择准则的导数。为了避免如此高的存储成本,我们可以返回到更早的方法,最初是由于Gu(1992). 这只是使用GCV,C类第页或REML在PIRLS算法的每个步骤中选择工作线性模型的平滑参数。顾(1992,2002)称之为“面向性能的迭代”,它与Breslow和Clayton的非常相似(1993)受惩罚的准自杀。该方法通常收敛,尽管不能保证收敛,但促进收敛问题的病态条件往往会随着时间的增加而减少n个.

3.1、。面向性能的大型数据集迭代

可以使用二维码-矩阵附录B的更新方法WX公司,在PIRLS算法的每个步骤。这意味着形成模型矩阵(的子矩阵)所需的计算必须在每个步骤中重复,但这些计算是O(运行)(净现值)而不是O(运行)(净现值2)的二维码-分解,所以对于许多平滑基,操作成本并不重要。形式上,算法如下。

3.1.1. 初始化

x个表示与响应变量相关的协变量向量哪里= 1,…,n个,并将整数从1除以n个进入之内M(M)非重叠子集γ1,…,γM(M)大小近似相等的(soŞγ={1,…,n个}和γj个γ=∅j个).M(M)选择此选项是为了避免计算机内存不足。η¯=(+ξ)(带有ξ如前一节所述)。设置PIRLS迭代索引q个=0和D类=0(或任何常数)。执行为平滑项设置基础所需的任何初始化步骤。

3.1.2. 迭代

  • 步骤1:设置D类古老的=D类,R(右)为0×第页矩阵,(f)0向量,D类=0和第页= 0.

  • 第2步:对重复以下步骤(a)-(f)k个= 1,…,M(M).

    • (a)

      设置(f)0=(f)R(右)0=R(右).

    • (b)

      形成模型子矩阵X(X)k个对于协变量集{x个:γk个}.

    • (c)

      如果q个>0表单η^=X(X)k个β^λ否则η^=η¯γk个.

    • (d)

      表格μ^=负极1(η^),z(z)=(μ^)(负极μ^)+η^w个=V(V)(μ^)负极1/2(μ^)负极1γk个.让z(z)是包含这些的向量z(z)-值和W公司是对应的对角矩阵w个-值。

    • (e)

      设置第页第页+‖Wz公司2,计算当前数据子集的偏差残差,并将其平方和加到D类.

    • (f)
      表格
      R(右)=(R(右)0W公司X(X)k个)
      (f)=T型((f)0W公司z(z))
      然后丢弃.
  • 步骤3:套‖第页2=第页——————(f)2.

  • 第4步:如果q个>通过比较电流偏差进行0收敛测试D类与之前的偏差D类古老的.如果已达到收敛,则停止(或q个已超过某些预先确定的极限,表明发生了故障)。

  • 第5步:估计λ通过优化V(V)第页/C类第页与第节完全相同2。这也会产生β^λ.

  • 第6步:q个q个+1.

在汇合处,最终β^λλ是系数和平滑参数估计值。进一步的推断是基于贝叶斯近似最有用的

βN个{β^λ,(R(右)T型R(右)+S公司λ)负极1ϕ},

如果有ϕ-作为REML优化的一部分或作为

ϕ^=(f)负极R(右)β^λ2+第页2n个负极信托收据(如果λ)

(参见木材示例(2006)更多详细信息)。

注意,步骤2可以按照附录B中的描述进行并行化γ被分组为大小相等的非重叠集,以分配给不同的处理器。步骤2在每个处理器上运行γ,导致R(右)(f)来自每个处理器。应用附录B中的表达式(7),得出所需的R(右)(f)此外,在步骤5中,计算了Wood平滑参数优化方法的操作数(2011)可能由于以下原因而大幅下降V(V)在平滑参数优化中未加权。

基于Choleski的替代方案只需替换迭代的三个步骤,如下所示。

  • 步骤1:设置D类古老的=D类,R(右)成为一名第页×第页0s矩阵,(f)第页-向量0s,D类=0和第页= 0.

  • 步骤2(f):设置R(右)=R(右)0+X(X)k个T型W公司X(X)k个(f)=(f)0+X(X)k个T型Wz公司.

  • 步骤3:替换R(右)(其中真正包含X(X)T型WX公司)通过其Choleski分解,然后替换(f)通过R(右)−1(f).套装‖第页2=第页——————(f)2.

同样,在这种情况下,步骤2很容易并行化。

3.1.3. 初始值的子采样

使用完整的数据集运行PIRLS算法的早期步骤在计算上是浪费的,因为这相当于浪费精力来精确拟合已知错误的工作模型。因此,通常明智的做法是,首先在5–10%的随机数据子样本上估计模型,然后使用结果βλ-估计值作为拟合完整数据的起始值。实际上,在拟合完整数据时,此技巧通常节省PIRLS算法的一到两个步骤。

3.2. 调整平滑度选择步骤

应用GCV或C类第页对于工作模型,在PIRLS迭代的每个步骤:这些标准所需的假设适用于工作模型。

REML(或ML)不如工作数据那么简单z(z)可能与推导所需的常态相去甚远V(V)第页。但是,对于具有n个第页中心极限定理意味着(f)=T型z(z)(其中z(z)现在指的是整个工作数据n个-向量)将倾向于N个(R(右)β,ϕ)分配。基于密度的REML得分(f)就是那个时候

V(V)第页*(λ)=(f)负极R(右)β^λ2+β^λT型S公司λβ^λ2ϕ+第页负极M(M)第页2日志(2πϕ)+日志|R(右)T型R(右)+S公司λ|负极日志|S公司λ|+2

哪里|S公司λ|+是的正特征值的乘积S公司λ其中有M(M)第页零特征值。所以V(V)第页*具有确切的形式V(V)第页参见附录A,但带有‖第页2≡0和n个设置为第页.对于任何固定ϕ,V(V)第页V(V)第页*明显被相同的λ然而,如果ϕ是未知的,那么必须以某种方式进行估计。优化V(V)第页*关于ϕ显然不是一个好的选择,因为缺席学期‖第页2携带有关的信息ϕ然而,

ϕ^=(f)负极R(右)β^λ2+第页2+β^λT型S公司λβ^λn个负极M(M)第页

是的估计量ϕ可以通过类比REML估计来激发ϕ或者作为一个简单的矩估计量(通常,估计量的分子是皮尔逊统计量加上平滑惩罚)。现在,因为ϕ^很容易被视为V(V)第页,而V(V)第页V(V)第页*以相同的方式最小化λ,然后最小化V(V)第页*关于λ使用时ϕ^作为ϕ与查找相同λ^ϕ^通过最小化V(V)第页(带有R(右),(f)和‖第页2由PIRLS计算,如第节所示3.1.).

莱斯和奥格登最近的工作(2009)表明V(V)第页不太容易出现多个局部极小值V(V)(或者可能是密切相关的C类第页). 在以性能为导向的迭代环境中,这表明使用REML可能会促进迭代收敛,因为迭代在多个优化之间循环的范围较小。

3.3. 这个第页<n个假设

当对如此大的数据集使用GAM模型时,计算可行性取决于使用降秩平滑的能力,以便第页大大小于n个。但一个明显的问题是,这样假设是否合理第页将比n个,正如这里给出的方法所做的那样,隐式地。显然,如果平滑被用作随机场以消除剩余相关性,或者当附加数据扩展时间轴时,对于平滑时间,这种假设是不合理的。否则,所有理论证据都表明,平滑基只需要随着样本大小缓慢增长:例如,考虑具有等距节点的三次回归样条曲线的情况。众所周知,三次样条曲线的平均平方偏差为O(运行)(小时8)=O(运行)(k个−8)其中小时是结间距k个节数(例如de Boor(1978))。根据基本回归理论,样条曲线的平均方差为O(运行)(k个/n个)其中n个是数据的数量。为了避免偏差或方差支配另一方,并给出次优均方误差n个→∞, 我们应该选择k个将平方偏差和方差的阶数相等,即。k个n个第1页,共9页此比率也可用于惩罚,因为在任何有限的样本量下,我们选择惩罚程度来减少相对于纯回归的均方误差,尽管此论点本身并不表示替代比率在惩罚下可能不是最佳的。确实,在惩罚之下,顾和金(2002)建议基本维度的比例应为n个2/9也就是说,100万个数据只需要大约1000个所需系数数量的五倍。

4.短期电网负荷预测

如第节所述1,图。1图中显示了法国国家电网的负荷(单位为千兆瓦)。EDF已成功建立了一天一次的电网负荷预测模型,该模型基于每天半小时将数据拆分,并将GAM拟合到48个子集中的每一个子集。尽管使用现有方法进行估算是可行的,但使用48个独立模型存在三个实际缺点。

  • (a)

    它无法有效地利用信息,因为没有利用相邻半小时内数据之间的相关性。

  • (b)

    由于半小时期间的模型连续性没有得到强制执行,因此存在解释困难,这有点不合实际。

  • (c)

    运营预测模型必须在统计上是稳定的,同时模型的预测目的表明,此类模型的平滑度估计应使用预测误差标准,如GCV。然后出现了一个困难,即已知GCV和相关标准会产生一小部分实质性飞越,随着样本量的减少,问题也会增加(例如,见Reiss和Ogden(2009))。将单独的模型拟合到48个数据子集会加剧此问题的暴露,从而降低模型稳定性,并且每次更新模型估计值时,运营预测员都要承担大量的模型检查任务。将一个模型拟合到所有数据,可以大大减少过拟合的范围,同时将相关的模型检查任务减少到可管理的程度。

本文所述方法开发的主要动机是允许使用单个模型来代替48个单独的模型,此处给出了此类模型的一个版本。

图。1这些数据来自EDF收集的数据集,其中还包括气象数据(法国气象局对气温(摄氏度)和八分之一天空云量的1天预报)、日历和电价信息。虽然负荷的决定因素相对复杂,因此必须进行统计预测,但EDF采取了明智的预防措施,坚持预测模型必须具有物理意义上的可解释影响。这在一定程度上有助于在模型拟合数据集中通常适用的条件之外,对异常事件进行操作性预测。

以往的预测经验和探索性数据分析表明,基于提前24小时的电网负荷、预测温度、提前24小时和48小时的实际温度、预测云量和一年中的时间,可以获得良好的预测性能。一般来说,负载最低约为15C、 当供热需求随温度降低而增加时,增幅低于该值,而当温度高于15℃时,增幅较小这是由于空调。用于预测的温度是法国的平均温度,以用电量加权。对除云层外相似的天数进行比较表明,云层的另一个重要直接影响,这可能归因于建筑物的被动太阳能加热效应和照明效应。

另一个重要影响是,EDF通过减少对大用户的经济刺激,管理了冬季预期的超高需求日。提前1天,EDF Services在下午7点对此类优惠政策适用的“特殊关税日”的预期需求减少进行预测,该预测也包括在内。基于这些考虑、探索性分析、每天每半小时使用单独模型的经验以及一些模型选择后,建议使用以下模型:

L(左)=γj个()+(f)k个()(,L(左)负极48)+1(t吨)+2(,玩具)+(T型,)+4(T型.24,T型.48)+5()+装货单小时()+e(电子)
(6)

如果观察到是从星期几开始的j个、和日班k.符号j个()表示索引j个本身是索引的函数误差项被建模为AR(1):e(电子)=ρe−1+ɛɛN个(0,σ2). 在这里L(左)是指第半小时;是一天的半小时周期(1到48之间的整数;此后为“瞬间”);t吨是自2002年9月1日起的时间;玩具是一年中从0到8760的小时数;云是法国上空云量的指数;T型是温度,T.24和T.48温度测量值分别滞后24和48小时。日班是hh、hw、ww或hw的一种,取决于所讨论的一天是假日之后的假日,还是假日之后的工作日,等等。其想法是,一天的负荷取决于前一天负荷的方式在很大程度上取决于所涉及的日是工作日还是假日。在这种模式下,周末是假日。ST是预测的特殊关税减载,正常日为0,每个特殊关税日的整体为单个数字。这个k个(f)j个所有平滑函数都表示为惩罚回归样条。四个人(f)k个是三次回归样条(循环);2是相似的张量积,但每一个都是秩120;4是秩45张量积样条;15是三次样条曲线和小时是一个非规范化函数,相当于一天中半小时的48级因子变量。

模型(6)中的个别术语主要基于EDF经验预期的影响,在某些情况下,如温度影响先验的理由;然而,有一些实质性的结构假设并不明显先验的一个考虑因素是,这些影响是否最好被视为加法或乘法,这可以通过使用或不使用对数链估计模型来解决。完成后,假设ρ=0,估计值第页2这两个模型的三个显著数字相同,预测性能没有明显差异。因此,我们决定使用加性版本,这使得相关性的处理比其他版本更容易。另一个问题是季节性的处理方式。我们在一年中的这个时候使用了循环效应,但另一种优雅的季节性方法是Eilers的变系数方法. (2008),其中通过以下形式的截断傅里叶级数处理季节性

k个=1K(K)(f)2k个负极1(t吨,)(2k个π/T型)+(f)2k个(t吨,)余弦(2k个π/T型)

其中(f)j个(t吨,)术语是时间和一天中瞬间的平滑变化函数,它控制一个周期的相位和振幅T型循环。我们用这样的术语来处理季节性的模型进行了试验,但对于这些数据,在以下方面,我们只能产生比模型(6)稍差的性能R(右)2以及预测误差性能。

为了说明模型的有效性,使用截至2008年8月31日的数据对模型进行了估计,然后使用附录B中的在线更新方法用每一天的数据更新模型估计值,对下一年进行了为期一天的预测。在操作预测中,处理银行假日特别的并且通常被排除在常规的自动预测之外,所以我们在这里也这样做了。鉴于建模的预测性质,GCV用于平滑参数和AR参数估计。请注意,在进行初始模型拟合时,有必要以允许的方式设置底座1(t吨)将一个域运行到2009年9月,以便在预测期结束之前继续保持适当的基础。设置起来没有问题,但如果我们使用基于Choleski的更新方案作为初始方案,则可能会导致数值问题X(X)T型X(X)然后可能是等级不足,或者几乎如此。

图。2显示了初始拟合数据的残差图,表明主要误差出现在周一上午的预测中,较大的异常值主要出现在白天,而残差在冬季负荷较高时往往更大。图。2(d)显示了最后一年数据的1天预测残差。预测的残差图与初始拟合的残差曲线相似。AR参数估计为0.98,图。2还示出了参数设置为0的简化模型的预测残差。模型的平均绝对百分比误差MAPE和平方根误差RMSE图(表1). 图。3(a)直接显示了模型的拟合荷载与图中所示期间的观测荷载叠加。1(b).

 (a) 每天半小时“瞬间”的残差,(b)剩余平均绝对百分比误差MAPE与瞬间(,周一;,周二;,周三;,周四;,周五;,周六;,周日)的残差残差)
图2。

(a) 每天半小时“瞬间”的残差,(b)残差平均绝对百分比误差MAPE(图解的,星期一;图解的,星期二;图解的,星期三;图解的,星期四;图解的,星期五;图解的,周六;图解的,星期日),(c)时间残差和(d)时间预测残差(图解的模型忽略残差自校正;图解的,带AR(1)残差的完整模型)

 (a) 最后100天数据的拟合半小时千兆瓦负荷()覆盖在观测负荷上(),(b)忽略相关性时的自相关函数,以及(c)AR(0.98)误差模型标准化残差的自相关功能
图3。

(a) 最后100天数据的半小时千兆瓦负荷(图解的)覆盖在观测荷载上(图解的),(b)忽略相关性时的自相关函数,以及(c)AR(0.98)误差模型中标准化残差的自相关功能

表1

有相关性和无相关性估计的模型(6)的模型拟合和预测性能的比较

模型变量RMSE(Mw)MAPE(%)
ρ=0拟合数据8311.17
ρ=0预测12201.87
ρ=0.98拟合10241.46
ρ=0.98预测11561.62
模型变量RMSE(Mw)MAPE(%)
ρ=0拟合数据8311.17
ρ=0预测12201.87
ρ=0.98适合10241.46
ρ=0.98预测11561.62
表1

模型(6)的模型拟合和预测性能比较(有相关性和无相关性估计)

模型变量RMSE(Mw)MAPE(%)
ρ=0拟合数据8311.17
ρ=0预测12201.87
ρ=0.98适合10241.46
ρ=0.98预测11561.62
模型变量RMSE(Mw)MAPE(%)
ρ=0拟合数据8311.17
ρ=0预测12201.87
ρ=0.98适合10241.46
ρ=0.98预测11561.62

AR模型的预测误差较低,拟合集和预测集之间的失配也低得多。更大的不匹配ρ当忽略相关性时,=0几乎肯定是由于过拟合造成的:具有相关性的模型的有效自由度是无相关性模型的83%。这强调了第节中开发的方法的实际重要性2.1.图3(b)3(c)显示模型残差的自相关函数ρ=0,对于标准化残差,当ρ= 0.98. 显然AR(1)模型带来了实质性的改进,但还有进一步完善的空间。模型(6)的性能与Pierrot和Goude的性能具有竞争力(2011)但它具有上面列出的单模型拟合的三个实际优点,同时还可以在新数据可用时简化模型估计的更新过程。为了理解后一个优点,我们使用一台运行LINUX(即零售价低于600美元的个人电脑)的机器,使用3.1-GHz Intel i3 540处理器和3.7 GB随机访问内存,比较了计算时间。与模型(6)等价的48个独立模型,在没有AR残差的情况下,使用包mgcv中的gam估计大约需要1小时,而使用本文的方法对模型(6ρ使用附录B进行后续更新只需不到2分钟,而之前需要对所有48款车型进行完全改装。

一个明显的问题是,当模型已经适合如此长的数据集时,是否需要每天进行更新。这样做的原因是,即使在相当长的数据运行中,也有许多条件组合可能无法很好地采样,预测变量本身也高度相关。这意味着,如果条件在几天内不寻常,那么最后几天的信息可能会在这些条件下构成有关负载的信息的不可忽略的比例,因此可能会对这些不寻常协变量附近的平滑函数的估计产生不可忽视的影响。由于这个原因,在用于预测的模型拟合中排除最新的数据通常是不希望的。

虽然我们开发了本文所述的方法,直接回应了使用现有GAM估计方法解决该问题的困难,但模型(6)可以通过Wood方法进行估计(2011)提供约8 GB内存。然而,这样做的速度还不到这里提出的方法的十分之一,只有通过完全改装才能进行更新。对于在模型开发过程中考虑的一些更复杂的模型,使用更大的基数和更复杂的日类型分类,将需要更多的内存,如果没有新方法,按区域建模某些影响的未来目标将完全无法实现。在线支持材料分析了第节中讨论的空气污染数据1并提供了一个示例,该示例在没有这里考虑的方法的情况下,远远超出了标准个人计算硬件。

为了说明可解释性(f)k个如图所示。4注意,一天中任何时候,一个工作日的负荷与前一个工作天的负荷都是相对线性关系,与平均值之间存在明显的回归。节假日呈现出不同的模式,前一天的低功耗并不能很好地预测,但在50 Gw以上,线性关系明显。节假日后的工作日以周一为主,与节假日的工作日呈现出相反的形状,在低负荷时影响最为明显。

 针对第一天时间(t)的滞后负荷效应,针对(a)工作日-工作日,(b)工作日–节假日,(c)节假日-工作日和(d)节假日–工作日过渡
图4。

滞后负荷对一天中时间的影响(t吨)对于(a)工作日-工作日,(b)工作日–节假日,(c)节假日-工作日和(d)节假日–工作日过渡

总之,本文开发的方法使我们能够通过同时拟合所有可用数据来提高EDF预测模型的稳定性和可解释性,从而允许在新数据可用时进行自相关建模和有效估计更新。

5.结论

出于提高EDF所用预测模型的速度和稳定性的实际需要,我们展示了GAM如何应用于比迄今为止通常可能应用的更大的数据集。我们的方法的一个特别优点是,它可以通过对现有方法进行相对直接的扩展来实现,同时在可以建模的数据集大小和某些情况下的拟合速度方面都有了很大的改进。平滑泊松回归空气污染示例,在第节中介绍1并进一步开发了在线支持材料,为实际改进提供了一个鲜明的例子。在这种情况下,如果一次性形成,仅模型矩阵就需要7 GB以上的存储空间,但我们可以使用不到1 GB的存储空间来适应模型。我们知道,没有其他公开可用的方法可以适用于我们用于空气污染数据集的广泛结构模型。估计(在第节中描述的廉价计算机上花费不到12分钟4)如果存储没有问题,它的速度也比现有方法预计的管理速度快100倍左右。

除了推动这项工作的领域外,遥感、基因技术、金融和信息学的发展都在产生越来越多的超大数据集。有时,此类数据需要全新的分析方法,但有时,成熟的现有模型类也很有用,所需的只是使其在计算上可行。本文中提供的方法可以做到这一点,因此除了这里给出的直接应用程序之外,它应该具有广泛的适用性。

本文讨论的方法在R包mgcv(Wood,2009)作为函数bam。

鸣谢

作者感谢两位审稿人,副编辑和联合编辑的评论,这有助于改进本文。西蒙·伍德(Simon Wood)由工程和物理科学研究委员会(Engineering and Physical Sciences Research Council)拨款EP/K005251/1资助。

工具书类

贝茨
,
D类
.和
梅克勒
,
M(M)
. (
2013
)
矩阵:稀疏和稠密矩阵类和方法
.
R软件包版本1.0-12
。(可从http://CRAN.R-project.org/package=矩阵

德布尔
,
C类
. (
1978
)
花键实用指南
.
纽约
:
施普林格
.

布雷斯洛
,
东北。
克莱顿
,
D.G.公司。
(
1993
)
广义线性混合模型中的近似推理
.
《美国统计杂志》。助理。
,
88
,
9
25
.

克雷文
,
第页。
瓦赫巴
,
G.公司。
(
1979
)
用样条函数平滑噪声数据:用广义交叉验证方法估计平滑的正确程度
.
数字。数学。
,
31
,
377
403
.

艾尔斯
,
P.H.C.公司。
,
游戏
,
J。
,
马克思
,
出生日期。
,
对。
(
2008
)
季节时间序列和关联表的调制模型
.
统计师。医学。
,
27
,
3430
3441
.

艾尔斯
,
P.H.C.公司。
马克思
,
出生日期。
(
1996
)
使用B样条和惩罚进行灵活平滑
.
统计师。科学。
,
11
,
89
121
.

法尔迈尔
,
L。
,
科尼布
,
T。
冗长的
,
美国。
(
2004
)
时空数据的惩罚结构加性回归:一个贝叶斯的视角
.
统计师。罪。
,
14
,
731
761
.

格鲁布
,
G.H公司
货车贷款
,
成本加运费
. (
1996
)
矩阵计算
,第3版。
巴尔的摩
:
约翰霍普金斯大学出版社
.

绿色
,
P.J公司
.和
西尔弗曼
,
B.宽
. (
1994
)
非参数回归与广义线性模型
.
伦敦
:
查普曼和霍尔
.

,
C、。
(
1992
)
交叉验证非高斯数据。
J.计算图表。统计师。
,
1
,
169
179
.

,
C类
. (
2002
)
平滑样条方差分析模型
.
纽约
:
施普林格
.

,
C类
.和
基姆
,
Y.-J.公司。
(
2002
)
惩罚似然回归:一般公式和有效逼近
.
可以。J.统计师
.,
30
,
619
628
.

哈斯蒂
,
T。
提比什拉尼
,
对。
(
1986
)
广义可加模型(讨论)
.
统计师。科学。
,
1
,
297
318
.

哈斯蒂
,
T型
提比什拉尼
,
R(右)
. (
1990
)
广义可加模型
.
伦敦
:
查普曼和霍尔
.

锦葵
,
C.L公司
. (
1973
)
关于的一些评论C类第页
.
技术计量学
,
15
,
661
675
.

马克思
,
出生日期。
艾尔斯
,
P.H.公司。
(
1999
)
采样信号和曲线的广义线性回归:一种P样条方法
.
技术计量学
,
41
,
1
13
.

帕克
,
R(右)
.和
大米
,
J型
. (
1985
)
关于“非参数回归曲线拟合的样条平滑方法的某些方面”的讨论(B.W.Silverman)
.
J.R.统计学会
,
47
,
40
42
.

,
钢筋混凝土。
韦尔蒂
,
洛杉矶。
(
2004
)
NMMAPSdata包
.
R新闻
,
4
,
10
14
.

皮埃罗
,
A类
.和
古德
,
Y(Y)
. (
2011
)
基于广义可加模型的短期电力负荷预测
.英寸程序。第16届国际电力系统智能系统应用大会,第页。
410
415
。纽约:电气和电子工程师协会。

赖斯
,
体育
.和
奥格登
,
R.T公司
. (
2009
)
一类半参数线性模型的平滑参数选择
.
J.R.统计。Soc.B公司
,
71
,
505
523
.

木材
,
序号。
(
2000
)
具有多重二次惩罚的建模和平滑参数估计
.
J.R.统计。Soc.B公司
,
62
,
413
428
.

木材
,
序号。
(
2004
)
广义可加模型的稳定有效多重平滑参数估计
.
《美国统计杂志》。助理。
,
99
,
673
686
.

木材
,
序号
. (
2006
)
广义可加模型:R引言
.
博卡拉顿
:
查普曼和霍尔——CRC
.

木材
,
序号
. (
2009
)
mgcv公司
.
R包版本1.6-0
。(可从http://CRAN.R-project.org/package=mgcv

木材
,
序号。
(
2011
)
半参数广义线性模型的快速稳定限制极大似然和边缘似然估计
.
J.R.统计。Soc.B公司
,
73
,
36
.

附录A:二维码-基于建模细节

本附录提供了模型矩阵的进一步详细信息X(X)可以替换为二维码-因子R(右)用于平滑度选择和系数估计。

A.1。平滑度选择标准

除GCV外,第节中还介绍了GCV2,马尔洛的C类第页和REML标准也可用于平滑度选择。

如果ϕ已知,则尝试最小化预测误差的估计值会导致最小化Mallows(1973)C类第页(克雷文和瓦巴氏(1979)“无偏风险估计”)
C类第页(λ)=负极X(X)β^λ2+2ϕ信托收据(如果λ)
关于平滑参数。
或者我们可以治疗S公司j个作为高斯随机效应的精度矩阵,并估计λj个通过边际可能性或REML,即我们寻求λj个最小化
V(V)第页(λ)=负极X(X)β^λ2+β^λT型S公司λβ^λ2ϕ+n个负极M(M)第页2日志(2πϕ)+日志|X(X)T型X(X)+S公司λ|负极日志|S公司λ|+2
哪里|S公司λ|+是的严格正特征值的乘积S公司λ、和M(M)第页是等级缺失的(形式)程度S公司λ.
两个Mallows的C类第页REML标准可以重新表述为(f), ‖第页2R(右)来自节2如下:
C类第页(λ)=(f)负极R(右)β^λ2+第页2+2ϕ信托收据(如果λ)
V(V)第页(λ)=(f)负极R(右)β^λ2+第页2+β^λT型S公司λβ^λ2ϕ+n个负极M(M)第页2日志(2πϕ)+日志|R(右)T型R(右)+S公司λ|负极日志|S公司λ|+2.

A.2、。备选方案

第节中建议的方法的两个备选方案2值得评论。

  • (a)

    原则上,我们可以根据工作线性模型进行所有推断(f)=R(右)β+ɛ,其中ɛN个(0,ϕ),在这种情况下‖第页2不会再起作用了。如果ϕ是已知的,则推理对这种变化是不变的,但如果ϕ未知,然后掉落‖第页2通常会导致较差的结果,因为‖第页2包含关于ϕ.

  • (b)

    避免高存储成本的一个表面上有吸引力的替代方案是为组件使用稀疏表示(f)j个,所以X(X)是一个稀疏矩阵,产生较低的存储成本(例如,通过使用Bates和Maechler的矩阵库(2013))。然而,在测试中,与稀疏计算相关的开销意味着我们无法在速度或存储要求方面为这种方法产生实际优势,并且它大大限制了可以使用的平滑器的范围。无论如何X(X)不结转至R(右)(R(右)也是X(X)T型X(X),当有多个时,通常接近稠密(f)j个).

附录B:获取R(右),(f)和‖第页2不成形X(X)

对于大型数据集和相当灵活的模型,X(X)可能变得太大而无法放入计算机内存。事实上,任何内存占用是X(X)对于较小的数据集大小,将耗尽内存。例如木材(2011)所需的存储n个×第页每个平滑参数的矩阵λj个.来自节2很明显,如果R(右),(f)和‖第页2无需形成即可获得X(X)一气呵成。有两种可能的方法。

B.1、。二维码-更新

考虑构造一个二维码-分区的分解X(X)。假设
X(X)=(X(X)0X(X)1),
和类似的
=(01).
X(X)00两者都有n个0行,而X(X)11两者都有n个1排。n个0+n个1=n个.现在形成二维码-分解X(X)0=0R(右)0
(R(右)0X(X)1)=1R(右).
这是例行检查X(X)=二维码哪里
=(000)1
(n个1×n个1此处)和
T型=1T型(0T型01).

重复使用这种结构可以R(右)(f)只考虑一个子块X(X)一次。有了足够的块,内存占用可以减少到所需内存的一小部分,如果X(X)形成整体。如果X(X)n个×第页而且有M(M)很容易看出,这种方法的运算次数是O(运行)(净现值2+百万英镑)与…相反O(运行)(净现值2)用于显式形成X(X),即何时n个大大大于第页间接费用较小。‖第页2=‖2−‖(f)2提供估计所需的其余成分。各种类型的二维码-更新被广泛使用,包括在普通回归建模中:参见Golub和van Loan(1996),节12.5,用于讨论二维码-更新方法。

这种简单方法的一个优点是,大多数工作“令人尴尬地并行”。数据可以在处理器,每个处理器累加R(右)(f)用于其数据子集。最后,按要求R(右)(f)由一个进一步的QR分解得出:
(R(右)1R(右))=˜R(右),(f)=˜T型((f)1(f)).
(7)

B.2节。Choleski和X(X)T型X(X)-更新

如果X(X)=二维码那么很明显X(X)T型X(X)=R(右)T型R(右),所以R(右)是Choleski因子X(X)T型X(X).如果X(X)按行划分为块X(X)1,X(X)2,……然后
X(X)T型X(X)=j个X(X)j个T型X(X)j个
可以用来积累X(X)T型X(X)无需形成X(X)整体。一次X(X)T型X(X)可用,R(右)可以通过Choleski分解计算,而(f)=R(右)−TX(X)T型和之前一样,‖第页2=‖2−‖(f)2.X(X)T型可以同时累加X(X)T型X(X)同样,大多数过程都是令人尴尬的并行过程。这种方法的优点是,它大约将领先订单操作数减半,与二维码-更新。缺点是数值稳定性降低,相对于二维码-更新(例如参见木材(2004))以及在等级不足的情况下进行正规化的必要性X(X).

B.3节。在线更新

对于用于预测的大型模型,如Section的电网负荷预测模型4,的二维码-更新方法提供了一种在新数据可用时更新模型的明显方法。如果新数据不强制更改用于表示平滑项的任何基数,则以下过程是可能的。

  • 步骤1:使用新数据进行更新R(右),(f)和‖第页2,正如上面在处理模型矩阵的新计算的块时所示。

  • 第2步:重新估计平滑参数λ和系数βλ基于第1步中更新的项的平滑度选择标准,并使用平滑参数和ϕ-原始拟合的估计值作为此优化的起始值。

这是一篇根据http://creativecommons.org/licenses/by/4.0/许可证,允许在任何媒体上使用、分发和复制,前提是正确引用了原始作品。