物质扩散模型

September 27, 2015 zr9558 Leave a comment

引言

在推荐系统中，常常用二部图来表示用户和物品的关系：把用户（Users）看成一类点，把物品（Objects）看成另一类点。如果用户购买了某种物品，那么用户和物品之间就存在一条连线，也就是一条边。如果用户没有购买某种物品，那么用户和物品之间就没有连线。但是用户和用户之间，物品和物品之间是不存在连线的。换句话说就是同一类点之间是不存在连线的。这样的图结构就叫做二部图。电子商务中间的物品推荐，可以看成是二部图的边的挖掘问题。扩散过程可以用来寻找二部图中两个节点的关联强度。经典的扩散过程分成两类：一种叫做物质扩散（Probabilistic Spreading），它满足能量守恒定律；另一种叫做热传导（Heat Spreading），一般由一个或多个热源（物品）驱动，不一定满足能量守恒定律。

物质扩散（Probabilistic Spreading）算法流程

如图所示，用黄色的点来表示用户，用红色的点表示物品。用户有三位，物品有四个。第一个用户和第一个物品有连线表示第一个用户购买了第一个物品，而与第二个物品没有连线表示第一个用户没有购买第二个物品。根据图像(a)所示，第一个物品和第三个物品有能量1，第二个物品和第四个物品有能量0.对于物质扩散过程，第一个物品需要把自己的能量1平均分配给购买过它的用户一和用户二。第三个物品也需要把自己的能量1平均分配给购买过它的用户一和用户三。用户的能量就是所有用物品得到的能量的总和。正如图像(b)所示，第一个用户此时具有能量1，第二个和第三个用户都具有能量1/2.接下来，每个用户再把自身的能量平均分配给所有购买过的物品，物品的能量则是从所有用户得到的能量的总和。所以第一个物品的能量就是第一个用户的能量乘以0.5加上第二个用户的能量乘以1/3，所得的能量则是2/3.其他物品的能量用类似的方法计算。

物质扩散（Probabilistic Spreading）数学模型

假设有m个用户（Users）和n个物品（Objects），可以构造矩阵 $A_{m\times n},$ 如果 $a_{ij}=1$ ，那么表示用户i购买了物品j；如果 $a_{ij}=0$ , 那么表示用户i没有购买物品j。用 $(c_{1},\cdot\cdot\cdot, c_{n})$ 来表示物品1到物品n的能量。用 $k(U_{\alpha})=\sum_{j=1}^{n}a_{\alpha j}$ 表示用户 $\alpha$ 的度，也就是该用户购买的物品数量。用 $k(O_{j})=\sum_{i=1}^{m}a_{ij}$ 表示物品j的度，也就是该物品被多少个用户购买过。根据物质扩散的算法描述，第一步需要计算出用户从物品那里得到的能量，此时用户 $i$ 得到的能量用 $b_{i}$ 表示，那么

$b_{i}=\sum_{\ell=1}^{n} a_{i\ell}\cdot c_{\ell}/k(O_{\ell}) \text{ for all } 1\leq i\leq m.$

