一维数据转二维的方法 1D to 2D

这篇文章好久之前写完了，一直放在草稿箱今天翻到了。之前研究了一下一维数据转二维的方法整理了一下，记录下来方便以后查阅。

格拉姆角场（GAF）

基础数学知识

点积（Dot product）

内积是两个向量之间的运算，用于衡量它们的“相似性”。

在最简单的情况（2D 空间）中，两个向量 $u$ 和 $v$ 之间的内积定义为：

也可以写为：

因此，如果 $u$ 和 $v$ 的范数为 1，我们得到：

因此，在处理单位向量时，它们的内积仅由 $u$ 和 $v$ 之间的角度 $\theta$（以弧度表示）来表征。并且，结果值位于 $[-1,1]$ 范围内。

注意：在欧几里德设置（维度 $n$）中，两个向量 $u$ 和 $v$ 的内积正式定义为：

格拉姆矩阵（Gram Matrix）

格拉矩阵是线性代数和几何中的一个有用工具。它经常用于计算一组向量的线性依赖关系。

定义： 一组 $n$ 个向量的 Gram 矩阵是由每对向量的点积定义的矩阵：

同样，假设所有二维向量的规范都是 1，可以得到：

其中 $\Phi (i,j)$ 是向量 $i$ 和 $j$ 之间的角度。

要点：为什么使用格拉姆矩阵？

格拉姆矩阵保留了时间依赖性。随着时间的增加，编码的数据从矩阵的左上角移动到右下角。

疑问：
时序数据的每个数据点并不是一个向量而是一个标量，那么Gram Matrix如何处理时序数据呢？
确实如此，GAF并不直接应用Gram Matrix。

格拉姆角场计算步骤

• 第一步：使用“最小-最大缩放器”将系列缩放到 [-1, 1]：

为了使内积不偏向具有最大值的观测值，我们需要缩放数据，这里需要使用MinMaxScaler：

其处理过程如下图所示：

StandardScaler 在此用例中不是合适的候选者，因为其输出范围和生成的内积都可能超过 [-1, 1]。

标准差标准化（StandardScale）使得经过处理的数据符合标准正态分布，即均值为0，标准差为1，其转化函数为：

其中$\mu$为所有样本数据的均值，$\sigma$为所有样本数据的标准差。

转化到这里，让我们进一步分析一下采用Gram矩阵的原因。
Gram矩阵表示了向量集合中的向量两两之间的内积，可以理解为 $GramMatrix$表示了向量之间的相关程度，而这种相关程度实际上是由向量之间的夹角来决定的。前面也提到过，时序数据点并不是向量，那么如何解决这个问题呢？
是否可以将其转化为带有角度的表示呢？答案是肯定的：极坐标。

• 第二步：将缩放时间序列转换为“极坐标”

需要考虑两个量：时间序列的值及其相应的时间戳。这两个变量将分别用角度和半径来表示。

假设时间序列由 $N$ 个时间戳 $t$ 和对应的值 $x$ 组成，则：

1. 使用 $arccos(x)$ 计算角度。它们位于 $[0,\pi ]$ 范围内。（$x$经过上一步的归一化，数值范围限定在了[-1, 1]，正好满足$arccos()$的输入条件）

2. 半径变量的计算方法是：首先，我们将区间 $[0,1]$ 分成 $N$ 等份。因此，我们得到 $N+1$ 个分界点 $0，...，1$。然后，我们舍弃 $0$，将这些点连续关联到时间序列中。

其转换公式为：

其中，$t_i\in N$代表了点 $x_i$的时间戳，$N$ 是时序数据中所包含的所有时间点的个数。

每一个时序点数据包含两个信息：一个是该数据点的规范化值 $\widetilde{x_i}$；另一个是其所在的时序位置$t_i$。上面的极坐标转换编码把这两个信息都包含了进来，且没有损失任何信息。

从数学上看是一个双射函数，也就是说此函数的自变量和因变量具有一一对应的关系，正反两方向都是如此。极轴$r_i$保留了时间上的关系；极角 ${\varphi }_i$保留了数值上的关系。

通过极坐标转换我们将时序数据中的每个时间点都转换成了一个向量，那么我们时候可以直接使用该向量计算内积来得到Gram Matrix？答案是不能的！

二维极空间中的内积有几个限制，因为每个向量的范数都根据时间依赖性进行了调整。更确切地说：

1. 两个不同观察值之间的内积将偏向于最近的一个（因为范数随着时间的推移而增加）。
2. 当计算观测值与其自身的内积时，所得范数也会有偏差。

因此，如果存在类似内积的运算，它应该仅取决于角度。

• 第三步：自定义内积计算Gram Matrix

把每个时序点转换为极坐标编码后任然没有解决如何表示Gram矩阵中的角度定义相关程度问题。

