[机器学习] 奇异谱分析(SSA)原理及Python实现
最近做时间序列分析的时候需要用到奇异谱分析,发现网上可以查到的资料很有限,看paper的时候发现大部分也说得有些简略,所以这里看完之后总结一下。
奇异谱分析(Singular Spectrum Analysis, SSA)是一种处理非线性时间序列数据的方法,通过对所要研究的时间序列的轨迹矩阵进行分解、重构等操作,提取出时间序列中的不同成分序列(长期趋势,季节趋势,噪声等),从而进行对时间序列进行分析或去噪并用于其他一些任务。
奇异谱分析主要包括四个步骤:嵌入——分解——分组——重构。
1. 嵌入
SSA的分析对象是有限长一维时间序列 [x1,x2,...,xN][x_1, x_2,...,x_N][x1,x2,...,xN],NNN 为序列长度。首先需要选择合适的窗口长度 LLL 将原始时间序列进行滞后排列得到轨迹矩阵:
X=[x1x2⋯xN−L+1x2x3⋯xN−L+2⋮⋮⋮xLxL+1⋯xN]\boldsymbol{X}=\left[\begin{array}{cccc}{x_{1}} & {x_{2}} & {\cdots}& {x_{N- L+1}} \\ {x_{2}} & {x_{3}} & {\cdots} & {x_{N-L+2}} \\ {\vdots} & {\vdots} & {} & {\vdots} \\ {x_{L}} & {x_{L+1}} & {\cdots} & {x_{N}}\end{array}\right]X=⎣⎢⎢⎢⎡x1x2⋮xLx2x3⋮xL+1⋯⋯⋯xN−L+1xN−L+2⋮xN⎦⎥⎥⎥⎤通常情况下取 L<N/2L<N/2L<N/2。令 K=N−L+1K =N-L+1K=N−L+1,则轨迹矩阵X\boldsymbol{X}X 为L×KL\times{K}L×K的矩阵
X=[x1x2⋯xKx2x3⋯xK+1⋮⋮⋮xLxL+1⋯xN]\boldsymbol{X}=\left[\begin{array}{cccc}{x_{1}} & {x_{2}} & {\cdots}& {x_{K}} \\ {x_{2}} & {x_{3}} & {\cdots} & {x_{K+1}} \\ {\vdots} & {\vdots} & {} & {\vdots} \\ {x_{L}} & {x_{L+1}} & {\cdots} & {x_{N}}\end{array}\right]X=⎣⎢⎢⎢⎡x1x2⋮xLx2x3⋮xL+1⋯⋯⋯xKxK+1⋮xN⎦⎥⎥⎥⎤
2. 分解
接下来对轨迹矩阵进行奇异值分解,注意,这里是对轨迹矩阵进行SVD分解。看资料的时候就是在奇异值分解这里困惑了很久,具体来说就是将 X\boldsymbol{X}X 分解为以下形式:
X=UΣVT\boldsymbol{X}=\boldsymbol{U} \boldsymbol{\Sigma} \boldsymbol{V}^{T}X=UΣVT 其中 U\boldsymbol{U}U 称为左矩阵;Σ\boldsymbol{\Sigma}Σ 仅在主对角线上有值,就是奇异值,其他元素均为零; V\boldsymbol{V}V 称为右矩阵。此外 U、V\boldsymbol{U}、\boldsymbol{V}U、V 均为单位正交阵,满足 UUT=I,VVT=I\boldsymbol{U}\boldsymbol{U}^T=\boldsymbol{I}, \boldsymbol{V}\boldsymbol{V}^T=\boldsymbol{I}UUT=I,VVT=I。
由于直接对轨迹矩阵分解比较困难,因此首先计算轨迹矩阵的协方差矩阵:
S=XXT\boldsymbol{S} = \boldsymbol{X}\boldsymbol{X}^TS=XXT 接下来对 S\boldsymbol{S}S 进行特征值分解得到特征值 λ1>λ2>⋯>λL⩾0\lambda_{1}>\lambda_{2}>\cdots>\lambda_{L} \geqslant 0λ1>λ2>⋯>λL⩾0 和对应的特征向量 U1,U2,⋯,ULU_{1}, U_{2}, \cdots, U_{L}U1,U2,⋯,UL。此时U=[U1,U2,⋯,UL]\boldsymbol{U} =[U_{1}, U_{2}, \cdots, U_{L}]U=[U1,U2,⋯,UL],λ1>λ2>⋯>λL⩾0\sqrt{\lambda_{1}}>\sqrt{\lambda_{2}}>\cdots>\sqrt{\lambda_{L}} \geqslant 0λ1>λ2>⋯>λL⩾0为原序列的奇异谱 。并且有
X=∑m=1LλmUmVmT,Vm=XTUm/λm,m=1,2,...,L\boldsymbol{X}=\sum_{m=1}^{L} \sqrt{\lambda_{m}} U_{m} V_{m}^{T}, \quad V_{m}=\boldsymbol{X}^{\mathrm{T}} U_{m} / \sqrt{\lambda_{m}}, \quad m=1,2,...,L X=m=1∑LλmUmVmT,Vm=XTUm/λm,m=1,2,...,L 这里 λi\lambda_{i}λi 对应的特征向量 UiU_{i}Ui 反映了时间序列的演变型,称为时间经验正交函数(T-EOF)。
实际上python已经提供了奇异值分解的函数np.linalg.svd()
可以很方便的计算。关于奇异值分解更详细的介绍可以看这篇博客。
3. 分组
关于分组,文献中很常见的叙述是下面这样:
简单来说将所有的 LLL 个成分分为 ccc 个不相交的组,代表着不同的趋势成分。这样接下来选择主要的成分进行重构得到重构序列。Emmm。。。。这样介绍可真是太简洁明了导致动手实现的时候真是一脸懵。
因此在实现的时候参考了另一个版本,这里将分组和重构放到一块吧。。。。。这个版本有助于实现但是ran半天ran不清哪里是分组,被自己菜哭。。。。。。。。。。
4. 重构
所以这里接分解步。首先计算迟滞序列 XiX_iXi 在 UmU_mUm 上的投影:
aim=XiUm=∑j=1Lxi+jUm,j,0≤i≤N−La_{i}^{m}=\boldsymbol{X}_{i} U_m=\sum_{j=1}^{L} x_{i+j} U_{m,j}, \quad 0\leq{i}\leq{N-L} aim=XiUm=j=1∑Lxi+jUm,j,0≤i≤N−L XiX_iXi 表示轨迹矩阵 X\boldsymbol{X}X 的第 iii 列,aima_{i}^{m}aim 是 Xi\boldsymbol{X}_{i}Xi 所反映的时间演变型在原序列的xi+1,xi+2,…,xi+Lx_{i +1} , x_{i +2} ,…, x_{i +L}xi+1,xi+2,…,xi+L时段的权重, 称为时间主成分(TPC)。看到这里应当发现了,由aima_{i}^{m}aim 构成的矩阵实际上就是没有归一化的右矩阵, 即 λmVm\sqrt{\lambda_{m}}V_{m}λmVm !
接下来就可以通过时间经验正交函数和时间主成分来进行重建,具体重构过程如下:
xik={1i∑j=1iai−jkUk,j,1⩽i⩽L−11L∑j=1Lai−jkUk,j,L⩽i⩽N−L+11N−i+1∑j=i−N+LLai−jkEk,j,N−L+2⩽i⩽Nx_{i}^{k}=\left\{\begin{array}{l}{\frac{1}{i} \sum_{j=1}^{i} a_{i-j}^{k} U_{k, j}, \quad 1 \leqslant i \leqslant L-1} \\ \\{\frac{1}{L} \sum_{j=1}^{L} a_{i-j}^{k} U_{k, j}, \quad L \leqslant i \leqslant N-L+1} \\ \\ {\frac{1}{N-i+1} \sum_{j=i-N+L}^{L} a_{i-j}^{k} E_{k, j}, \quad N-L+2 \leqslant i \leqslant N}\end{array}\right. xik=⎩⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎧i1∑j=1iai−jkUk,j,1⩽i⩽L−1L1∑j=1Lai−jkUk,j,L⩽i⩽N−L+1N−i+11∑j=i−N+LLai−jkEk,j,N−L+2⩽i⩽N 这样,所有重构序列的和应当等于原序列,即
xi=∑k=1Lxiki=1,2⋯,Nx_{i}=\sum_{k=1}^{L} x_{i}^{k} \quad i=1,2 \cdots, N xi=k=1∑Lxiki=1,2⋯,N 通常情况下我们使用SSA只是为了提取原序列的主要成分,以去噪为例,我们只需要根据奇异值的大小选择前 k(k≤L)k(k \leq L)k(k≤L) 个贡献大的成分重构原序列即可。
python程序
#!/usr/bin/python3
# -*- coding: utf-8 -*-'''
@Date : 2019/11/11
@Author : Rezero
'''import numpy as np
import matplotlib.pyplot as pltpath = "xxxx" # 数据集路径series = np.loadtxt(path)
series = series - np.mean(series) # 中心化(非必须)# step1 嵌入
windowLen = 20 # 嵌入窗口长度
seriesLen = len(series) # 序列长度
K = seriesLen - windowLen + 1
X = np.zeros((windowLen, K))
for i in range(K):X[:, i] = series[i:i + windowLen]# step2: svd分解, U和sigma已经按升序排序
U, sigma, VT = np.linalg.svd(X, full_matrices=False)for i in range(VT.shape[0]):VT[i, :] *= sigma[i]
A = VT# 重组
rec = np.zeros((windowLen, seriesLen))
for i in range(windowLen):for j in range(windowLen-1):for m in range(j+1):rec[i, j] += A[i, j-m] * U[m, i]rec[i, j] /= (j+1)for j in range(windowLen-1, seriesLen - windowLen + 1):for m in range(windowLen):rec[i, j] += A[i, j-m] * U[m, i]rec[i, j] /= windowLenfor j in range(seriesLen - windowLen + 1, seriesLen):for m in range(j-seriesLen+windowLen, windowLen):rec[i, j] += A[i, j - m] * U[m, i]rec[i, j] /= (seriesLen - j)rrr = np.sum(rec, axis=0) # 选择重构的部分,这里选了全部plt.figure()
for i in range(10):ax = plt.subplot(5,2,i+1)ax.plot(rec[i, :])plt.figure(2)
plt.plot(series)
plt.show()
运行程序结果如下,左边是原始序列,右边是按奇异值排序的前十个成分序列,可以看到除了前几个剩余的基本都可以视为噪声序列。
如果取前五个序列重构,最后重构出的序列如下
相比原序列可以看到重构出的序列明显比原序列平滑,但是同时保持了总体的变化情况。
参考资料
https://www.cnblogs.com/endlesscoding/p/10033527.html
基于SSA的GPS坐标序列去噪及季节信号提取
[机器学习] 奇异谱分析(SSA)原理及Python实现相关推荐
- 【机器学习】Weighted LSSVM原理与Python实现:LSSVM的稀疏化改进
[机器学习]Weighted LSSVM原理与Python实现:LSSVM的稀疏化改进 一.LSSVM 1.LSSVM用于回归 2.LSSVM模型的缺点 二.WLSSVM的数学原理 三.WLSSVM的 ...
- 【机器学习】RBF神经网络原理与Python实现
[机器学习]RBF神经网络原理与Python实现 一.RBF神经网络原理 1. RBF神经网络结构与RBF神经元 2. RBF神经网络求解 2.1 正向传播:计算误差 2.2 反向传播:调整参数 二. ...
- 机器学习 | 早期停止法原理及Python实现
文章目录 1. 早期停止法 1.1 Python 实现 参考文献 相关文章: 机器学习 | 目录 机器学习 | 梯度下降原理及Python实现 1. 早期停止法 对于梯度下降这一类迭代学习的算法,还有 ...
- 机器学习笔记九——线性模型原理以及python实现案例
线性模型 1.线性模型概述 2 .广义线性模型 3.用于回归的线性模型 3.1 线性回归(又名普通最小二乘法) 3.1.1 单变量线性回归 3.1.2 多变量线性回归 3.2 岭回归(ridge re ...
- 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 ...
- 机器学习强基计划4-5:详解半朴素贝叶斯分类TAN原理(附Python实现)
目录 0 写在前面 1 条件互信息 2 最大带权生成树 3 TAN算法原理 4 Python实现 4.1 计算条件互信息 4.2 构造属性最大生成树 4.3 属性依赖关系可视化 4.4 预测 0 写在 ...
- 时间序列分析ARMA模型原理及Python statsmodels实践(上)
目录 1. 时间序列及相关基本概念 1.1. 时间序列分解 1.2. 时间平稳序列 1.3. 自相关与自相关函数(ACF) 1.4. 白噪声及Ljung-Box检验 1.4.1. 白噪声 1.4.2. ...
- 冲量(momentum)的原理与Python实现
冲量(momentum)的原理与Python实现 前言 参考:https://www.jianshu.com/p/58b3fe300ecb 梯度下降法(Gradient Descent)是机器学习中最 ...
- 监督学习 | ID3 决策树原理及Python实现
文章目录 1. 信息熵 Information Entropy 1.1 信息熵公式推导 2. 信息增益 Information Gain 2.1 信息增益最大化 2.1.1 利用离散特征进行分类 2. ...
最新文章
- js实现图片上传本地预览
- 霍尔开关YS1382检测速度 以及对 智能车竞赛节能组的影响
- ASP.NET MVC Music Store教程(1):概述和新项目
- 【Flutter】Dart 面向对象 ( 抽象类 | 抽象方法 )
- 网络工程师成长日记421-某银行技术支持
- wxWidgets:wxComboCtrl 示例
- python islower函数_python字符串是否是小写-python 字符串小写-python islower函数-python islower函数未定义-嗨客网...
- java 切换后台程序_将 Android 程序切换到后台及从后台切换到前台实现
- pascal闪电入门系列目录
- eclipse for php开发环境,eclipse for php 开发环境配置
- 永远不要因为这个工作不好而辞职、、、、
- java访问控制关键字_Java 访问控制关键字
- 区块链学习(二)以太坊私有链搭建
- C# PDF 静默打印
- 高中计算机会考vb教程,高中会考计算机vb知识点
- 网络层 详解,网络层功能,网络层协议,网络层设备。
- 设计模式实例学习-策略模式
- 【Qt5 for VS】关于出现 Qt platform plugin windows 运行错误的解决方案
- uniapp安卓端禁止截屏允许截屏
- Ubuntu 安装 php8.1