后面需要通过用户的能量来更新物品的能量，假设物品 $j$ 的新能量用 $c_{j}^{'}$ 表示，那么

$c_{j}^{'}=\sum_{i=1}^{m}a_{ij}\cdot b_{i}/ k(U_{i}) \text{ for all } 1\leq j \leq n.$

把 $b_{i}$ 的公式带入上面可以得到对于 $1\leq j\leq n,$ 我们有

$c_{j}^{'}=\sum_{\ell=1}^{n}\left(\frac{1}{k(O_{\ell})}\cdot\sum_{i=1}^{m}\frac{a_{ij}a_{i\ell}}{k(U_{i})}\right)\cdot c_{\ell}=\sum_{\ell=1}^{n}w_{j\ell}^{P}\cdot c_{\ell}.$

在这里对于 $1\leq j\leq n, 1\leq \ell \leq n,$

$w_{j\ell}^{P}=\frac{1}{k(O_{\ell})}\cdot\sum_{i=1}^{m}\frac{a_{ij}a_{i\ell}}{k(U_{i})}.$

假设 $n\times n$ 矩阵 $W^{P}=(w_{j\ell}^{P}),$ 列向量 $\vec{C}=(c_{1},\cdot\cdot\cdot,c_{n})^{T}$ 和 $\vec{C^{'}}=(c_{1}^{'},\cdot\cdot\cdot,c_{n}^{'})^{T}$ ，那么 $\vec{C^{'}}=W^{P}\vec{C}.$

性质1. 矩阵每列的和都是1，i.e. 对于所有的都成立。

证明. $\sum_{j=1}^{n}w_{j\ell}^{P}=\sum_{j=1}^{n}\sum_{i=1}^{m}\frac{1}{k(O_{\ell})}\cdot\frac{a_{ij}a_{i\ell}}{k(U_{i})}=\sum_{i=1}^{m}\frac{a_{i\ell}}{k(O_{\ell})}=1.$

性质2. 对于所有的都成立。换言之，矩阵不一定是对称矩阵，只有在所有物品的度都完全一致的情况下才是对称矩阵。

证明. 根据 $w_{j\ell}^{P}$ 和 $w_{\ell j}^{P}$ 的定义可以得到结论。

性质3. 矩阵对角线上的元素是该列的最大值。i.e. 对于所有的都成立。

证明. 根据矩阵 $A$ 的定义，对于所有的 $1\leq j,\ell\leq n,$ 可以得到

$w_{jj}^{P}=\frac{1}{k(O_{j})}\cdot \sum_{i=1}^{m}\frac{a_{ij}^{2}}{k(U_{i})}\geq\frac{1}{k(O_{j})}\cdot\sum_{i=1}^{m}\frac{a_{ij}a_{i\ell}}{k(U_{i})}=w_{\ell j}^{P}.$

性质4. 能量守恒定律：

证明. 根据前面的性质，已经知道 $W^{P}$ 矩阵每一列的和都是1，i.e. $\sum_{j=1}^{n}w_{j\ell}^{P}=1.$ 可以直接计算得到： $\sum_{j=1}^{n}c_{j}^{'}=\sum_{j=1}^{n}\sum_{\ell=1}^{m}w_{j\ell}^{P}c_{\ell}=\sum_{\ell=1}^{n}c_{\ell}\cdot\sum_{j=1}^{n}w_{j\ell}^{P}=\sum_{\ell=1}^{m}c_{\ell}.$

性质5. 如果物品的用户都是物品的用户，换句话说，如果 则有那么对于所有的都成立。反之，如果那么物品的用户都是物品的用户。

证明. 根据 $w_{j\ell}^{P}$ 的定义，可以得到

$w_{jj}^{P}=\frac{1}{k(O_{j})}\cdot \sum_{i=1}^{m}\frac{a_{ij}^{2}}{k(U_{i})}=\frac{1}{k(O_{j})}\cdot \sum_{\{1\leq i\leq m, a_{ij}=1\}}\frac{1}{k(U_{i})},$

$w_{\ell j}^{P}=\frac{1}{k(O_{j})}\cdot \sum_{i=1}^{m}\frac{a_{i\ell}a_{ij}}{k(U_{i})}=\frac{1}{k(O_{j})}\cdot \sum_{1\leq i\leq m, a_{ij}=1}\frac{a_{i\ell}}{k(U_{i})}.$

根据条件，如果 $a_{ij}=1,$ 则有 $a_{i\ell}=1,$ 那么上面两个式子相等。i.e. $w_{jj}^{P}=w_{\ell j}^{P}$ 对于所有的 $1\leq \ell \leq n$ 都成立。

根据上面性质，可以定义变量

$I_{j}=\sum_{\ell=1}^{n}(w_{\ell j}^{P}/w_{jj}^{P})^{2},$

其中 $w_{\ell j}^{P}/w_{jj}^{P}$ 可以看成物品 $O_{j}$ 与物品 $O_{j}$ 的独立程度。如果 $w_{\ell j}^{P}/w_{jj}^{P}=1,$ 那么说明物品 $O_{j}$ 的用户都是物品 $O_{\ell}$ 的用户。如果 $w_{\ell j}^{P}/w_{jj}^{P}=0,$ 那么说明物品 $O_{j}$ 的用户与物品 $O_{\ell}$ 的用户则没有交集。换言之，物品 $O_{j}$ 与物品 $O_{\ell}$ 则相对独立。根据这个性质，如果 $I_{j}$ 越小，那么物品 $O_{j}$ 与其他物品就相对独立。

性质6. 对于所有的都是成立的。

证明. 根据以上性质可以得到 $w_{\ell j}^{P}\leq w_{jj}^{P}$ 对于 $1\leq \ell \leq n$ 都成立，并且 $I_{j}\geq (w_{jj}^{P}/w_{jj}^{P})^{2}=1,$ 所以上面性质成立。

基于物质扩散算法的用户相似度

回顾一下基于用户的协同过滤算法：用户 $U_{\alpha}, U_{\beta}$ 的相似度是

$S_{\alpha \beta}=(\sum_{j=1}^{n}a_{\alpha j}\cdot a_{\beta j})/\min\{k(U_{\alpha}),k(U_{\beta})\}$

或者使用余弦表达式

$S_{\alpha \beta}=(\sum_{j=1}^{n}a_{\alpha j}\cdot a_{\beta j})/\sqrt{k(U_{\alpha})\cdot k(U_{\beta})}.$

如果用户 $U_{i}$ 并没有选择物品 $O_{j},$ 那么预测分数则是

$v_{ij}=(\sum_{1\leq \ell \leq m, \ell\neq i} S_{i\ell}\cdot a_{j\ell}) / \sum_{1\leq \ell \leq m, \ell\neq i} S_{i\ell}.$

这里的求和指的把除了用户 $i$ 的其他用户做相似度的加权，看用户 $U_{\ell}$ 是否购买了物品 $O_{j}.$ 然后对于任意的用户和物品组成的对 $(i,j),$ 只要 $a_{ij}=0$ （意思是用户 $U_{i}$ 并没有购买过物品 $O_{j}$ ），根据预测分数 $v_{ij}$ 从大到小对 $1\leq j\leq n$ 排序，推荐 $v_{ij}$ 值大的物品给用户 $U_{i}$ 即可。以上就是基于用户的协同过滤算法。

下面来介绍基于物质扩散方法的算法。与上面的方法类似，第一步就是定义用户相似度 $S_{\alpha \beta}^{P}$ 如下：

$S_{\alpha \beta}^{P}=\frac{1}{k(U_{\beta})}\cdot \sum_{j=1}^{n}\frac{a_{\alpha j}\cdot a_{\beta j}}{k(O_{j})}.$

上面定义表示物品 $O_{j}$ 购买的人越多，用户的相似度 $S_{\alpha \beta}^{P}$ 就会减少。原因就是比方说像新华字典这种购买的人非常多的物品，并不能说明大多数人都是相似的，反而是购买数学分析这种物品的人相似性比较大。也就是说物品的热门程度 $k(O_{j})$ 对计算用户相似度的时候有抑制作用。此时的评分系统则是：如果用户 $U_{\alpha}$ 之前没有购买物品 $O_{i},$ 那么其预测分数则是

$v_{\alpha i}=(\sum_{1\leq \beta\leq m, \beta\neq\alpha}S_{\alpha \beta}^{P}\cdot a_{\beta i})/\sum_{1\leq \beta\leq m, \beta\neq\alpha}S_{\alpha \beta}^{P}$

对于物质扩散的算法，可以用以下简单方法进行推广。增加参数 $d>0,$ 定义用户 $U_{\alpha}$ 和用户 $U_{\beta}$ 的相似度 $S_{\alpha \beta}^{P}(d)$ 如下：

$S_{\alpha \beta}^{P}(d)=\frac{1}{k(U_{\beta})}\cdot \sum_{j=1}^{n}\frac{a_{\alpha j}\cdot a_{\beta j}}{(k(O_{j}))^{d}}.$

当 $d>1$ 时，减弱了热门物品对用户相似度的影响；当 $d\in (0,1)$ 时，增加了热门物品对用户相似度的影响。某篇论文显示基于某些数据， $d=1.9$ 是最佳的参数。

Time Series

时间序列模型之相空间重构模型

September 14, 2015 zr9558 Leave a comment

一般的时间序列主要是在时间域中进行模型的研究，而对于混沌时间序列，无论是混沌不变量的计算，混沌模型的建立和预测都是在所谓的相空间中进行，因此相空间重构就是混沌时间序列处理中非常重要的一个步骤。所谓混沌序列，可以看作是考察混沌系统所得到的一组随着时间而变化的观察值。假设时间序列是 $\{ x(i): i=1,\cdot \cdot \cdot, n\},$ 那么吸引子的结构特性就包含在这个时间序列之中。为了从时间序列中提取出更多有用的信息，1980年Packard等人提出了时间序列重构相空间的两种方法：导数重构法和坐标延迟重构法。而后者的本质则是通过一维的时间序列 $\{x(i)\}$ 的不同延迟时间 $\tau$ 来构建 $d$ 维的相空间矢量

$y(i)=(x(i),...,x(i+(d-1)\tau)), 1\leq i\leq n-(d-1)\tau.$

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

1. 相空间重构

相空间重构技术有两个关键的参数：嵌入的维数 $d$ 和延迟时间 $\tau.$ 在Takens嵌入定理中，嵌入维数和延迟时间都只是理论上证明了其存在性，并没有给出具体的表达式，而且实际应用中时间序列都是有噪声的有限序列，嵌入维数和时间延迟必须要根据实际的情况来选取合适的值。

关于嵌入维数 $d$ 和延迟时间 $\tau$ 的取值，通常有两种观点：第一种观点认为 $d$ 和 $\tau$ 是互不相关的，先求出延迟时间之后再根据它求出合适的嵌入维数。求出延迟时间 $\tau$ 比较常用的方法主要有自相关法，平均位移法，复自相关法和互信息法等，关键的地方就是使得原时间序列经过时间延迟之后可以作为独立的坐标来使用。同时寻找嵌入维数的方法主要是几何不变量方法，虚假最临近法(False Nearest Neighbors)和它改进之后的Cao方法等。第二种观点则是认为延迟时间和嵌入维数是相关的。1996年Kugiumtzis提出的时间窗长度则是综合考虑两者的重要参数。1999年，Kim等人提出了C-C方法，该方法使用关联积分同时估计出延迟时间和时间窗。

(1) 延迟时间 $\tau$ 的确定：

如果延迟时间 $\tau$ 太小，则相空间向量

$y(i)=(x(i),\cdot\cdot\cdot, x(i+(d-1)\tau), 1\leq i\leq n-(d-1)\tau$

中的两个坐标分量 $x(i+j\tau)$ 与 $x(i+(j+1)\tau)$ 在数值上非常接近，以至于无法相互区分，从而无法提供两个独立的坐标分量；但是如果延迟时间 $\tau$ 太大的话，则两个坐标分量又会出现一种完全独立的情况，混沌吸引子的轨迹在两个方向上的投影就毫无相关性可言。因此需要合适的方法来确定一个合适的延迟时间 $\tau,$ 从而在独立和相关两者之间达到一种平衡。

(1.1) 自相关系数法：

自相关函数是求延迟时间 $\tau$ 比较简单的一种方法，它的主要理念就是提取序列之间的线性相关性。对于混沌序列 $x(1), x(2),\cdot \cdot \cdot, x(n),$ 可以写出其自相关函数如下：

$R(\tau)=\frac{1}{n}\sum_{i=1}^{n-\tau}x(i)x(i+\tau).$

此时就可以使用已知的数据 $x(1),\cdot \cdot\cdot x(n)$ 来做出自相关函数 $R(\tau)$ 随着延迟时间 $\tau$ 变化的图像，当自相关函数下降到初始值 $R(0)$ 的 $1-e^{-1}$ 时，i.e. $R(\tau)=(1-e^{-1})R(0),$ 所得到的时间 $\tau$ 也就是重构相空间的延迟时间。虽然自相关函数法是一种简便有效的计算延迟时间的方法，但是它仅仅能够提取时间序列的线性相关性。根据自相关函数法可以让 $x(i), x(i+\tau)$ 以及 $x(i+\tau), x(i+2\tau)$ 之间不相关，但是 $x(i), x(i+2\tau)$ 之间的相关性可能会很强。这一点意味着这种方法并不能够有效的推广到高维的研究。而且选择下降系数 $1-e^{-1}$ 时，这一点可能有点主观性，需要根据具体的情况做适当的调整。因此再总结了自相关法的不足之后，下面介绍一种判断系统非线性关系性的一种方法：交互信息法。

(1.2) 交互信息法：

在考虑了以上方法的局限性之后，Fraser和Swinney提出了交互信息法(Mutual Information Method)。假设两个离散信息系统 $\{ s_{1},\cdot \cdot \cdot, s_{m}\}, \{q_{1},\cdot\cdot\cdot, q_{n}\}$ 所构成的系统 $S$ 和 $Q$ 。通过信息论和遍历论的知识，从两个系统中获得的信息熵分别是：

$H(S)=-\sum_{i=1}^{m}P_{S}(s_{i})log_{2}P_{S}(s_{i}),$

$H(Q)=-\sum_{j=1}^{n}P_{Q}(q_{j})log_{2}P_{Q}(q_{j}).$

其中 $P_{S}(s_{i}), P_{Q}(q_{i})$ 分别是 $S$ 和 $Q$ 中事件 $s_{i}$ 和 $q_{i}$ 的概率。交互信息的计算公式是：

$I(S,Q)=H(S)+H(Q)-H(S,Q),$

其中 $H(S,Q)=-\sum_{i=1}^{m}\sum_{j=1}^{n}P_{S,Q}(s_{i},q_{j})log_{2}P_{S,Q}(s_{i},q_{j}).$

这里的 $P_{S,Q}(s_{i},q_{j})$ 称为事件 $s_{i}$ 和事件 $q_{j}$ 的联合分布概率。交互信息标准化就是

$I(S,Q)=I(S,Q)/\sqrt{H(S)\times H(Q)}.$

现在，我们通过信息论的方法来计算延迟时间 $\tau.$ 定义 $(S,Q)=(x(i), x(i+\tau)), 1\leq i\leq n-\tau,$ 也就是 $S$ 代表时间 $x(i),$ $Q$ 代表时间 $x(i+\tau).$ 那么 $I(S,Q)$ 则是关于延迟时间 $\tau$ 的函数，可以写成 $I(\tau).$ $I(\tau)$ 的大小表示在已知系统 $S$ ，也就是 $x(i)$ 的情况下，系统 $Q$ 也就是 $x(i+\tau)$ 的确定性的大小。 $I(\tau)=0$ 表示 $x(i)$ 和 $x(i+\tau)$ 是完全不可预测的，也就是说两者完全不相关。而 $I(\tau)$ 的第一个极小值，则表示了 $x(i)$ 和 $x(i+\tau)$ 是最大可能的不相关，重构时使用 $I(\tau)$ 的第一个极小值最为最优的延迟时间 $\tau.$

交互信息法的关键在于怎么计算联合概率分布 $P_{S,Q}(s_{i},q_{j})$ 以及 $S$ 和 $Q$ 系统的概率分布 $P_{S}(s_{i})$ 和 $P_{Q}(q_{j}).$ 这里采取的方法是等间距格子法，其方法简要概述如下。

令 $(S,Q)=(x(i), x(i+\tau)), 1\leq i \leq n-\tau,$ 在 $(S,Q)$ 平面用一个矩形包含上面所有的点。将矩形 $S$ 方向上等分成 $M_{1}$ 份， $Q$ 方向等分成 $M_{2}$ 份（注： $M_{1}, M_{2}$ 取值100~200之间即可）。那么在 $S$ 方向上格子的长度就是 $\epsilon_{1},$ $Q$ 方向上格子的长度就是 $\epsilon_{2}.$ 假设 $(a,b)$ 是 $(S,Q)$ 平面矩形左下角的顶点坐标。

如果 $(i-1)\epsilon_{1}\leq s-a <i\epsilon_{1},$ 那么 $s$ 在第 $i$ 个格子中，对 $Row[i]$ 做一次记录；

如果 $(j-1)\epsilon_{2}\leq s-b <j\epsilon_{2},$ 那么 $q$ 在第 $j$ 个格子中，对 $Col[j]$ 做一次记录。 $x(i)$ 和 $x(i+\tau)$ 序列的总数都是 $n-\tau,$ 那么

$P_{S}(i)=Row[i]/(n-\tau), 1\leq i\leq M_{1},$

$P_{Q}(j)=Col[j]/(n-\tau), 1\leq j\leq M_{2}.$

如果 $(i-1)\epsilon_{1}\leq s-a<i\epsilon_{1}$ 且 $(j-1)\epsilon_{2}\leq q-b < j\epsilon_{2},$ 则 $(s,q)$ 在标号为 $(i,j)$ 的格子中，对 $Together[i][j]$ 做一次记录。那么

$P_{S,Q}(i,j)=Together[i][j]/(n-\tau)^{2}, 1\leq i\leq M_{1}, 1\leq j\leq M_{2}.$

根据以上信息论的公式，可以得到：

$H(S)=-\sum_{i=1}^{M_{1}}P_{S}(i)log_{2}P_{S}(i),$

$H(Q)=-\sum_{j=1}^{M_{2}}P_{Q}(j)log_{2}P_{Q}(j),$

$H(S,Q)=-\sum_{i=1}^{M_{1}}\sum_{j=1}^{M_{2}}P_{S,Q}(i,j)log_{2}P_{S,Q}(i,j).$

$I(S,Q)=H(S)+H(Q)-H(S,Q),$

$I(S,Q)=I(S,Q)/\sqrt{H(S)\times H(Q)}.$

交互信息曲线 $I(\tau)=I(S,Q)$ 第一次下降到极小值所对应的延迟时间 $\tau$ 则是最佳延时时间。

(2) 嵌入维数 $d$ 的确定：

(2.1) 几何不变量法：

为了确定嵌入维数 $d,$ 实际应用中通常的方法是计算吸引子的某些几何不变量（如关联维数，Lyapunov指数等）。选择好延迟时间 $\tau$ 之后逐渐增加维数 $d,$ 直到他们停止变化为止。从Takens嵌入定理分析可知，这些几何不变量具有吸引子的几何性质，当维数 $d$ 大于最小嵌入维数的时候，几何结构已经被完全打开，此时这些几何不变量与嵌入的维数无关。基于此理论，可以选择吸引子的几何不变量停止变化时的嵌入维数 $d$ 作为重构的相空间维数。

(2.2) 虚假最临近点法：

从几何的观点来看，混沌时间序列是高维相空间混沌运动的轨迹在一维空间的投影，在这个投影的过程中，混沌运动的轨迹就会被扭曲。高维相空间并不相邻的两个点投影到一维空间上有的时候就会成为相邻的两点，也就是虚假邻点。重构相空间，实际上就是从混沌时间序列中恢复混沌运动的轨迹，随着嵌入维数的增大，混沌运动的轨道就会被打开，虚假邻点就会被逐渐剔除，从而整个混沌运动的轨迹得到恢复，这个思想就是虚假最临近点法的关键。

在 $d$ 维相空间中，每一个矢量

$y_{i}(d)=(x(i),\cdot \cdot \cdot, x(i+(d-1)\tau), 1\leq i\leq n-(d-1)\tau$

都有一个欧几里德距离的最邻近点 $y_{n(i,d)}(d), (n(i,d)\neq i, 1\leq n(i,d)\leq n-(d-1)\tau)$ 其距离是

$R_{i}(d)=||y_{i}(d)-y_{n(i,d)}(d)||_{2}.$

当相空间的维数从 $d$ 变成 $d+1$ 时，这两个点的距离就会发生变化，新的距离是 $R_{i}(d+1),$ 并且

$(R_{i}(d+1))^{2}=(R_{i}(d))^{2}+||x(i+d\tau)-x(n(i,d)+d\tau)||_{2}^{2}.$

如果 $R_{i}(d+1)$ 比 $R_{i}(d)$ 大很多，那么就可以认为这是由于高维混沌吸引子中两个不相邻得点投影到低维坐标上变成相邻的两点造成的，这样的临近点是虚假的，令

$a_{1}(i,d)=||x(i+d\tau)-x(n(i,d)+d\tau)||_{2}/R_{d}(i).$

如果 $a_{1}(i,d)>R_{\tau} \in [10,50],$ 那么 $y_{n(i,d)}(d)$ 就是 $y_{i}(d)$ 的虚假最临近点。这里的 $R_{\tau}$ 是阀值。

对于实际的混沌时间序列，从嵌入维数的最小值2开始，计算虚假最临近点的比例，然后逐渐增加维数 $d,$ 直到虚假最临近点的比例小于5%或者虚假最临近点的不再随着 $d$ 的增加而减少时，可以认为混沌吸引子已经完全打开，此时的 $d$ 就是嵌入维数。在相空间重构方面，虚假最临近点法（FNN）被认为是计算嵌入维数很有效的方法之一。

(2.3) 虚假最临近点法的改进-Cao方法：

虚假最临近点法对数据的噪声比较敏感，而且实际操作中需要选择阀值 $R_{\tau},$ 具有很强的主观性。此时Cao Liangyue教授提出了改进的FNN方法，此方法计算时只需要延迟时间 $\tau$ 一个参数，用较小的数据量就可以求的嵌入维数 $d.$

假设我们有时间序列 $x(1),\cdot\cdot\cdot, x(n).$ 那么基于延迟时间 $\tau$ 所构造的向量空间就是： $y_{i}(d)=(x(i),\cdot\cdot\cdot x(i+(d-1)\tau)), i=1,2, \cdot\cdot\cdot, n-(d-1)\tau,$ 这里的 $d$ 是作为嵌入维数， $\tau$ 是作为嵌入时间。 $y_{i}(d)$ 则表示第 $i$ 个 $d$ 维重构向量。与虚假最临近点法类似，定义变量

$a(i,d)=||y_{i}(d+1)-y_{n(i,d)}(d+1)||_{\infty}/||y_{i}(d)-y_{n(i,d)}(d)||_{\infty}, i=1,2,\cdot\cdot\cdot n-d\tau,$

这里的 $||\cdot ||_{\infty}$ 是最大模范数。其中 $y_{i}(d+1)$ 是第 $i$ 个 $d+1$ 维的重构向量，i.e. $y_{i}(d+1)=(x(i),\cdot\cdot\cdot, x(i+d\tau)).$ $n(i,d) \in \{1,\cdot\cdot\cdot, n-d\tau\}$ 是使得 $y_{n(i,d)}(d)$ 在 $d$ 维相空间里面，在最大模范数下与 $y_{i}(d)$ 最近的向量，并且 $n(i,d)\neq i$ ，显然 $n(i,d)$ 由 $i$ 和 $d$ 所决定。

定义

$E(d)=(n-d\tau)^{-1}\cdot\sum_{i=1}^{n-d\tau} a(i,d).$

$E(d)$ 与延迟时间 $\tau$ 和嵌入维数 $d$ 有关。为了发现从 $d$ 到 $d+1$ ，该函数变化了多少，可以考虑 $E1(d)=E(d+1)/E(d).$ 当 $d$ 大于某个值 $d_{0}$ 时， $E1(d)$ 不再变化，那么 $d_{0}+1$ 则是我们所寻找的最小嵌入维数。除了 $E1(d),$ 还可以定义一个变量 $E2(d)$ 如下。令

$E^{*}(d)=(n-d\tau)^{-1}\cdot \sum_{i=1}^{n-d\tau}|x(i+d\tau)-x(n(i,d)+d\tau)|,$

这里的 $n(i,d)$ 和上面的定义一样，也就是使得 $y_{n(i,d)}(d)$ 与 $y_{i}(d)$ 最近的整数，并且 $n(i,d)\neq i.$ 定义 $E2(d)=E^{*}(d+1)/E^{*}(d).$ 对于随机序列，数据之间没有相关性， $E2(d)\equiv 1.$ 对于确定性的序列，数据点之间的关系依赖于嵌入维数 $d$ 的变化，总存在一些值使得 $E2(d)\neq 1.$

2. 相空间的预测

通过前面的相空间重构过程，一个混沌时间序列 $x(1),\cdot \cdot \cdot, x(n)$ 可以确定其延迟时间 $\tau$ 和嵌入维数 $d,$ 形成 $n-(d-1)\tau$ 个 $d$ 维向量:

$\vec{y}_{1}=y_{1}(d)=(x(1),\cdot\cdot\cdot, x(1+(d-1)\tau)),$

$\vec{y}_{2}=y_{2}(d)=(x(2),\cdot\cdot\cdot, x(2+(d-1)\tau)),$

$\cdot \cdot \cdot$

$\vec{y}_{i}=y_{i}(d)=(x(i),\cdot\cdot\cdot, x(i+(d-1)\tau)),$

$\cdot \cdot \cdot$

$\vec{y}_{n-(d-1)\tau}=y_{n-(d-1)\tau}(d)=(x(n-(d-1)\tau),\cdot\cdot\cdot, x(n)).$

这样我们就可以在 $d$ 维欧式空间 $\mathbb{R}^{d}$ 建立一个动力系统的模型如下：

$\vec{y}_{i+1}=F(\vec{y}_{i}),$ 其中 $F$ 是一个连续函数。

(1) 局部预测法：

假设 $N=n-(d-1)\tau,$ 根据连续函数的性质可以知道：如果 $\vec{y}_{N}$ 与 $\vec{y}_{j}$ 非常接近，那么 $\vec{y}_{N+1}$ 与 $\vec{y}_{j+1}$ 也是非常接近的，可以用 $x(j+1+(d-1)\tau)$ 作为 $x(N+1+(d-1)\tau)=x(n+1)$ 的近似值。

(1.1) 局部平均预测法：

考虑向量 $\vec{y}_{N}$ 的 $k$ 个最临近向量 $\vec{y}_{t_{1}},\cdot\cdot\cdot, \vec{y}_{t_{k}}.$ i.e. 也就是从其他的 $N-1$ 个向量中选取前 $k$ 个与 $\vec{y}_{N}$ 最临近的向量，此处用欧几里德范数 $||\cdot ||_{2}$ 或者最大模范数 $||\cdot ||_{\infty}.$ 根据局部预测法的观点，可以得到 $x(n+1)$ 的近似值：

$x(n+1)\approx k^{-1}\cdot \sum_{j=1}^{k}x(t_{j}+1+(d-1)\tau) .$

也可以引入权重的概念来计算近似值：

$x(n+1)\approx \sum_{j=1}^{k}x(t_{j}+1+(d-1)\tau)\cdot \omega(\vec{y}_{t_{j}},\vec{y}_{N}).$

其中 $\omega(\vec{y}_{t_{j}}, \vec{y}_{N})=K_{h}(||\vec{y}_{t_{j}}-\vec{y}_{N}||) / \sum_{j=1}^{k} K_{h}(||\vec{y}_{t_{j}}-\vec{y}_{N}||), 1\leq j \leq k,$

$K(x)=exp(-x^{2}/2),$

$K_{h}(x)=h^{-1}K(h^{-1}x)=h^{-1}exp(-x^{2}/(2*h^{2})).$

(1.2) 局部线性预测法：

假设 $N=n-(d-1)\tau,$ 则有 $y_{N}(d)=(x(N),\cdot\cdot\cdot, x(n)).$

局部线性预测模型为

$\hat{x}(n+1)=c_{0}x(N)+c_{1}x(N+\tau)+\cdot\cdot\cdot+c_{d-1}x(N+(d-1)\tau)+c_{d},$

其中 $c_{i}, 0\leq i\leq d$ 是待定系数。假设向量 $\vec{y}_{N}$ 的 $k$ 个最临近向量是 $\vec{y}_{t_{1}},\cdot\cdot\cdot,\vec{y}_{t_{k}}.$ 此处可以用欧几里德范数 $||\cdot ||_{2}$ 或者最大模范数 $||\cdot ||_{\infty}.$ 这里使用最小二乘法来求出 $c_{i}, 0\leq i\leq d$ 的值。也就是求 $c_{i}$ 使得

$||Ab-Y||_{2}^{2}=\sum_{j=1}^{k}(\hat{x}(t_{j}+1+(d-1)\tau)-x(t_{j}+1+(d-1)\tau))^{2}$

$=\sum_{j=1}^{k}(c_{0}x(t_{j})+c_{1}x(t_{j}+\tau)+\cdot\cdot\cdot+c_{d-1}x(t_{j}+(d-1)\tau)+c_{d}-x(t_{j}+1+(d-1)\tau))^{2}$

取得最小值，其中 $Y=(x(t_{1}+1+(d-1)\tau), x(t_{2}+1+(d-1)\tau),\cdot\cdot\cdot x(t_{k}+1+(d-1)\tau))^{T},$ $b=(c_{0},c_{1},\cdot\cdot\cdot, c_{d})^{T},$ 矩阵 $A$ 的第 $j$ 行是 $d+1$ 维向量 $(x(t_{j}),\cdot\cdot\cdot, x(t_{j}+(d-1)\tau), 1), 1\leq j\leq k.$ 根据最小二乘法可以得到待定系数向量 $b=(A^{T}A)^{-1}A^{T}Y$ 是最小二乘解。

(1.3) 局部多项式预测法：

(2) 全局预测法：

(2.1) 神经网络

(2.2) 小波网络

(2.3) 遗传算法

3. 实验数据

这里使用Lorenz模型来作为测试数据，

$dx/dt=\sigma(y-x),$

$dy/dt=x(\rho-z)-y,$

$dz/dt=xy-\beta z,$

其中 $x, y, z$ 是系统的三个坐标， $\sigma, \rho, \beta$ 是三个系数。在这里，我们取 $\rho=28, \sigma=10, \beta=8/3.$ 在解这个常微分方程的时候，使用了经典的Runge-Kutta数值方法。

测试数据1：

使用1300个点，预测未来100个点，绿色曲线表示Lorenz模型在z坐标上的实际数据，红线右侧表示开始预测，蓝色曲线表示使用相空间重构模型所预测的数据。根据统计学分析可以得到：在红色线条的右侧，蓝色曲线的点和绿色曲线的点的Normalized Root Mean Square Error<8%, Mean Absolute Percentage Error<5.3%, 相关系数>97%.

把红线附近的图像放大可以看的更加清楚：

测试数据2：

使用1800个点，预测未来100个点，绿色曲线表示Lorenz模型在x坐标上的实际数据，红线右侧表示开始预测，蓝色曲线表示使用相空间重构模型所预测的数据。根据统计学分析可以得到：在红色线条的右侧，蓝色曲线的点和绿色曲线的点的Normalized Root Mean Square Error<1%, Mean Absolute Percentage Error<2%, 相关系数>99%.

把红线附近的图像放大：

使用同样的数据量，也就是n=1800，预测未来100个点，绿色曲线表示Lorenz模型在z坐标上的实际数据，红线右侧表示开始预测，蓝色曲线表示使用相空间重构模型所预测的数据。根据统计学分析可以得到：在红色线条的右侧，蓝色曲线的点和绿色曲线的点的Normalized Root Mean Square Error<1%, Mean Absolute Percentage Error<1%, 相关系数>99%.

把红线附近的图像放大：

参考文献：

Time Series

时间序列模型之灰度模型

September 14, 2015 zr9558 Leave a comment

灰度预测法是一种对含有不确定因素的系统的预测方法，灰色系统是位于白色系统和黑色系统之间的一种系统。

白色系统指的是一个系统内部的特征是完全已知的，使用者不仅知道系统的输入-输出关系，还知道实现输入-输出的具体方式，譬如函数表达式，微分方程的变化公式，或者物理学的基本定律。比方说牛顿第二定律F=ma, 使用者只需要知道物体的质量和加速度，就可以通过牛顿第二定律求出所使用的力F的具体值。或者说当物体的质量固定之后，已知不同的加速度，就可以求出不同的力的值，力和加速度之间有明确的关系表达式。这种系统就称为白色系统。

黑色系统和白色系统就有鲜明的对比，使用者完全不清楚黑色系统的内部特征，只是知道一些输入和相应的输出。使用者需要通过这一系列的输入和输出的对应关系来判断这个黑色系统的内部特征，对于没有办法打开的黑箱，使用者需要通过输入和输出来建立模型来了解这个黑箱。比方说：高校在进行招生的时候，其实对考生是几乎不了解的，这个时候就需要通过高考或者自主招生考试来判断学生的知识掌握程度和运用知识解决问题的能力。通过测试的结果来判断学生的综合能力的大小，进而判断是不是符合高校的专业需求。

而灰度系统恰好位于白色系统和黑色系统之间，灰度系统内部一部分是已知的，但是另一部分却是未知的。比方说：经济学中的模型，通过数学或者统计学的方法，可以判断某些因素之间确实有关系，但是却不能够完全的判断整个系统的情况，对于很多未知的情况，这些模型就有一定的局限性。于是，灰度系统就需要通过判断系统各个因素之间的发展规律，进行相应的关联分析。通过原始的数列来生成规律性更强的数列，通过生成的数列来建立相应的微分方程模型，从而预测数列的未来发展趋势。灰度模型使用等时间距观测到的数据值来构造灰色预测模型，从而达到能够预测未来某一时刻的数值的目的。

尽管灰度系统的现象不够清楚，内部的结构并不为人所知，数据是杂乱的，但是却是互相关联的，有整体功效的，因而可以对它的变化过程进行预测。其中灰度模型是灰度系统理论的重要组成部分，将杂乱的原始数据整理成规律性较强的数据，然后再利用离散的灰度方程建立连续的微分模型，从而对系统的发展做出全面的观察分析，并做出长期预测。

下面来介绍下灰度模型的几种分类：

已知序列 $x^{(0)}(1),...,x^{(0)}(n)$ ，其中n是数组的长度。通过分析这组序列的内在结构，从而建立相应的微分方程模型，达到预测未来序列 $x^{(0)}(n+1), x^{(0)}(n+2),...$ 的目的。

(1) 灰度模型GM(1,1)

第一步需要通过累加形成新序列 $x^{(1)}(1),..., x^{(1)}(n)$ ，其中

$x^{(1)}(1)=x^{(0)}(1),$

$x^{(1)}(2)=x^{(0)}(1)+x^{(0)}(2), \cdot \cdot \cdot$

$x^{(1)}(n)=x^{(0)}(1)+\cdot\cdot\cdot+x^{(0)}(n).$

对上面的累加序列取均值，就可以得到序列

$z^{(1)}(1)=x^{(1)}(1),$

$z^{(1)}(2)=(x^{(1)}(1)+x^{(1)}(2))/2, \cdot \cdot \cdot$

$z^{(1)}(n)=(x^{(1)}(n-1)+x^{(1)}(n))/2.$

那么GM(1,1)灰度模型相应的微分方程就是 $\frac{dx^{(1)}}{dt}+a x^{(1)}(t)=b$ , 它的离散方程是 $x^{(0)}(k)+az^{(1)}(k)=b,$ 这里 $2\leq k\leq n.$ 其中a成为发展系数，b称为控制系数。通过最小二乘法的推导，可以得到 $(a,b)^{T}=(B^{T}B)^{-1}B^{T}Y.$ 其中B和Y都是矩阵，定义如下：

$B=\begin{pmatrix} -z^{(1)}(2) & 1 \\ -z^{(1)}(3) & 1 \\ . & . \\ -z^{(1)}(n) & 1 \end{pmatrix}$ ,

$Y= (x^{(0)}(2),\cdot \cdot \cdot, x^{(0)}(n))^{T}.$

通过常微分方程的求解，可以得到预测模型的解析表达式。这里分两种情况：

第一种情况： $a\neq 0$

$x^{(1)}_{p}(k)=(x^{(0)}(1)-\frac{b}{a})e^{-a(k-1)}+\frac{b}{a}, k \geq 1,$ 从而得到序列的预测值：

$x^{(0)}_{p}(k)=x^{(1)}(k)-x^{(1)}(k-1)=(x^{(0)}(1)-\frac{b}{a})\cdot e^{-a(k-1)}\cdot (1-e^{a}), k\geq 2,$ $x^{(0)}_{p}(1)=x^{(0)}(1).$

第二种情况： $a=0$

$x^{(1)}_{p}(k)=bk+x^{(0)}(1)-b, k\geq 1,$ 从而得到的序列的预测值：

$x^{(0)}_{p}(k)=x^{(0)}(1).$

(2) 灰度Verhulst模型

灰度Verhulst模型和前面的灰度模型GM(1,1)相似，也是计算累加生成数列，然后计算中值数列。不过Verhulst模型是建立非线性的常微分方程，所以相应的矩阵B和Y就需要做更改。微分方程是 $\frac{dx^{1}(t)}{dt}+ax^{(1)}(t)=b(x^{(1)}(t))^{2}.$ 离散方程则是 $x^{(0)}(k)+az^{(1)}(k)=b(z^{(1)}(k))^{2} \text{ for } 2\leq k\leq n.$ 通过最小二乘法的推导，可以得到 $(a,b)^{T}=(B^{T}B)^{-1}B^{T}Y.$ 其中B和Y都是矩阵，定义如下：

$B=\begin{pmatrix} -z^{(1)}(2) & (z^{(1)}(2))^{2} \\ -z^{(1)}(3) & (z^{(1)}(3))^{2} \\ . & . \\ -z^{(1)}(n) & (z^{(1)}(n))^{2} \end{pmatrix}$ ,

$Y= (x^{(0)}(2),\cdot \cdot \cdot, x^{(0)}(n))^{T}.$

通过微分方程的求解，可以得到

$x_{p}^{(1)}(k)=\frac{ax^{(0)}(1)}{bx^{(0)}(1)+(a-bx^{(0)}(1))e^{a(k-1)}} \text{ for } k\geq 1.$ 从而可以得到预测的序列 $x_{p}^{(0)}(k)=x_{p}^{(1)}(k)-x_{p}^{(1)}(k-1) \text{ for }k\geq 2,$ $x_{p}^{(0)}(1)=x^{(0)}(1).$

注：如果 $a<0$ , 那么 $\lim_{k\rightarrow \infty}x_{p}^{(1)}(k)=a/b, \lim_{k\rightarrow \infty}x_{p}^{(0)}(k)=0.$

(3) 灰度模型DGM(2,1)

这个DGM(2,1)灰度模型并不需要累加生成序列，直接使用原序列就可以进行操作。微分方程是 $\frac{d^{2}x^{(1)}}{dt^{2}}+a\frac{dx^{(1)}}{dt}=b,$ 通过最小二乘法的推导，可以得到 $(a,b)^{T}=(B^{T}B)^{-1}B^{T}Y.$ 其中B和Y都是矩阵，定义如下：

$B=\begin{pmatrix} -x^{(0)}(2) & 1 \\ -x^{(0)}(3) & 1 \\ . & . \\ -x^{(0)}(n) & 1 \end{pmatrix}$ ,

$Y= (x^{(0)}(2)-x^{(0)}(1),\cdot \cdot \cdot, x^{(0)}(n)-x^{(0)}(n-1))^{T}.$

通过微分方程的求解，可以得到： $x^{(1)}_{p}(k)=(\frac{b}{a^{2}}-\frac{x^{(0)}(1)}{a})e^{-a(k-1)}+\frac{b}{a}k+(x^{(0)}(1)-\frac{b}{a})\frac{1+a}{a} \text{ for } k\geq 1.$

求解原序列的预测值就是： $x^{(0)}_{p}(k)=(\frac{b}{a^{2}}-\frac{x^{(0)}(1)}{a})(1-e^{a})e^{-a(k-1)}+\frac{b}{a} \text{ for } k\geq 2,$ $x^{(0)}_{p}(1)=x^{(0)}(1).$

案例1：灰度模型在电力行业的应用。

某市1984~1990年用电量数据如下:

年份 1984 1985 1986 1987 1988 1989 1990

用电量/GWh 2783.20 3028.26 3290.55 3477.77 3685.02 3935.09 4210.29

这里用n=4，也就是前四个点建立模型，然后后面三个点用来测试模型的正确性。注：这里使用了一些根据上面模型进行改进的灰度模型，进行了模型的融合。 $x^{(0)}(1)=2783.20, x^{(0)}(2)=3028.26, x^{(0)}(3)=3290.55, x^{(0)}(4)=3477.77.$ 通过计算累积生成序列，可以进一步求的系数a和b的值，从而根据上面的模型进行模拟。通过模型可以算出

$x^{(0)}_{p}(1)=2783.20, x^{(0)}_{p}(2)=3042.19,$ $x^{(0)}_{p}(3)=3258.00, x^{(0)}_{p}(4)=3489.12,$ $x^{(0)}_{p}(5)=3736.64, x^{(0)}_{p}(6)=4001.72, x^{(0)}_{p}(7)=4285.60.$

通过平均绝对百分比误差(Mean Absolute Percentage Error)的定义得到：前四项的MAPE<0.005, 所有七项的MAPE<0.01. 通过标准均方根误差(Normalized root mean square error)的定义可以求出前四项的NRMSE<0.027, 所有七项的NRMSE<0.031.

作图如下：

其中绿色线条表示实际的用电量数据，蓝色线条表示用灰度模型模拟的用电量数据。红线左侧表示用来测试的前四项，红线右侧表示预测的未来三年的走势。

案例2：某油藏原始数据

年份 1972 1973 1974 1975 1976 1977 1978 1979

综合含水 31.8 39.1 43.2 48.6 49.8 53.3 58.6 61.7

这里用n=5，也就是前五个点建立模型，然后后面三个点用来测试模型的正确性。

$x^{(0)}(1)=31.8, x^{(0)}(2)=39.1, x^{(0)}(3)=43.2, x^{(0)}(4)=48.6, x^{(0)}(5)=49.8.$ 通过计算累积生成序列，可以进一步求的系数a和b的值，从而根据上面的模型进行模拟。通过模型可以算出

$x^{(0)}_{p}(1)=31.8, x^{(0)}_{p}(2)=39.72,$ $x^{(0)}_{p}(3)=43.10, x^{(0)}_{p}(4)=46.78,$ $x^{(0)}_{p}(5)=50.77, x^{(0)}_{p}(6)=55.10, x^{(0)}_{p}(7)=59.80, x^{(0)}_{p}(8)=64.89.$

通过平均绝对百分比误差(Mean Absolute Percentage Error)的定义得到：前五项的MAPE<0.015, 所有八项的MAPE<0.023. 通过标准均方根误差(Normalized root mean square error)的定义可以求出前五项的NRMSE<0.054, 所有八项的NRMSE<0.052.

作图如下：

其中绿色线条表示实际的综合含水数据，蓝色线条表示用灰度模型模拟的综合含水数据。红线左侧表示用来测试的前五项，红线右侧表示预测的未来三年的走势。

参考文献：

ZHANG RONG

Monthly Archives: September 2015

物质扩散模型

引言

物质扩散（Probabilistic Spreading）算法流程

物质扩散（Probabilistic Spreading）数学模型

基于物质扩散算法的用户相似度

时间序列模型之相空间重构模型

1. 相空间重构

(1) 延迟时间 $\tau$ 的确定：

(1.1) 自相关系数法：

(1.2) 交互信息法：

(2) 嵌入维数 $d$ 的确定：

(2.1) 几何不变量法：

(2.2) 虚假最临近点法：

(2.3) 虚假最临近点法的改进-Cao方法：

2. 相空间的预测

(1) 局部预测法：

(1.2) 局部线性预测法：

(1.3) 局部多项式预测法：

(2) 全局预测法：

(2.1) 神经网络

(2.2) 小波网络

(2.3) 遗传算法

3. 实验数据

测试数据1：

测试数据2：

时间序列模型之灰度模型

zr9558's Blog

引言

物质扩散（Probabilistic Spreading）算法流程

物质扩散（Probabilistic Spreading）数学模型

基于物质扩散算法的用户相似度

1. 相空间重构

(1) 延迟时间的确定：

(1.1) 自相关系数法：

(1.2) 交互信息法：

(2) 嵌入维数的确定：

(2.1) 几何不变量法：

(2.2) 虚假最临近点法：

(2.3) 虚假最临近点法的改进-Cao方法：

2. 相空间的预测

(1) 局部预测法：

(1.2) 局部线性预测法：

(1.3) 局部多项式预测法：

(2) 全局预测法：

(2.1) 神经网络

(2.2) 小波网络

(2.3) 遗传算法

3. 实验数据

测试数据1：

测试数据2：

zr9558's Blog

(1) 延迟时间 $\tau$ 的确定：

(2) 嵌入维数 $d$ 的确定：