状态空间模型起源于平稳时间序列分析, 其自被提出以来受到了广泛的关注. 时间序列观测通常是一个不可直接观测动态过程的间接观测, 该不可直接观测的动态过程往往是研究者比较感兴趣的对象. 这里不可直接观测通常来源于测量过程中的污染干扰或者动态过程本身就是不可观测的. 由于其能够处理隐藏变量的内在能力, SSMs模型对于经济学和生物学研究具有非常大的吸引力^{[1, 2]}. 例如, 股票价格通常是包含噪声的一种股票真实价值的间接表现. 一般地, 状态空间模型^[3]通过两个方程来定义: 一个是状态转移方程, 其规定了随机过程随时间的变化规则; 另一个是观测方程, 其描述了潜在状态或隐藏过程如何被观测

(1.1) $ \begin{equation} \begin{array}{ll} x_t = h_t({\bf{x}}_{t-1}, \epsilon_t)\quad \mbox{或}\quad x_t\sim q_t(\cdot|{\bf{x}}_{t-1}), \\ y_t = g_t({\bf{x}}_t, \delta_t)\quad \mbox{或} \quad y_t\sim f_t(\cdot|{\bf{x}}_t), \end{array} \end{equation} $

这里$ x_t $是潜藏状态变量, $ y_t $为观测结果, 且黑体$ {\bf{x}}_j=(x_1, \cdots, x_j)' $代表了截至时间$ j $为止的状态$ x_t $变量的整个历史. 在大部分参考文献中, 潜在的状态转移过程方程$ h_t $以及噪声变量$ \epsilon_t $是分布已知的. 即, 状态变量$ x_t $通过条件分布$ q_t(\cdot) $实现转换. 一般地, 状态转移方程$ h_t $和噪声$ \epsilon_t $被假设为已知^[4-6], 除开可能的未知参外. 同样地, 观测函数$ g_t $以及观测噪声$ e_t $的分布被假设为已知, 除开可能的未知参数外. 其中, 状态转移噪声$ \epsilon_t $和观测噪声$ \delta_t $被假设为相互独立.

从通用SSM模型(1.1)可以看到, 状态转移的马尔可夫性并不是必需的, 当前状态可能依赖于状态过去的历史轨迹, 甚至整个历史过程. 类似地, 在时间节点$ t $地观测值同样可能依赖截至时间$ t $为止的整个状态轨迹. 在我们的研究中, 我们假设状态转移函数仅依赖于前一状态, 以及观测方程仅依赖于当前状态, 即

(1.2) $ \begin{equation} \begin{array}{ll} \mbox{状态转移方程}: x_t=h_t(x_{t-1}, \epsilon_t)&\quad \mbox{或}\quad x_t\sim q_t(\cdot|x_{t-1}), \\ \mbox{观测方程}: y_t = g_t(x_t, \delta_t)&\quad \mbox{或} \quad y_t\sim f_t(\cdot|x_t). \end{array} \end{equation} $

本文中称之为马尔可夫状态空间模型. 在此假设下, 我们致力于提供一种基于再生核希尔伯特空间方法的通用性函数重构方法. 其中, 转移噪声$ \epsilon_t $和观测噪声$ \delta_t $被假设为分别独立抽取自分布$ \rho^x_{x_t} $和$ \rho^o_{x_t} $. 选取马尔可夫状态空间模型作为研究重点主要基于两点考虑: 第一, 马尔可夫结构能够提供极大的简便性, 在记号、建模以及参数推断中; 第二, 可以通过变量重新标注的方式, 很容易得将有限步滞后得状态转移过程转变为单一步转台转移过程, 即马尔可夫化. 例如, 假设$ h_t(\cdot) $和$ g_t(\cdot) $仅依赖于有限步滞后变量$ x_{t-1}, \cdots, x_{t-j} $, 那么可以通过引入新的$ j $维度变量$ x^{\star}_t=(x_t, x_{t-1}, \cdots, x_{t-j+1}) $. 重新构建$ h_t, g_t $, $ x^{\star}_t=h^{\star}_t(x^{\star}_{t-1}, \epsilon_t) $和$ y_t=g^{\star}_t(x^{\star}_t, \delta_t) $, 从而将其转化为单一步马尔可夫状态空间模型. 当$ h_t, g_t $被假设具有线性或广义线性结构时, 马尔可夫状态空间模型被证明拥有较高的稳定性和极高的效率, 与之相对应的一系列估计算法被实践证明具有较高的价值, 包括卡尔曼滤波(Karman Filter)^[7], 序列门德卡洛方法(S equential Monte Carlo Methods)^[8], 独立粒子滤波(Independent Particle Filters)^[9].