GAF定义了自己的特殊的内积：$⟨x_1,x_2⟩=\cos ({\varphi }_1+{\varphi }_2)$ ，就是说两个时序点之间的"内积"是这两个时序点的极坐标转化后的极角之和的余弦，形如：

上面这个带有自定义内积的矩阵就是GAF的结果。我们可以看到这是一个$n\times n$的方阵，时间被编码成方阵的几何维度 $1\cdots n$ 。

当时间滚动时，GAF把时序数据的每一个点转化成了这个点与其他点的相关关系，这种关系使用GAF自定义的内积来度量的。

再进一步拆解一下:

所以，格拉姆角场可以理解为带惩罚项的内积，请见下式:

补充一下，在这种内积定义下的GAF称为GASF（Gramian Angular Summation Fields）；而另外一种称为GADF（Gramian Angular Difference Fields），其内积定义为:

效率问题： 对于长度为$n$的序列数据，转换后的GASF或GADF尺寸为$[n,n]$的矩阵，可以采用PAA(分段聚合近似)先将序列长度减小，然后在转换。
所谓的PAA就是：将序列分段，然后通过平均将每个段内的子序列压缩为一个数值。

格拉姆角场计算流程图解

下图为格拉姆角场的计算流程：

下图为两种格拉姆角场的计算流程：

GAF的计算为什么不直接使用Gram Matrix

或者说GAF的计算Gram Matrix为什么不直接使用向量的内积？下面进行说明：

现在时间序列已经缩放，我们可以计算成对的点积并将它们存储在 Gram 矩阵中：

让我们通过检查 Gram 的值来深入了解此成像：

我们观察到两件事：

1. 输出似乎遵循以 0 为中心的高斯分布。
2. 生成的图像有噪声。

前者解释了后者，因为数据的高斯性越高，就越难将其与高斯噪声区分开来。

高斯分布并不令人惊讶。当查看内积值 $z$ 的 3D 图时，对于 $(x,y)\in R^2$ 的每种可能组合，我们得到：

假设时间序列的值遵循均匀分布 [-1, 1]，则 Gram 矩阵值遵循类高斯分布。以下是不同时间序列长度 $n$ 的 Gram 矩阵输出的直方图：

自定义内积计算Gram Matrix，现在让我们绘制格拉姆角场值的密度：

从上图中我们可以看到，格拉姆角场要稀疏得多。这是由于新构造的内积操作是对应于传统内积的惩罚版本。让我们对这个惩罚的作用有一些了解。我们先看一下完整操作的3D图：

我们可以看到：

• 惩罚将平均输出转向 $-1$；
• $x$和$y$越接近0，惩罚越大；
• 对于 $x=y$，它被转换为 $-1$；
• 输出很容易与高斯噪声区分开。

GASF 和 GADF 的区别

GASF是一种对称矩阵，‌其定义基于角度的和，‌能够保留时间序列中的周期性信息。‌GASF通过将时间序列转换为极坐标表示，‌其中时间步长转换为半径，‌序列中的值转换为角度，‌从而将一维时间序列转换为二维特征图。‌这种转换不仅保持了原始时间序列的双射性质，‌即可以通过角度和半径重构原始序列，‌而且通过半径的坐标保存了原始序列中的时间依赖性。‌GASF的对称性质使其在处理对称性问题时具有优势。‌

GADF则是一种非对称矩阵，‌其定义基于角度的差，‌能够捕捉时间序列中的变化模式。‌与GASF不同，‌GADF不具有对称性，‌它通过计算相邻数据点之间的角度差来反映时间序列中数据点之间的关系变化。‌这种差异编码方式对于检测时间序列中的突变或趋势变化非常有效。‌

具体来说，GASF强调的是时间序列中信号的整体变化趋势，而GADF则更关注信号间的局部差异和变化。因此，在实际应用中，可以根据具体问题和数据集的特点选择合适的方法。

Python应用示例

API: pyts.image.GramianAngularField

import numpy as np
import matplotlib.pyplot as plt
from pyts.datasets import load_gunpoint
from pyts.image import GramianAngularField
from sklearn.preprocessing import MinMaxScaler
from mpl_toolkits.axes_grid1 import make_axes_locatable

X, _, _, _ = load_gunpoint(return_X_y=True)
gasf = GramianAngularField(method='summation')
X_gasf = gasf.transform(X)
gadf = GramianAngularField(method='difference')
X_gadf = gadf.transform(X)

scaler = MinMaxScaler()

# 直线图及极坐标图
plt.figure(figsize=(12, 4))
plt.suptitle('gunpoint_index_' + str(0))

ax1 = plt.subplot(121)
ax1.plot(np.arange(len(X[0])), scaler.transform(np.expand_dims(X[0], axis=1)))
plt.title('rescaled time series')

