多维奇异谱分析(Multivariate Singular Spectrum Analysis,MSSA)
多维奇异谱分析(Multivariate Singular Spectrum Analysis,MSSA)
- 1. HMSSA
- 1. 1 分解
- 1.1.1 第一步 嵌入
- 1.1.2 第二步 SVD
- 1.2 重构
- 1.2.1 第三步 分组
- 1.2.2 第四步 对角平均
- 1.3 预测
- 1.3.1 recurrent forecast
- 1.3.2 vector forecast
- 2. VMSSA
- 2.1 分解
- 2.1.1 第一步 嵌入
- 2.1.2 第二步 SVD
- 2.2 重构
- 2.2.1 第三步 分组
- 2.2.2 第四步 对角平均
- 2.3 预测
- 2.3.1 recurrent forecast
- 2.3.2 vector forecast
在实际运用中,通常需要对多个时间序列进行分析,而每个时间序列都有内部结构,且序列之间也会存在一定的依赖关系。因此,基本的单变量奇异谱分析(SSA)可以扩展到多元情况,也就是多维奇异谱分析(MSSA)。
轨迹矩阵是求解SSA的主要工具,对于多变量情况,轨迹矩阵可以有不同的定义方式。主要可以分为水平和垂直两种方式。而每一种定义方式下,又可以有recurrent和vector两种预测方法,如下图所示:
下面就主要梳理一下MSSA的内容。和SSA一样,MSSA也由两个互补阶段组成:分解和重建。第一阶段对序列进行分解,第二阶段对无噪声序列进行重构。
1. HMSSA
1. 1 分解
1.1.1 第一步 嵌入
考虑M个时间序列,先根据SSA的方法,生成第i(i=1,...,M)i(i=1,...,M)i(i=1,...,M)个序列的轨迹矩阵X(i)X^{(i)}X(i),而后,水平地将这些轨迹矩阵拼接在一起:
X=[X(1):X(2):...:X(M)]X=[X^{(1)}:X^{(2)}:...:X^{(M)}]X=[X(1):X(2):...:X(M)]
要顺利的实现拼接,假设第iii个序列的轨迹矩阵X(i)X^{(i)}X(i)是Li∗KiL_i*K_iLi∗Ki的矩阵,其中,Li(2≤Li≤Ni−1)L_i(2\le L_i\le N_i-1)Li(2≤Li≤Ni−1)是每个长度为NiNiNi的时间序列的窗口长度,Ki=Ni−Li+1,i=1,...,MK_i=N_i-L_i+1, i=1,...,MKi=Ni−Li+1,i=1,...,M,则有L1=...=LM=L,L_1=...=L_M=L,L1=...=LM=L,而Ki,NiK_i,N_iKi,Ni可以不同。
1.1.2 第二步 SVD
这一步,对前面得到的XHX_HXH进行SVD操作,用λH1,...,λHL\lambda_{H_1},...,\lambda_{H_L}λH1,...,λHL表示XHXHTX_HX_H^TXHXHT的特征值,同样,有λH1≥,...,≥λHL≥0\lambda_{H_1} \ge ,...,\ge\lambda_{H_L}\ge0λH1≥,...,≥λHL≥0,UH1,...,UHLU_{H_1},...,U_{H_L}UH1,...,UHL则表示对应的特征向量。
对XHX_HXH进行SVD的过程可以表示为:XH=XH1+...+XHLX_H=X_{H_1}+...+X_{H_L}XH=XH1+...+XHL,其中
因此,XHXHTX_HX_H^TXHXHT的结构等价于:
1.2 重构
1.2.1 第三步 分组
与SSA类似,这个步骤相当于分解矩阵XH1,...,XHLX_{H_1},...,X_{H_L}XH1,...,XHL分成几个不相交的组,并对每组中的矩阵求和。下标集合1,…,L{1,…,L}1,…,L被划分成不相交的子集I1,……,ImI_1,……,I_mI1,……,Im,与XH=XI1+...+XImX_H=X_{I_1}+...+X_{I_m}XH=XI1+...+XIm相对应。
1.2.2 第四步 对角平均
对角平均的目的是将重构的hankel矩阵XIjX_{I_j}XIj变换为一维的序列。对于每个重构的轨迹矩阵而言,这一步与SSA相同。
1.3 预测
1.3.1 recurrent forecast
对于固定的LLL,为每个序列YNi(i)(i=1,...,M)Y^{(i)}_{N_i}(i=1,...,M)YNi(i)(i=1,...,M)分别构造轨迹矩阵X(i)=[X1(i),...,XK(i)]=(xmn)m,n=1L,KiX^{(i)}= [X^{(i)}_{1},...,X^{(i)}_{K}]=(x_{mn})^{L,K_i}_{m,n=1}X(i)=[X1(i),...,XK(i)]=(xmn)m,n=1L,Ki
构造块轨迹矩阵XVX_VXV如下:
定义UVj=(Uj(1),…,Uj(M))T,U_{V_j}=(U_{j}^{(1)},…,U_{j}^{(M)})^T,UVj=(Uj(1),…,Uj(M))T, 作为一个XVXVTX_VX_V^TXVXVT的第jjj个特征向量。其中Uj(i)U_{j}^{(i)}Uj(i)的长度为L,对应于序列YNi(i)(i=1,...,M)Y^{(i)}_{N_i}(i=1,...,M)YNi(i)(i=1,...,M)
定义X^V=[X^1:...:X^K]=Σi=1rUViUViTXV\hat{X}_V=[\hat{X}_1:...:\hat{X}_K]=\Sigma_{i=1}^{r}U_{V_i}U^T_{V_i}X_VX^V=[X^1:...:X^K]=Σi=1rUViUViTXV,即
考虑矩阵
作为从上一步中获得的矩阵X^(i)\hat{X}^{(i)}X^(i)汉克尔化过程的结果令Uj(i)∇U^{(i)\nabla}_{j}Uj(i)∇是特征向量Uj(i)U_{j}^{(i)}Uj(i)的前L−1L-1L−1个元素构成的向量。πj(i)\pi_{j}^{(i)}πj(i)表示特征向量Uj(i)U_{j}^{(i)}Uj(i)的最后一个元素。(i=1,...,M)(i=1,...,M)(i=1,...,M)
选择r个重建阶段的特征三元组用于预测。
定义矩阵U∇M=(U1∇M,...,Ur∇M)U^{\nabla M}=(U^{\nabla M}_1,...,U^{\nabla M}_r)U∇M=(U1∇M,...,Ur∇M),其中,Uj∇MU^{\nabla M}_jUj∇M为:
定义矩阵WWW:
如果矩阵(IM×M−WWT)−1(I_{M×M}−WW^T)^{−1}(IM×M−WWT)−1存在,且r≤Lsum−Mr≤L_{sum}−Mr≤Lsum−M,则存在h步VMSA预测,公式如下:
其中,y~ji(m)\tilde{y}^{(m)}_{j_i}y~ji(m)是重构序列mmm的第jij_iji个观察值,m=1,...,M.m=1,...,M.m=1,...,M.
1.3.2 vector forecast
HMSSA-V的过程非常类似于它的单变量版本SSA-V和HMSSA-R.算法的第一部分与HMSSA-R一样,继续考虑以下矩阵:
- 定义向量Zj(i),(i=1,...M)Z_j^{(i)},(i=1,...M)Zj(i),(i=1,...M)如下:
其中,X~j(i)\tilde{X}^{(i)}_{j}X~j(i)是分组并保留噪声分量后,第i个序列轨迹矩阵的重构列。 - 通过构造矩阵Z(i)=[Z1(i),...,Zki+h+L−1(i)]Z^{(i)}= [Z^{(i)}_1,...,Z^{(i)}_{ k_i+h+L−1}]Z(i)=[Z1(i),...,Zki+h+L−1(i)],并求对角平均,我们可以得到一个新的序列y^1(i),...,y^Ni+h(i)\hat{y}_1^{(i)},...,\hat{y}_{N_i+h}^{(i)}y^1(i),...,y^Ni+h(i),其中,y^N+1(i),...,y^N+h(i)\hat{y}_{N+1}^{(i)},...,\hat{y}_{N+h}^{(i)}y^N+1(i),...,y^N+h(i)则是通过HMSSA-V方法获得的hhh个预测值。
2. VMSSA
2.1 分解
2.1.1 第一步 嵌入
与HMSSA一样,考虑M个时间序列,先根据SSA的方法,生成第i(i=1,...,M)i(i=1,...,M)i(i=1,...,M)个序列的轨迹矩阵X(i)X^{(i)}X(i)
而后,垂直地将这些轨迹矩阵拼接在一起,得到XVX_VXV:
要顺利的实现拼接,假设第iii个序列的轨迹矩阵X(i)X^{(i)}X(i)是Li∗KiL_i*K_iLi∗Ki的矩阵,其中,Li(2≤Li≤Ni−1)L_i(2\le L_i\le N_i-1)Li(2≤Li≤Ni−1)是每个长度为NiNiNi的时间序列的窗口长度,Ki=Ni−Li+1,i=1,...,MK_i=N_i-L_i+1, i=1,...,MKi=Ni−Li+1,i=1,...,M,与HMSSA不同的是,VMSSA要求K1=...=KM=K,K_1=...=K_M=K,K1=...=KM=K,而Li,NiL_i,N_iLi,Ni可以不同。
2.1.2 第二步 SVD
这一步,对前面得到的XVX_VXV进行SVD操作,用λV1,...,λVL\lambda_{V_1},...,\lambda_{V_L}λV1,...,λVL表示XVXVTX_VX_V^TXVXVT的特征值,同样,有λV1≥,...,≥λVL≥0\lambda_{V_1} \ge ,...,\ge\lambda_{V_L}\ge0λV1≥,...,≥λVL≥0,UV1,...,UVLU_{V_1},...,U_{V_L}UV1,...,UVL则表示对应的特征向量。
XVXVTX_VX_V^TXVXVT的结构等价于:
2.2 重构
2.2.1 第三步 分组
与HMSSA类似,这个步骤相当于分解矩阵XV1,...,XVLX_{V_1},...,X_{V_L}XV1,...,XVL分成几个不相交的组,并对每组中的矩阵求和。下标集合1,…,L{1,…,L}1,…,L被划分成不相交的子集I1,……,ImI_1,……,I_mI1,……,Im,与XV=XI1+...+XImX_V=X_{I_1}+...+X_{I_m}XV=XI1+...+XIm相对应。
2.2.2 第四步 对角平均
对角平均的目的是将重构的hankel矩阵XIjX_{I_j}XIj变换为一维的序列。对于每个重构的轨迹矩阵而言,这一步与SSA相同。
2.3 预测
2.3.1 recurrent forecast
考虑MMM个序列YNi(i)=(y1(i),...,yNi(i))Y^{(i)}_{N_i}=(y^{(i)}_{1},...,y^{(i)}_{N_i})YNi(i)=(y1(i),...,yNi(i))和对应的窗口长度Li,2≤Li≤Ni−1,i=1...,ML_i,2≤Li≤Ni−1,i=1...,MLi,2≤Li≤Ni−1,i=1...,M,VMSSA-R预测算法如下所示。
对于固定的KKK,为每个序列YNi(i)(i=1,...,M)Y^{(i)}_{N_i}(i=1,...,M)YNi(i)(i=1,...,M)分别构造轨迹矩阵X(i)=[X1(i),...,XK(i)]=(xmn)m,n=1L,KiX^{(i)}= [X^{(i)}_{1},...,X^{(i)}_{K}]=(x_{mn})^{L,K_i}_{m,n=1}X(i)=[X1(i),...,XK(i)]=(xmn)m,n=1L,Ki
构造块轨迹矩阵XHX_HXH如下:
定义UHj=(u1j,…,uLj)T,U_{H_j}=(u_{1j},…,u_{Lj})^T,UHj=(u1j,…,uLj)T, 长度为L,作为一个XHXHTX_HX_H^TXHXHT的第jjj个特征向量。
定义Hankel算子Xi^\hat{X_i}Xi^,考虑由r个分量获得的重构矩阵XH^=Σi=1rUHiUHiTXH\hat{X_H}=\Sigma_{i=1}^{r}U_{H_i}U^T_{H_i}X_HXH^=Σi=1rUHiUHiTXH
考虑矩阵
作为从上一步中获得的矩阵X^(i)\hat{X}^{(i)}X^(i)汉克尔化过程的结果令UHj∇U^\nabla_{H_j}UHj∇是特征向量UHjU_{H_j}UHj的前L−1L-1L−1个元素构成的向量。πHj\pi_{H_j}πHj表示特征向量UHjU_{H_j}UHj的最后一个元素。(j=1,...,r)(j=1,...,r)(j=1,...,r)
定义ν2=Σj=1rπHj2\nu^2=\Sigma_{j=1}^{r}\pi_{H_j}^2ν2=Σj=1rπHj2
定义线性系数向量:
如果ν2<1\nu^2<1ν2<1,则存在HMSSA的H步预测,计算公式如下:
其中,i=1,...,M.i=1,...,M.i=1,...,M.
2.3.2 vector forecast
此算法的前10部分与VMSSA-R相同,继续考虑矩阵:
其中,
可以证明,矩阵Π\PiΠ是U∇U^{\nabla}U∇的列空间上的正交投影算子。令
其中,Π(i)\Pi^{(i)}Π(i)的维度为(Li−1)×(Lsum−M)(L_i−1)×(L_{sum}−M)(Li−1)×(Lsum−M),而R(i)(i=1,...,M)\mathscr{R}^{(i)}(i=1,...,M)R(i)(i=1,...,M)的长度为Lsum−ML_{sum}−MLsum−M,对应于序列YNi(i)Y_{N_i}^{(i)}YNi(i).然后,定义线性运算:
其中,YΔT=(YΔ(1),...,YΔ(M))Y_\Delta^T=(Y_\Delta^{(1)},...,Y_\Delta^{(M)})YΔT=(YΔ(1),...,YΔ(M)),而YΔ(i),(i=1,...,M)Y_\Delta^{(i)},(i=1,...,M)YΔ(i),(i=1,...,M)表示YiY_iYi(长度为LiL_iLi)的最后Li−1L_i-1Li−1个元素。使用上面的符号,VMSSA-V预测算法如下:
- 定义向量ZiZ_iZi如下:
其中,X~i\tilde{X}_{i}X~i是分组并保留噪声分量后,第i个序列轨迹矩阵的重构列。Lmax=max{L1,...,LM}L_{max}=max\{L_1,...,L_M\}Lmax=max{L1,...,LM}
2. 构造矩阵Z=[Z1,...,ZK+h+Lmax−1]Z= [Z_1,...,Z_{ K+h+L_{max}−1}]Z=[Z1,...,ZK+h+Lmax−1],并求对角平均,我们可以得到一个新的序列y^1(i),...,y^N+h+Lmax(i)(i=1,...,M)\hat{y}_1^{(i)},...,\hat{y}_{N+h+L_{max}}^{(i)}(i=1,...,M)y^1(i),...,y^N+h+Lmax(i)(i=1,...,M)
3. 元素y^Ni+1(i),...,y^Ni+h(i)(i=1,...,M)\hat{y}_{N_i+1}^{(i)},...,\hat{y}_{N_i+h}^{(i)}(i=1,...,M)y^Ni+1(i),...,y^Ni+h(i)(i=1,...,M)则是通过VMSSA-V方法获得的hhh个预测值。
多维奇异谱分析(Multivariate Singular Spectrum Analysis,MSSA)相关推荐
- matlab 谱分析函数,科学网—经典谱分析(Power Spectrum Analysis) - 刘磊的博文
频谱分析(或功率谱分析)大家可能都不陌生,然而细究起来,恐怕还是有很多模糊的地方.博主很早就知道并且学会使用各种数学软件来计算功率谱,可是长期以来总是知其然而不知其所以然,一些道理和概念总是含含糊糊过 ...
- [机器学习] 奇异谱分析(SSA)原理及Python实现
最近做时间序列分析的时候需要用到奇异谱分析,发现网上可以查到的资料很有限,看paper的时候发现大部分也说得有些简略,所以这里看完之后总结一下. 奇异谱分析(Singular Spectrum A ...
- matlab 奇异谱分析,SSA寻北-奇异谱分析
% function [y,r,vr]=ssa(x1,L) %使用方法:[y,r,vr]=ssa(x1,L) clear all clc close all t=0:0.1:1000; x1=[sin ...
- 微软服务器模式表格多维,在表格模式下安装 Analysis Services
在表格模式下安装 Analysis Services 05/08/2013 本文内容 如果您要安装 Analysis Services 以便使用新增的表格建模功能,则必须在支持这种模型的服务器模式下安 ...
- KDD2021|小红书在推荐多样化的实践——SSD
|company:小红书: |stage:Re-rank: |model:SDD: 前言 本次要分享的是小红书在2021KDD上发表的工业推荐论文--""Sliding Spect ...
- 非平稳信号的自适应分解算法:EMD、SSA、ITD、VMD以及其变体之间的总结与对比
目录 一.EMD 1.算法步骤: 2.算法优点: 3.算法缺点及其解决方法 4.EMD算法的变体:主要针对EMD的模态混叠问题提出的 二.SSA 1.算法步骤: 2.算法优点: 3.算法缺点及其解决方 ...
- MATLAB不同时频信号处理方法介绍及效果对比
本文欢迎非商业目的的学习分享转载,转载请附上原文链接及作者ID 本文为作者自身的一个学习总结,大部分内容在相关教材上也可以找到,有空的也会不定期更新.本身也在学习的过程中,出现错误在所难免,欢迎大家在 ...
- 文献综述--------山东某地区基于深度学习神经网络的配电网负荷预测研究
摘 要:地区电网负荷预测是供电企业在电网建设.运营过程中一项十分要的基础性的工作.小到一个企业的负荷预测,大到全国性电网的负荷预测研究,它的应用结果都会对适用范围内的企业经营管理.电力设施(电网)的 ...
- ML.NET 1.1 发布,模型构建器升级和新的异常检测算法
https://www.toutiao.com/a6703297861997560328/ 2019-06-17 08:59:15 ML.NET 1.1 已发布.ML.NET 是一个跨平台的机器学习框 ...
- 信号处理:希尔伯特黄变换
目录: 目录: 前言 简介 基本原理 经验模态分解 希尔伯特变换 特点 (1)HHT能分析非线性非平稳信号. (2)HHT具有完全自适应性. (3)HHT不受Heisenberg测不准原理制约--适合 ...
最新文章
- 实战1--应用EL表达式访问JavaBean的属性
- 数据库schema 是什么
- boost::math::non_central_chi_squared用法的测试程序
- GDCM:将PDF文件转换为DICOM / PDF文件的测试程序
- CF886E Maximum Element(dp、组合数学)
- ubuntu下定时任务的执行
- 轻量级的移动开发JavaScript框架-zepto.js
- 苹果手机速度慢_安卓手机用户也想体验一下MagSafe充电器?还是算了吧!
- 人工智能语音训练数据的制作方式?
- 自动控制系统中的典型环节
- 【Google Paper】对比学习用于解决推荐系统长尾问题
- marvel 1548 phy芯片调式
- java ssm人体健康体检信息管理系统-
- 怎么进行用户体验与可用性测试?
- HTML期末大作业—— 游戏网页(5个页面) ~ 全屏游戏美术大赛作品征集网页 HTML+CSS+JS ~ web课程设计网页规划与设计
- 奥特曼系列ol如何进老服务器,《奥特曼系列OL》新手攻略
- 商务办公软件应用与实践【7】
- python socket recv非阻塞_socket非阻塞recv大坑
- 课程设计:通讯录系统(数据库)
- 从光波叠加到条纹分布的matlab仿真,基于Matlab仿真算法的光源空间相干性研究
热门文章
- 如何有效的阅读开源代码
- excel使用教程_Excel教程大合集:史上最全面的Excel视频教程合集+模板,免费送...
- 宏碁传奇14 Swift 指纹模块失效解决
- 伺服速度控制模式接线图_PLC采用转矩、位置、速度模式控制伺服电机的方法
- 三菱plc控制步进电机实例_自动化工程师必掌握的PLC控制步进电机逻辑思路
- loki日志收集系统部署
- 使用MySQL创建数据库,实现基本SQL语句
- 洛杉矶马拉松本周末开跑!27000名选手正常参赛
- POS58票据热敏打印机,怎么用ESC/POS命令控制打印文字大小?
- c语言求正弦余弦正切,公式( 正弦 余弦 正切 余切 正割 余割 )