本文中, 观测函数$ g_t(\cdot) $被假设为具有线性结构, 同时状态转移函数$ h_t(\cdot) $被假设为再生核希尔伯特空间$ {\cal H}_K $中的一个函数. 最后利用正则化最小二乘法实现转移函数$ h_t(\cdot) $的重构. 为此, 利用核方法估计得到的状态空间模型被称为核状态空间模型(KSSM). 其形式为

(1.3) $ \begin{equation} \begin{array}{ll} x_t= h_t(x_{t-1}) + \epsilon_t, \\ y_t=Hx_t + \delta_t, \end{array} \end{equation} $

其中, 假设$ H $为已知可逆矩阵且$ x_t, y_t\in X $来源于同一维度空间. 该假设是合理的, 因为在大部分问题中, 观测变量$ y_t $是状态变量$ x_t $的包含噪声的直接观测. 至此为止, 任务转变为如何从再生核希尔伯特空间中找到转移方程$ h_t(\cdot) $的估计, 一个一般地估计步骤是急需的.

本文的主要贡献为

● 核状态空间模型极大地增加了状态转移方程的学习能力, 从而使得马尔可夫状态空间模型能够被应用到更加复杂的场景.

● 一个正则化最小二乘法被设计用来实现核状态转移方程的估计. 我们证明了该方法的解的存在性和一致性问题, 同时给出了该方法估计下的误差估计.

● 模拟数据实验被设计证明了该算法的有效性. 通过真实数据实验, 我们发现非参数状态空间模型(核状态空间模型)能够更好得适应现实状态转移规律, 相比于传统状态空间模型在机场能见度预测数据中拥有更加良好得表现.

设$ \bar{y}=\{y_t\}_{t\in\{1, \cdots, T\}} $为观测到的数据点, $ \bar{x}=\{y_t\}_{t\in\{1, \cdots, T\}} $为状态变量, 其中状态变量不可以被直接测量. 假设KSSM (1.3)正确揭示了系统的转移规律. 这里假设$ H $已知, $ h_t(\cdot) $为带估计函数, 噪声变量$ \epsilon_t, \delta_t $分别独立抽取自分布$ \rho^x_{x_t} $和$ \rho^o_{x_t} $. 这里对任一时间状态$ x_t\in X $, 噪声变量$ \rho^x_{x_t} $和$ \rho^o_{x_t} $期望为$ 0 $, 且方差$ \sigma^{x, 2}_{x_t}, \sigma^{o, 2}_{x_t} $满足$ \sigma^{x, 2}:=\sum\limits_{t}\sigma^{x, 2}_{x_t}<\infty $, $ \sigma^{o, 2}:=\sum\limits_{t}\sigma^{o, 2}_{x_t}<\infty $. 进一步, 一个通用的函数重构方法是急需的, 其中函数来源于基于再生核希尔伯特空间$ {\cal H}_K $. 这里, 再生核希尔伯特空间被定义为:

定义2.1  假设$ {\cal X}\subset {{\Bbb R}} ^n $为闭集, $ \bar{x} $为该闭集的离散子集. 例如, $ {\cal X} $为紧集, 而$ \bar{x} $为其有限子集. 通常, 观测数据点可以被用来当作该离散子集$ \bar{x} $. 对任意给定Mercer核$ K:{\cal X}\times{\cal X}\rightarrow {{\Bbb R}} $, 基于$ \bar{x} $的特征函数记为: $ K_z(\cdot)=K(\cdot, z), z\in\bar{x} $. 在此特征函数构成的函数空间上, $ {\cal H}_{K, \bar{x}}=\{\sum\limits_{z\in\bar{x}}C_zK_z(\cdot)|C_z\in l^2_{(\bar{x})}\} $, 定义内积