ax2 = plt.subplot(122, polar=True)
r = np.array(range(1, len(X[0]) + 1)) / len(X[0])
theta = np.arccos(scaler.transform(np.expand_dims(X[0], axis=1))) * 2 * np.pi  # radian -> Angle
ax2.plot(theta, r, color='r', linewidth=3)

plt.title('polar system')
plt.show()

# GASF 和 GADF的热力图
plt.figure(figsize=(12, 4))
plt.suptitle('gunpoint_index_' + str(0))

ax1 = plt.subplot(121)
plt.imshow(X_gasf[0])
plt.title('GASF')
divider = make_axes_locatable(ax1)
cax = divider.append_axes("right", size="5%", pad=0.2) # Create an axes at the given *position*=right with the same height (or width) of the main axes
plt.colorbar(cax=cax)

ax2 = plt.subplot(122)
plt.imshow(X_gadf[0])
plt.title('GASF')
divider = make_axes_locatable(ax2)
cax = divider.append_axes("right", size="5%", pad=0.2)  # Create an axes at the given *position*=right with the same height (or width) of the main axes

plt.colorbar(cax=cax)
plt.show()

马尔可夫转换场（MTF）

基础数学知识

马尔科夫链

马尔科夫链（Markov Chains）以俄罗斯数学家安德烈-马尔科夫（Andrey Markov）命名，他是 20 世纪初研究马尔科夫链的先驱。

马尔科夫链的一个决定性特征是马尔科夫特性：过渡到任何未来状态的概率只取决于当前状态，而不取决于之前的状态序列。这一特性减少了对历史数据的依赖，从而简化了对随机过程的分析。

从数学上讲，如果我们将 $S$ 视为一组状态，将 $P$ 视为从状态 $i$ 到状态 $j$ 的转移概率，则对于任何给定的状态 $i,j\in S$，转移概率 $P_{ij}$ 决定从状态 $i$ 转移到状态 $j$ 的可能性。

马尔科夫链的核心三要素：

马尔科夫转移矩阵

定义：（Markov Status Transition Matrix）
假设系统有$N$个离散状态，状态转移矩阵 $M$是一个$N\times N$的矩阵，其中每个元素$M_{ij}$表示从状态$i$转移到状态$j$的转移概率。状态转移矩阵满足以下性质：
（1）非负性：对于所有的$$i和$j$，$M_{ij}>0$；
（2）行和为 $1$：状态转移矩阵每一行的元素之和等于$1$，即对于所有的$i$，${\sum }_{j=1}^N{M}_{ij}\ge 0$。如果我们用$Q$分位数来定义这个矩阵，那么其维度就是$Q\times Q$的。

马尔可夫转换场计算步骤

基于上述定义，我们可以实现MTF。

Step 1:
假设输入一维时间序列信号为$X=(x_t,t=1,2,\cdots ,T)$，将$X$分成$Q$个分位数段，分别标记为$1,2,\cdots ,Q$，每个分位数段内的数据量相同。将$X$中每一个数据更改为其对应的分位数段的序号。

Step 2:
构造马尔可夫状态转移矩阵$W$，其$i$行$j$列的元素$w_{ij}$表示分位数段$i$转移到分位数段$j$的频率：

通常情况下，我们会采用最大似然法来估计转移概率，简单来说$w_{ij}$可以用从状态$i$到$j$的计数除以状态$i$的次数或者是计数矩阵进行规范化。

可以看到，从原时序数据转化来的马尔可夫状态转移矩阵对于原数据的分布不太敏感，并且丢失了时间信息，这并不是一个好事。所以MTF浮出水面。

Step 3:
构造马尔可夫转换场，记作$M$ ,是一个$N\times N$矩阵，$N$为时序长度：
$$M_{kl}=w_{q_k{q}_l}$$

其中，$q_k$是$x_k$的分位桶，$q_l$ 是$x_l$的分位桶，$x$为时序数据。

​MTF 形如:

注意这里的$M$和$A$是不一样的 ,$A$中的下标是状态，而$M$中的下标是时序数据中的时间。相对于$A$的意义来讲，$M_{kl}$是$x_k$所在的分位桶转移到$x_l$所在的分位桶的概率。

MTF表示了时序数据中的任意两个时间点的数据之间的关系，相对的给出了它们之间从状态上看是否经常相邻。

整体的计算流程如下图所示：

效率问题： 原因与GAF类似，为了提高效率，设法减小M的尺寸，思路与PAA类似，将M网格化，然后每个网格中的子图用其平均值替代。

Python应用示例

API: pyts.image.MarkovTransitionField

import matplotlib.pyplot as plt
from pyts.datasets import load_gunpoint
from pyts.image import MarkovTransitionField
from skimage.measure import block_reduce
from mpl_toolkits.axes_grid1 import make_axes_locatable