(2.1) $ \begin{equation} {\nonumber} \langle K_{z_1}, K_{z_2}\rangle_K=K(z_1, z_2). \end{equation} $

称$ {\cal H}_{K, \bar{x}} $为由特征核$ \{K_z\}_{z\in\bar{x}} $生成的希尔伯特空间. 其中$ l^2_{(\bar{x})}=\{a:{\cal X}\rightarrow {{\Bbb R}} |\sum\limits_{z\in\bar{x}}a^2_z<\infty\} $为$ \bar{x} $上的可积空间.

注2.1  设$ K:{\cal X}\times {\cal X}\rightarrow {{\Bbb R}} $为连续对称半正定函数, 例如对有限点集$ \{x_1, \cdots, x_l\}\subset {\cal X} $, 矩阵$ (K(x_i, x_j))_{i, j=1}^l $为半正定矩阵. 固定第二个参数, $ K_s(\cdot)=K(\cdot, o_s) $可以看作关于$ \cdot $的一元函数, 且被挑选出的点$ o_s $唯一决定.

根据定义2.1, 观测点集$ \bar{y} $通常被挑选为唯一确定特征函数的特征点. 设$ K $为一Mercer核, 则$ {\cal H}_{K, \bar{y}}=\{\sum\limits_{z\in\bar{y}}C_zK_z(\cdot)|C_z\in l^2(\bar{y})\} $可以看作一再生核空间. 进一步, 在本文中, 观测方差可以被简化为

(2.2) $ \begin{equation} y_t=x_t+\delta_t, \end{equation} $

其中, 观测方差的矩阵$ H $为恒等变换矩阵. 这种简化主要有两个方面的好处: 第一, 有助于统计性质的提升; 第二, 一般地, 可以通过变换$ x^{\star}_t=Hx_t $将非恒等变换矩阵$ H $吸收进状态变量$ x^{\star} $. 由于$ H $是可逆, 可以通过此变换逆推出所关注的状态变量.

考虑以下核状态空间模型

(2.3) $ \begin{equation} \begin{array}{ll} x_t=f(x_{t-1}) + \epsilon_t, \\ y_t=x_t + \delta_t. \end{array} \end{equation} $

设$ {\cal H}_{K, \bar{y}}=\{\sum\limits_{z\in\bar{y}}C_zK_z(\cdot)|C_z\in l^2(\bar{y})\} $为连续函数的再生核希尔伯特空间, 其中内积按照定义2.1的方式规定. 且状态转移函数$ f(\cdot)\in{\cal H}_{K, \bar{y}} $与时间$ t $无关. 由于$ y_t=x_t+\delta_t $, 我们不妨假设$ y_t $和$ x_t $来源于同一空间$ {\cal X} $. 在给出具体的性质之前, 此处给出抽样算子的定义.

定义2.2  设$ f\in {\cal H}_{K, \bar{y}} $. 抽样算子$ S_{\bar{y}}:{\cal H}_{K, \bar{y}}\rightarrow l^2(\bar{y}) $被定义为

(2.4) $ \begin{equation} {\nonumber} S_{\bar{y}}(f)=(f(y))_{y\in \bar{y}}. \end{equation} $

一般地, 假设$ S_{\bar{y}} $是有界算子. 当$ \bar{y} $为有限集时, 抽样算子的有界性是明显的. 进一步, 这里给出算子$ S_{\bar{y}} $的共轭算子, 记为$ S^T_{\bar{y}} $. 设$ c\in l^2(\bar{y}) $, 则有

$ \langle f, S^T_{\bar{y}}c\rangle_K=\langle S_{\bar{y}}f, c\rangle_{l^2(\bar{y})}=\sum\limits_{y\in\bar{y}}c_yf(y)=\langle f, \sum\limits_{y\in\bar{y}}c_yE_y\rangle_K, \quad \forall f\in {\cal H}_{K, \bar{y}}, $

其中, 最后一个等式是由于核函数的再生性, 即对任意$ {\cal H}_{K, \bar{y}} $中的有界函数$ f:{\cal X}\rightarrow {{\Bbb R}} $, 假设最大范数至多为$ \|J\| $, 存在$ E_y\in{\cal H}_{K, \bar{y}} $且$ \|E_y\|_K\leq \|J\| $, 满足

(2.5) $ \begin{equation} f(x)=\langle f, E_x\rangle_K, \quad \forall f\in {\cal H}_{K, \bar{y}}. \end{equation} $

从而$ S^T_{\bar{y}}c=\sum\limits_{y\in\bar{y}}c_yE_y, \quad \forall c\in l^2(\bar{y}) $.

设$ \gamma\geq 0 $, 给定时间序列样本$ \bar{y}=(y_t)_{t\in\bar{t}} $, 考虑下面最优化问题

(2.6) $ \begin{equation} \tilde{f}:=\mathop{\rm arg\min}\limits_{f\in{\cal H}_{K, \bar{y}}} \left\{\sum\limits_{t\in\bar{t}}(y_t-f(y_{t-1}))^2+\gamma\|f\|^2_K\right\}. \end{equation} $

对于该函数重构问题, 以下定理给出了求解算法和解的存在性.

定理2.1  $ \tilde{f} $存在且唯一, 当$ S^T_{\bar{y}}S_{\bar{y}}+\gamma I $是可逆的, 其解的算子表达式为

$ \tilde{f}=Ly, \quad L:(S^T_{\bar{y}}S_{\bar{y}}+\gamma I)^{-1}S^T_{\bar{y}}. $

证  注意到$ f\in{\cal H}_{K, \bar{y}} $, 因此有$ f=\sum\limits_{t\in\bar{t}}c_tE_{y_t}=S^T_{\bar{y}}c $. 从而该优化问题的内积形式为

$ \mathop{\rm arg\min}\limits_ {c\in l^2(\bar{t})} \langle f - y, f-y\rangle_{l^2(\bar{y})} +\gamma\langle S_{\bar{y}}^Tc, c\rangle_{l^2(\bar{t})}. $

记

$ \begin{eqnarray*} {\cal E}_{\bar{y}}(f)&=&\langle f-y, f-y\rangle_{l^2(\bar{y})}\\ &=&\langle S^T_{\bar{y}}c-y, S^T_{\bar{y}}c-y\rangle_{l^2(\bar{y})}\\ &=&\langle S^T_{\bar{y}}c, S^T_{\bar{y}}c\rangle_{l^2(\bar{y})}-2\langle S^T_{\bar{y}}c, y\rangle_{l^2(\bar{y})} + \langle y, y\rangle_{l^2(\bar{t})}\\ &=&\langle S^T_{\bar{y}}S^T_{\bar{y}}c, c\rangle_{l^2(\bar{t})} - 2\langle S^T_{\bar{y}}c, y\rangle_{l^2(\bar{y})} + \langle y, y\rangle_{l^2(\bar{t})}\\ &=&\langle S^T_{\bar{y}}S^T_{\bar{y}}c, c\rangle_{l^2(\bar{t})} - 2\langle S^T_{\bar{y}}y, c\rangle_{l^2(\bar{t})} + \langle y, y\rangle_{l^2(\bar{t})}. \end{eqnarray*} $

最后一个等式为

$ \langle S^T_{\bar{y}}c, y\rangle_{l^2(\bar{t})} =\sum\limits_{t\in\bar{t}}y_t\sum\limits_{t\in\bar{t}}c_tE_{y_t}=\sum\limits_{t\in\bar{t}}c_tS^T_{\bar{y}}y=\langle S^T_{\bar{y}}y, c\rangle_{l^2(\bar{t})}. $

将该优化问题关于$ c $求导, 并令导数为$ 0 $, 有

$ 2S^T_{\bar{y}}S^T_{\bar{y}}c - 2S^T_{\bar{y}}y+2\gamma S^T_{\bar{y}}c=0, $

得到

$ \hat{c}=\left[S^T_{\bar{y}}S^T_{\bar{y}}+\gamma S^T_{\bar{y}}\right]^{-1}S^T_{\bar{y}}y. $

根据$ \hat{f}=S^T_{\bar{y}}\hat{c} $, 可以得到其估计为

$ \hat{f}=S^T_{\bar{y}}\hat{c} =\left[S^T_{\bar{y}}S_{\bar{y}}+\gamma I\right]^{-1}S^T_{\bar{y}}y. $

定理2.1证毕.

到此为止, 优化问题(2.6)提供了一种函数重构方法, 当$ S^T_{\bar{y}}S_{\bar{y}}+\gamma I $可逆时.