X, _, _, _ = load_gunpoint(return_X_y=True)
mtf = MarkovTransitionField()
fullimage = mtf.transform(X)

# 使用block_reduce进行池化
pool_size_2 = 2 
pool_size_5 = 5 
pooled_image_2 = block_reduce(fullimage[0], (pool_size_2, pool_size_2), func=np.mean)
pooled_image_5 = block_reduce(fullimage[0], (pool_size_5, pool_size_5), func=np.mean)

# 绘图
plt.figure(figsize=(15, 4))
plt.suptitle('gunpoint_index_' + str(0))

ax1 = plt.subplot(131)
plt.imshow(fullimage[0])
plt.title('full image')
divider = make_axes_locatable(ax1)
cax = divider.append_axes("right", size="5%", pad=0.2)
plt.colorbar(cax=cax)

ax2 = plt.subplot(132)
plt.imshow(pooled_image_2)
plt.title('MTF with patch average (pool=2)')
divider = make_axes_locatable(ax2)
cax = divider.append_axes("right", size="5%", pad=0.2)
plt.colorbar(cax=cax)

ax3 = plt.subplot(133)
plt.imshow(pooled_image_5)
plt.title('MTF with patch average (pool=5)')
divider = make_axes_locatable(ax3)
cax = divider.append_axes("right", size="5%", pad=0.2)
plt.colorbar(cax=cax)

plt.show()

注： 关于马尔可夫转换场使用Python代码实现的更详细的讨论请参考这篇博客。

递归图（RP）

基础数学知识

相空间

‌相空间是数学和‌物理学中的一个重要概念，用于描述物理系统的所有可能状态。它是一个假想的空间，其中每个点代表系统的一个特定状态，如粒子的位置和动量。

相空间的维数取决于系统的自由度，即系统可以改变的方式。例如，一个质点在三维空间中运动，具有三个位置坐标和三个动量坐标，因此其相空间是六维的。通过相空间，我们可以理解和预测物理系统的演化和性质。‌

相空间是一个强大的工具，能够揭开物理系统的神秘面纱，让我们从不同角度去观察和理解自然界的现象。

相空间重构

基本概念

相空间重构的基本思想为：系统中任意变量的变化过程由与之相互作用的分量共同决定，其发展过程也隐含着其他变量的变化规律，故可从变量的变化过程中构建和恢复整个系统的变化规律。

重构的目的是为了挖掘整个时间序列更多的信息(包括相关性和动力学特性等)，找到与原始系统在某种意义上等价的另一个新系统.

相空间重构的主要是坐标延迟法，除此之外还有一种是导数重构法，但导数重构中会进行微分操作，对误差敏感且需要知道一些先验信息，所以一般不采用。

坐标延迟法

坐标延迟法本质是通过给定的一维时间序列 { x ( i ) : i = 1 , ⋅ ⋅ ⋅ , n } \{x(i):i=1,\cdot\cdot\cdot,n\} { x(i):i=1,⋅⋅⋅,n}的不同时间延迟来构建$m$维相空间矢量，即

$$\vec{y}(i)=(x(i),\dots ,x(i+(d-1)\tau )),\mathrm{where}1\le i\le n-(d-1)\tau$$

如此就将第$i$个序列值构建成了$d$维向量，其中$\tau$为延迟时间，且$i\in [1,n-(m-1)\tau ]$，所以$n$个这样的序列值就可以构成一个$n\times d$的矩阵相空间

$$Y=\left[\begin{array}{cccc}{x}_1& x_{1+\tau }& \cdots & x_{1+(d-1)\tau }\\ x_2& x_{2+\tau }& \cdots & x_{2+(d-1)\tau }\\ ⋮& ⋮& & ⋮\\ x_n& x_{n+\tau }& \cdots & x_{n+(d-1)\tau }\end{array}\right]$$

注：这里是一维时间序列构成的相空间

1981年 Takens 提出嵌入定理：对于无限长，无噪声的$d^{{}^{\prime }}$维混沌吸引子的一维标量时间序列 { x ( i ) : 1 ≤ i ≤ n } \{x(i):1\leq i\leq n\} { x(i):1≤i≤n} 都可以在拓扑不变的意义下找到一个$d$维的嵌入相空间，只要维数$d\ge 2d^{{}^{\prime }}+1$ 。
根据 Takens 嵌入定理，我们可以从一维混沌时间序列中重构一个与原动力系统在拓扑意义下一样的相空间，混沌时间序列的判定，分析和预测都是在这个重构的相空间中进行的，因此相空间的重构就是混沌时间序列研究的关键。

由坐标延迟法可以看出，有两个参数需要确定，一个是延迟时间 τ \tau