设$ f^{\star} $为真实状态转移函数, 定义

$ f^{\star}_{\bar{y}, \gamma}:=L(S_{\bar{y}}f^{\star}). $

考虑误差估计式$ \|\tilde{f}-f^{\star}_{\bar{y}, \gamma}\|^2_{K} $. 取$ \eta_t=\epsilon_t + \delta_t $, 假设其满足$ 0 $均值, 且其方差$ \sigma^2_t $满足有界方差性, 即$ \sigma^2:=\sum\limits_{t\in\bar{t}}\sigma^2_t<\infty $.

定理2.2  假设$ S^T_{\bar{y}}S_{\bar{y}}+\gamma I $可逆且关于$ \eta_t $的$ 0 $均值、有界方差性成立. 对$ \forall 0<\delta < 1 $, 以至少$ 1 - \delta $的概率, 有

$ \|\tilde{f} - f^{\star}_{\bar{y}, \gamma}\|^2_{K}\leq \frac{\|(S^T_{\bar{y}}S_{\bar{y}}+\gamma I)^{-1}\|^2\|J\|^2\sigma^2}{\delta}. $

假设$ \eta_t $有界且存在$ M $, 满足$ |\eta_t|\leq M $对$ \forall t\in\bar{t} $成立. 则对$ \forall \epsilon > 0 $, 有

$ \text{P}_y\left\{\|\tilde{f}-f^{\star}_{\bar{y}, \gamma}\|^2_{K}\leq \|L\|^2\sigma^2(1+\epsilon)\right\}\geq 1 - {\rm exp}\left\{-\frac{\epsilon\sigma^2}{2M^2}{\rm log}(1+\epsilon)\right\}. $

证  根据定理5.1, 将$ \|\tilde{f} - f^{\star}_{\bar{y}, \gamma}\|^2_{K} $展开为

$ \|L(y-S_{\bar{y}f^*})\|^2_{K}\leq \|(S^T_{\bar{y}}S_{\bar{x}} + \gamma I)^{-1}\|^2\|S^T_{\bar{y}}(y-S_{\bar{y}}f^{\star})\|^2_{K}. $

注意

$ S^T_{\bar{y}}(y - S_{\bar{y}}f^{\star})=\sum\limits_{t\in\bar{t}}\left(y_t-f^{\star}(y_{t-1})\right)E_{y_t}. $

因此

$ \|S^T_{\bar{y}}(y - S_{\bar{y}}f^{\star})\|^2_K=\sum\limits_{t\in\bar{t}}\sum\limits_{t'\in\bar{t}}\left(y_t-f^{\star}(y_{t-1}) \right)\left(y_{t'}-f^{\star}(y_{t'-1})\right)\langle E_{y_t}, E_{y_{t'}}\rangle_K. $

由于我们假设了$ \eta_t $的独立性以及$ E(y_t-f^{\star}_{y_{t-1}})=0 $, $ E\left\{(y_t-f^{\star}(y_{t-1}))^2\right\}=\sigma^2_{t} $, 两边取期望, 得到

$ E\left(\|S^T_{\bar{y}}(y-S_{\bar{y}}f^{\star})\|^2_K\right)=\sum\limits_{t\in\bar{t}}\sigma^2_t\langle E_{y_t}, E_{y_t}\rangle_K. $

根据$ {\cal H}_{K, \bar{y}} $为有界, 且最大范数小于等于$ \|J\|^2 $, $ \langle E_{y_t}, E_{y_t}\rangle_K\leq \|J\|^2 $. 因此误差期望的上界为

$ E\left(\|\tilde{f} - f^{\star}_{\bar{y}, \gamma}\|^2_{K}\right)\leq \|(S^T_{\bar{y}}S_{\bar{y}}+\gamma I)^{-1}\|^2\|J\|^2\sigma^2. $

根据Markov不等式, 得到第一个误差估计的概率.

对于第二个不等式, 根据定理2.3, 假设$ w\equiv 1 $. 由于$ \{\eta^2_t\}_{t\in\bar{t}} $且$ E(\eta^2_x)=\sigma^2_t $, 可知$ E(\eta^2_t)=\sigma^2_t $. 因此对$ \forall \epsilon > 0 $,

$ \mbox{P}_{y}\left\{\sum\limits_{t\in\bar{t}}\{\eta^2_t-\sigma^2_t\}> \epsilon\right\}\leq \mbox{exp}\left\{-\frac{\epsilon}{2M^2}\mbox{log}\left(1 + \frac{M^2\epsilon}{\sum\limits_{t\in\bar{t}}\sigma^2(\eta^2_t)}\right)\right\}, $

其中$ |\eta_t|\leq M $蕴含着$ |\eta^2_t-\sigma^2_t|\leq M^2 $, 因此

$ \sum\limits_{t\in\bar{t}}\sigma^2(\eta^2_t)\leq \sum\limits_{t\in\bar{t}}E(\eta^4_t)\leq M^2\sigma^2<\infty. $

根据

$ \|\tilde{f}-f^{\star}_{\bar{y}, \gamma}\|^2_{K}\leq \|L\|^2\|\{\eta_t\}\|^2_{l^2(t)}, $

根据(2.7)式取$ \epsilon \leftarrow \epsilon \sigma^2 $, 则第二个不等式成立.

定理2.3  设$ \{\eta_t\}_{t\in\bar{t}} $为$ {{\Bbb R}} $中的独立随机变量, 且方差为$ \{\sigma^2_t\}_t $. 设$ w_t\geq 0 $且有界$ \|w\|_{\infty}<\infty $. 设$ \sigma^2_w:=\sum\limits_{t\in\bar{t}}w_t\sigma^2_t<\infty $, 且对$ \forall \quad j $, $ |\eta_t-E(\eta_t)|\leq M $几乎处处成立. 对$ \forall \epsilon $, 有

(2.7) $ \begin{equation} \mbox{Prob}\left\{\sum\limits_{t\in\bar{t}}w_t[\eta_t-E(\eta_t)]>\epsilon\right\}\leq 2\mbox{exp}\left\{-\frac{\epsilon}{2\|w\|_{\infty}M}\mbox{log}\left(1 + \frac{M\epsilon}{\sigma^2_w}\right)\right\}. \end{equation} $

证明见附录.

针对多响应回归问题, 记$ y_t=(y_{t, 1}, \cdots, y_{t, p}), t\in\bar{t} $. 考虑以下正则化最小二乘问题

(2.8) $ \begin{equation} (\tilde{f}_j)^p_{j=1}:\mathop{\rm arg\min}\limits_{f_1, \cdots, f_p\in{\cal H}_{K, \bar{y}}} \left\{\sum\limits_{t\in\bar{t}}\mathop \sum \limits_{j = 1}^p (y_{t, j} - f_{j}(y_{t-1}))^2+\mathop \sum \limits_{j = 1}^p \gamma_j\|f_j\|^2_K\right\}, \end{equation} $

其等价于通过正则化最小二乘法分别构建$ f_j $, 即

$ \iff \mathop{\rm arg\min}\limits_{f_j\in{\cal H}_{K, \bar{y}}} \left\{\sum\limits_{t\in\bar{t}}(y_{t, j} - f_{j}(y_{t-1}))^2+\gamma_j\|f_j\|^2_K\right\}, \quad j=1, \cdots, p. $

能够将整体优化问题转化为分别优化$ \{f_j\}^p_{j=1} $, 其原因在于各响应之间的独立性. 算法1详细阐述了如何实现多响应回归问题.

算法1. 多变量情形

1) 将原始观测数据拓展为$ p $份独立观测数据$ (y_0, y_{1, j}) $, $ (y_1, y_{2, j}) $, $ \cdots $, $ (y_{T-1}, y_{T, j}) $, $ j=1, \cdots, p $.

2) 对于每一份观测数据, 运用正则化最小二乘法求解$ \tilde{f}_j(\cdot):=\sum\limits_{t\in\bar{t}}c_tK_t(\cdot) $.

3) 输出多状态估计函数式$ (\tilde{f}_j)^p_{j=1} $.

两个模拟实验被用来检验KSSM模型的有效性. 在所有实验中, 时间序列数据都来源于确定性方程, 外加混合一定的观测噪声.

例3.1  假设真实系统为

$ x_t=\sin(\frac{\pi}{12}t), \; \; \; y_t=x_t + \epsilon_t, $

其中噪声满足$ \epsilon_t\sim N(0, 0.1) $. 分别在时间点$ t=0, \cdots, 199 $进行$ 200 $次观测得到观测样本$ y_t $.

通过样本直观, 可以发现其具有明显的周期性特点. 且该系统不满足马尔可夫性, 因此尝试进行$ 2 $-步延迟转化

$ x^{\star}_t:=(x_t, x_{t-1}). $

给定适当的超参数, $ \gamma=0 $, 选取核函数$ \exp\left(-\frac{\|x_1 - x_2\|^2_2}{2* 0.2^2}\right) $. 随机抽取$ 20 $个样本点构造特征函数.

例3.2 (多响应回归)  假设系统为

$ x_t=\sin(x_{t-1} + \pi), \; \; \; y_t=x_t + \epsilon_t, $

其中$ \epsilon_t\sim N(0, 0.02) $, 且边界条件为$ x_0=1 $. 抽取指定时间节点$ t=0, \cdots, 51 $的$ 52 $个观测值, 构成时间序列样本$ y_t $.

这里, 首先考虑$ 2 $-步延时转换

$ x^{\star}_t:=(x_t, x_{t-1}). $

取$ \gamma=0.2 $且核函数为$ \exp\left(-\frac{\|x_1 - x_2\|^2_2}{2* 0.2^2}\right) $. 随机抽取$ 36 $个样本点作为特征函数的出发点.

从实验设计得到的样本案例, 可以发现核状态空间模型在周期性时间序列样本(如例3.1)和趋势性周期性时间序列样本(如例3.2)能够拥有良好的表现.

在这一部分, 对比KSSM模型和滑动平均模型(ARMA)以及Kalman Filter(Kalman)模型. 滑动平均模型起源于自回归模型, 可以追述到Akaike(1969), 用于识别白噪声条件下的结构化参数^{[10, 11]}. Kalman滤波是一种二次估计算法, 最先由Rudolf Kalman在1960提出^[12]. 现实实验数据由中国学位研究生发展中心提供. 该实验的目的是在给定相关影响因子的前提下, 预测机场的能见度水平. 具体的, 机场气压、着陆区域最高点气压以及海平面气压都会提供. 另外, 数据测量了每分钟的温度、湿度、风速等数据. 观测持续$ 2 $天时间. 我们将前面提到的两种估计算法和本文提出的状态空间模型都运用到该数据集, 比较了包括均方误差, 平方误差和总误差以及平均误差率在内的$ 4 $项指标进行对比实验. 结果见表 1.

我们对比了平均平方误差(MSE), 平方误差和以及平均绝对误差百分比(MAPE). 从表 1的对比结果, 可以看到本文提出的KSSM模型拥有更高的预测准确率. 进一步分析预测结果, 我们发现KSSM模型更加适应于规律化学习, 特别当数据具有明显的周期性时, 表现更加优秀. 这很容易理解, 是因为KSSM模型更加关注历史规则的学习. 至少在机场能见度预测数据中, 本文提出的核状态空间模型提供了更加准确的预测效果.

命题5.1  设$ X, Y $为线性空间, $ D $是$ X $的子集. 设算子$ T:D\rightarrow Y $为有界线性算子. $ D $为算子$ T $的定义域, 记为$ D(T) $, 则算子$ T $在$ D $上的范数定义为

$ \|T\|=\sup\left\{\frac{\|Tx\|}{\|x\|}:x\neq 0, x\in D(T)\right\}. $

如果$ T $为有界算子, 则有

$ \|Tx\|\leq\|T\|\|x\|, x\in D(T). $

这在很多泛函分析的书本中都可以找到. 下面是定理2.3的证明.

证  设$ E(\eta_t)=0 $, 记$ \eta_t $的方差为$ \sigma^2_t=E(\eta^2_t) $. 假设$ |\bar{t}|<\infty $, 则有

(5.1) $ \begin{equation} P:=\mbox{Prob}\left\{\sum\limits_{t\in\bar{t}}w_t\eta_t>\epsilon\right\}\leq \exp\left\{-\frac{\epsilon}{2\|w\|_{\infty}M}\mbox{log}\left(1 + \frac{M\epsilon}{\sigma^2_w}\right)\right\}. \end{equation} $

设$ c $为任意正常数, 根据独立性有

$ \begin{eqnarray*} P&=&\mbox{Prob}\left\{\exp\left\{\sum\limits_{t\in\bar{t}}cw_t\eta_t\right\}>e^{c\epsilon}\right\}\\ &\leq & e^{-c\epsilon}E\left(\exp\left\{\sum\limits_{t\in\bar{t}}cw_t\eta_t\right\}\right)\\ &=&e^{-c\epsilon}\Pi_{t\in\bar{t}}E\left(\exp\left\{cw_t\eta_t\right\}\right). \end{eqnarray*} $

因为$ |\eta_t|\leq M $几乎处处成立, 因而

$ E\left(\exp\left\{cw_t\eta_t\right\}\right)=1+\mathop \sum \limits_{k = 2}^{ + \infty } \frac{c^kw^k_tE(\eta^k_t)}{k!} \leq 1 + \mathop \sum \limits_{k = 2}^{ + \infty } \frac{c^kw^k_tM^{k-2}\sigma^2_t}{k!}. $

由于$ w_t\leq \|w\|_{\infty} $ and $ 1+t\leq e^t $, 可以得到

$ \begin{eqnarray*} E\left(\exp\left\{cw_t\eta_t\right\}\right)&\leq&\exp\left\{\mathop \sum \limits_{k = 2}^{ + \infty } \frac{c^k\|w\|^{k-1}_{\infty}M^{k-2}w_t\sigma^2_t}{k!}\right\}\\ &=&\exp\left\{\frac{e^{c\|w\|_{\infty}M}-1-c\|w\|_{\infty}M}{\|w\|_{\infty}M^2}w_t\sigma^2_t\right\}. \end{eqnarray*} $

因此

$ P\leq \exp\left\{ -c\epsilon + \frac{e^{c\|w\|_{\infty}M}-1-c\|w\|_{\infty}M}{\|w\|_{\infty}M^2}\sum\limits_{t\in\bar{t}}w_t\sigma^2_t \right\}. $

选择常量$ c $使得右式最小化

$ c=\frac{1}{\|w\|_{\infty}M}\log\left(1 + \frac{M\epsilon}{\sum\limits_{t\in\bar{t}}w_t\sigma^2_t}\right), $

其满足$ e^{c\|w\|_{\infty}M} - 1=\frac{M\epsilon}{\sigma^2_w} $. 从而

$ P\leq \exp\left\{ -\frac{\epsilon}{\|w\|_{\infty}M}\left\{ \left(1 + \frac{\sigma^2_w}{M\epsilon}\right)\log\left(1 + \frac{M\epsilon}{\sigma^2_w}\right) - 1 \right\} \right\}. $

设$ g(\lambda)=(1 + \lambda)\log(1+ \lambda) - \lambda $, $ \lambda \geq 0 $, 则有

$ P\leq \exp\left\{ -\frac{\sigma^2_w}{\|w\|_{\infty}M^2}g\left(\frac{M\epsilon}{\sigma^2_w}\right) \right\}. $

因为

$ g(\lambda)\geq\frac{\lambda}{2}\log(1+\lambda), \quad \forall \lambda\geq 0. $

为了得到该不等式, 考虑$ {{\Bbb R}} ^+ $上的函数

$ f(\lambda):=2\log(1+\lambda) - 2\lambda + \lambda\log(1+ \lambda), \quad \lambda\geq 0. $

容易得到$ f(0)=0, f'(0)=0 $, 且$ f^{{'}{'}}(\lambda)=\lambda(1+\lambda)^{-2} $, 只要$ \lambda\geq 0 $. 因此$ f(\lambda)\geq 0 $且

$ \log(1+\lambda) - \lambda \geq -\frac{1}{2}\lambda\log(1+\lambda), \quad \forall \lambda\geq 0. $

从而

$ g(\lambda)=\lambda\log(1+\lambda) + \log(1+\lambda) - \lambda\geq \frac{\lambda}{2}\log(1+\lambda), \quad \forall \lambda>0. $

至此, 我们证明了先前的声明. 不等式(5.1)成立, 只要令$ \lambda=\frac{M\epsilon}{\sigma^2_w} $. 定理2.3证毕.