状态空间MCMC论文笔记

论文一：

模拟案例

$$\begin{gathered} y_t=n_t\ast c_t \\ {n}_t=a+n_{{}_{t-1}}+w_{{}_t} \\ {c}_t=a_{{}_1}{c}_{{}_{t-1}}+a_2{c}_{{}_{t-2}}+v_t \end{gathered}$$

理论推导

（一）整体思想

令${\vec{y}}_{{}_T}=[y_1,y_2,...,y_n{]}^{\prime }$，${\vec{n}}_{{}_T}=[n_1,n_2,...,n_n{]}^{\prime }$,${\vec{c}}_{{}_T}=[c_1,c_2,...,c_n{]}^{\prime }$,$\vec{a}=[a_1,a_2]$
估计上述状态空间模型的目的是求出参数和状态在给定观测值${\vec{y}}_{{}_T}$下的条件分布，即要得的${\vec{n}}_{{}_T},{\vec{c}}_{{}_T}$以及$\vec{a}$的后验分布。
第一步，求状态在观测序列下的条件分布：
$$p({\vec{n}}_{{}_T}|{\vec{y}}_{{}_T})$$
直接求上述的后验分布比较麻烦，可以对上述的分布进行分解：
$$\begin{gathered} \pi ({\vec{n}}_{{}_T}|{\vec{y}}_{{}_T})=\pi (n_{{}_1},n_{{}_2},...,n_{{}_T}|{\vec{y}}_{{}_T}) \\ =\pi (n_{{}_T}|{\vec{y}}_{{}_T})\prod _{i=1}^T\pi (n_t|n_{{}_{t+1}},{\vec{y}}_{{}_t})\text{\hspace{1em}(8)} \end{gathered}$$
当t=T时，状态$n_{{}_t}$的抽样
利用(8)式，由右式第一项可知$\pi (n_{{}_T}|{\vec{y}}_{{}_T})$，即$n_{{}_T}$的后验分布仅仅依赖于观测序列${\vec{y}}_{{}_T}$，因此可以首先利用$\pi (n_{{}_T}|{\vec{y}}_{{}_T})$对$n_{{}_T}$进行抽样，得到$n_{{}_T}$的一个抽样值。
当t<T时，状态$n_{{}_t}$的抽样
根据(8)式，由右式第二项可知，$n_{{}_t}$的分布依赖于$n_{{}_{t+1}}$和${\vec{y}}_{{}_T}$，我们可以根据第一种情况，得到$n_{{}_T}$，然后由$n_{{}_T}$和${\vec{y}}_{{}_T}$得到$n_{{}_{T-1}}$，依次这样倒着抽样，就可以得到所有的 $\pi (n_t|n_{{}_{t+1}},{\vec{y}}_{{}_t})$
再根据t=T，t<T得出的抽样值，就可以计算出方程右侧的抽样值了。也就说，可以利用式（8）的右侧分解的抽样，来完成对左侧的抽样。

（二）条件分布的具体形式

（8）式右侧抽样涉及到两部分（t=T时，$\pi (n_{{}_T}|{\vec{y}}_{{}_T})$的分布形式，t<T时的$\pi (n_t|n_{{}_{t+1}},{\vec{y}}_{{}_t})$分布形式），由于假设状态空间模型的状态和随机干扰项都是高斯的，因为条件分布也是高斯分布，所以得到：
$$n_{{}_T}|{\vec{y}}_{{}_T}\sim N(E(n_T|{\vec{y}}_{{}_T}),cov(n_T|{\vec{y}}_{{}_T}))\text{\hspace{1em}(9)}$$
$$\pi (n_t|n_{{}_{t+1}},{\vec{y}}_{{}_t})\sim N(E(n_t|{\vec{y}}_{{}_t},n_{{}_{t+1}}),cov(n_t|{\vec{y}}_{{}_t},n_{{}_{t+1}})\text{\hspace{1em}(10)}$$
再来看（9）式代表的意义，该式表示要根据观测序列，估计状态$n_T$的后验分布，我们现在知道两种信息：1.观测值序列；2.状态转移方程；由这两个信息就可以得到$n_T$的分布估计，并且如果估计方法是卡尔曼滤波，则这个后验分布还是最佳估计（方差最小）。因此，可以用卡尔曼滤波得到（9）式。
假设（9）式的分布解决了，就可以抽一个样本值，将这个样本值带入到（10）式，可以得到当$t<T$时，分布参数的递推公式：
$$E(n_t|{\vec{y}}_{{}_t},n_{{}_{t+1}})=E(n_t|{\vec{y}}_{{}_t})+cov(n_t|{\vec{y}}_{{}_t})F^{\prime }[Fcov(n_t|{\vec{y}}_{{}_t})F^{\prime }+Q{]}^{-1}[n_{{}_{t+1}}-\mu -FE(n_t|{\vec{y}}_{{}_t})]\text{\hspace{1em}(11)}$$
问题：$\mu$是什么。
$$Cov(n_t|{\vec{y}}_{{}_t},n_{{}_{t+1}})=Cov(n_t|{\vec{y}}_{{}_t})+cov(n_t|{\vec{y}}_{{}_t})F^{\prime }[Fcov(n_t|{\vec{y}}_{{}_t})F^{\prime }+Q{]}^{-1}FCov(n_t|{\vec{y}}_{{}_t})\text{\hspace{1em}(12)}$$
按照(11),(12)得到$n_t(t=1,2,...,T)$后，需要更新参数的值（包括系数截距项$a$，另外一个状态变量$c_t$，以及两个随机干扰项的分布参数）。
第二步，更新参数$a$的值，根据式（6）得到a的估计值；再根据式（5）计算另一个状态$c_t$的估计值。
第三步，状态向量已知，参数$a$已知，现在对两个随机干扰项进行抽样。
（1）在状态$n_t$已知，参数$a$已知的情况下，对$w_t$的标准差进行采样
（2）在状态$c_t$已知，参数$a$已知的情况下，对$v_t$的标准差进行采样
（3）在状态和其他参数已知的情况下，对系数向量$\vec{a}=[a_1,a_2]$进行采样
完成采样后，就完成了参数的一轮更新，再把更新值作为初始值，重复第一步到第三步。
问题的关键，找出参数的满条件分布形式。

论文2：

状态方程：
$$\{\begin{array}{l}{y}_t=F_{{}_t}^{\prime }\overrightarrow{{\theta }_t}+v_t,v_t\sim N(0,V)\\ \overrightarrow{{\theta }_t}=G_t{\theta }_{{}_{t-1}}+w_t,w_t\sim N(0,W)\end{array}$$
其中，$V$是标量，$W=\left[\begin{array}{ccc}{w}_{11}& 0& 0\\ 0& w_{22}& 0\\ 0& 0& w_{33}\end{array}\right]$，其中$V,W$称为精度参数。

在精度参数$V,W$已知的情况下，对状态参数$\theta$采用FFBS抽样，根据Carter等提出的FFBS算法理论，可以对状态参数向量空间进行前项滤波，后向Gibbs抽样，从而获得时变参数的向量值。首先给定初始信息：
$${\theta }_1|D_1\sim N(\overrightarrow{m_1},C_1)$$
在时刻$t=1,2,...,T$上，求：
$${\theta }_t|D_t\sim N(\overrightarrow{m_t},C_t)\text{\hspace{1em}(12)}$$
如论文一中我们的分析，此时知道状态的初始值，知道状态方程，知道观测方程，可以用卡尔曼滤波法求出式（12）的具体表达形式，也就是说可以求出均值和协方差矩阵，由如下式子给出：
$$\{\begin{array}{l}{m}_t=E(\overrightarrow{{\theta }_t|D_t})=G_t\cdot E(\overrightarrow{{\theta }_{{}_{t-1}}|D_{t-1}})+R_t\cdot F_t\cdot Q_t^{-1}\cdot [y_t-f_t]\\ \overrightarrow{f_t}=F_t^{\prime }\cdot G_t\cdot m_{{}_{t-1}}\end{array}$$

应用：负荷预测天气影响因子分析

本文将采用FFBS算法，研究天气对负荷的影响效应，并建立变参数模型，研究天气的影响随时间的变化关系。本文的逻辑结构如下：

• 首先，建立状态方程，并用kalman滤波来估计变参数模型；
• 其次，指出kalman滤波在估计变参数时反应“迟钝”的问题，并给出初步解决思路；
• 再次，解决思路依赖于突变点的位置，因此，建立贝叶斯模型，推断突变点的位置；
• 最后，利用验证的模型框架，验证天气和负荷之间的动态关系。

一、模型假设

假设当天的负荷可以分解为天气影响的部分、非天气因子影响的部分以及随机干扰的部分，即天气增长率可以写成如下方程：
状态方程：
$$\{\begin{array}{l}{y}_t={\alpha }_{{}_t}\cdot k_t+{\beta }_t+v_t,v_t\sim N(0,V)\\ {\alpha }_{{}_t}={\alpha }_{{}_{t-1}}+w_{{}_{1t}},w_{{}_{1t}}\sim N(0,w_1)\\ {\beta }_{{}_t}={\beta }_{{}_{t-1}}+w_{{}_{2t}},w_{{}_{2t}}\sim N(0,w_2)\end{array}\text{\hspace{1em}(1)}$$
可以将上述方程写成矩阵的形式：
令$\overrightarrow{{\theta }_t}=[{\alpha }_t,{\beta }_t{]}^{\prime }$，$X_t=[k_t,1{]}^{\prime }$，$G_t=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]$，$W^{\prime }=[w_{{}_{1t}},w_{{}_{2t}}]$，则$W\sim N([0,0{]}^{\prime },\left[\begin{array}{cc}{w}_1& 0\\ 0& w_2\end{array}\right])$
则方程（1）可以写作：
$$\{\begin{array}{l}{y}_t=X_t^{\prime }\overrightarrow{{\theta }_t}+v_t,v_t\sim N(0,V)\\ \overrightarrow{{\theta }_t}=G_t\cdot {\vec{\theta }}_{{}_{t-1}}+W\end{array}\text{\hspace{1em}(2)}$$

二、实现算法

算法伪代码

初始化V，W
for i in range(迭代次数)：
	1.求theta_t的卡尔曼滤波;也就是theta_t的分布。
	2.时间倒序对theta进行采样：
	for j in range(T):
		求出theta的后验分布，并从中抽样；
	3.根据theata的采样，对V，W进行采样，更新V，W的值。

python代码实现

代码整体结构：

kalman滤波

本节将实现线性卡尔曼滤波器。卡尔曼滤波的本质是定义个可观测的变量，这个变量满足1.可测量；2.可通过隐含的不可观测的变量转换而来。也就是说这个可观测变量即可测量，又可预测，那么怎么把测量和预测的信息结合起来，就是kalman滤波要解决的问题，kalman滤波假设测量值服从一个正态分布，预测值也是一个正态分布，二者相乘，就得到了方差最小的可观测变量的估计值。具体代码如下：

import numpy as np
import matplotlib.pyplot as plt
def init_theta():
    # 状态方程参数1.初始状态的协方差；2.初始状态的干扰协方差
    # 初始状态分布的协方差
    P0 = np.mat([[2.5,0], [0, 2.5]])
    # 初始状态的随机干扰的协方差
    Q0 = np.eye(2)
    # 状态的初始值
    X0 = np.mat([[0.0, 0.0]])
    # 量测方程的随机干扰
    R0 = np.mat([[3.0]])
    # 量测方程的初始均值
    # R0 = np.mat(norm.rvs(loc=0, scale=1, size=2, random_state=1)).T
    # R0 = np.mat([[0.5,0], [0, 0.5]])
    # 状态转移矩阵:单位矩阵F
    F0 = np.eye(2)
    # 状态到观测的映射矩阵H
    return P0, Q0, R0, F0, X0

def kalman_filter(observations, H):
    """
    实现kalman滤波
    :param observations:观测值，也就是Z
    :param H: 状态到观测变量映射的系数即mu = HX
    :return:
    """
    # 初始化变量
    P0, Q0, R0, F0, X0 = init_theta()
    # 更新方程
    T = len(observations)
    # 均值的更新方程
    all_mean = [X0]
    all_cov = [P0]
    H = np.mat(H).T
    observations = np.mat(observations).T
    R_k = R0  # 观测误差协方差不变
    F_k = F0  # 状态转移矩阵不变
    P_K_1 = P0
    X_k_1 = X0.T
    for t in range(1, T):
        # 因为此时收敛了，所以协方差矩阵近乎为0，此时要出现系数突变
        H_k = H[:, t].T  # 解释变量,shape = [2,1]
        Z_k = observations[t,:]  # 观测值,shape=[1,1]
        # P_k的更新方程X_k = F_k*X_(k-1) ==> cov(X_k) = F_k*cov(X_(k-1))*F_k'
        P_k = F_k*P_K_1 * F_k.T
        #卡尔曼增益K
        K = (P_k*H_k.T*np.linalg.pinv(H_k*P_k*H_k.T+R_k)).T
        X_hat_k = F_k * X_k_1
        # 更新均值
        adjust_x_k = X_hat_k + K.T*(Z_k-H_k*X_hat_k)
        # 更新协方差
        adjust_p_k = P_k - K.T*H_k*P_k
        # 更新参数
        P_K_1 = adjust_p_k
        X_k_1 = adjust_x_k
        all_mean.append(adjust_x_k)
        all_cov.append(P_K_1)
    return all_mean, all_cov

kalman滤波的案例

下面将构造一个固定系数的一元线性回归模型，模拟生成如下数据：
$$y=100x+5+ϵ,x\in [0,1],ϵ\sim N(0,1)\text{\hspace{1em}(3)}$$
然后利用kalman滤波估计方程的斜率和截距。

 number = 2000
 x = np.random.rand(number)
 y = 100 * x + 5 + np.random.randn(number)
 plt.figure(figsize=(12,8))
 plt.plot(y)
 plt.legend('y')
 plt.show()

下面将利用第一节中定义的状态空间模型，估计截距和斜率。

 X = np.vstack([x, np.ones(number)]).T
 means,covs = kalman_filter(observations=y, H=X)
 means = np.hstack(means[1:])
 ax1 = plt.subplot(2, 1, 1, frameon=False) 
 plt.plot(means[0, :][0].tolist()[0], 'b--')
 plt.ylabel('k')
 ax2 = plt.subplot(2, 1, 2)
 plt.plot(means[1, :][0].tolist()[0], 'r--')
 plt.ylabel('b')
 plt.show()

如下图，给定初始值后，kalman滤波很快就收敛到了参数的真实值。

kalman滤波在变系数估计的局限性

天气因子在不同的季节影响效应是不同的，商业体在夏季空调是开启的，经过一段时间的过渡，进入秋冬季，此时空调可能大部分是关闭的，或者完全关闭，因此，天气对负荷的影响效应可能会有一个急剧的变化。本节将模拟一个变系数线性回归模型，在默认参数情况下，验证kalman滤波能否快速识别新的用电模式。模拟生成如下数据：

y t = { 12 x t + 5 + ϵ , t ∈ { 0 , 1 , . . . , 499 } 20 x t + 5 + ϵ , t ∈ { 500 , 501 , . . . , 1999 } (4) y_t=\begin{cases} 12x_t+5+\epsilon,t\in \{0,1,...,499\} \\ 20x_t+5+\epsilon,t \in \{500,501,...,1999\} \end{cases} \tag{4} yt​={12xt​+5+ϵ,t∈{0,1,...,499}20xt​+5+ϵ,t∈{500,501,...,1999}​(4)

number = 2000
x = np.random.randn(number)
# x.sort()
y = []
for i in range(number):
	if i<500:
	    y.append(12*x[i]+5 + np.random.randn())
	else:
	    y.append(20 * x[i] + 5 + np.random.randn())
y = np.array(y)
plt.figure(figsize=(12,8))
plt.legend("y")
plt.plot(y)
plt.show()

利用kalman滤波估计变系数的结果如下：

kalman滤波估计变系数存在的问题；我们人为的在t=499处设置了一个结构突变，来模拟天气效应的变化，发现kalman滤波更过很长一段时间仍然没有收敛到真实突变的水平，kalman滤波是一个缓慢上升的过程，原因在于kalman本质上是加权平均，在第一种数据生成的过程中，kalman滤波收敛到了真实的系数水平，此时的预测模型的协方差很小，在加权平均时，权重很大，导致后面新采集的数据在加权中，权重很小，因此，要经过很长一段时间的运行，才慢慢朝着真实的第二个系数逼近，为了验证这个推测，在运行kalman滤波时，在t=500处，重置一下预测模型的协方差，让协方差为单位矩阵，效果如下：

def kalman_filter(observations, H):
    """
    实现kalman滤波
    :param observations:观测值，也就是Z
    :param H: 状态到观测变量映射的系数即mu = HX
    :return:
    """
    # 初始化变量
    P0, Q0, R0, F0, X0 = init_theta()
    # 更新方程
    T = len(observations)
    # 均值的更新方程
    all_mean = [X0]
    all_cov = [P0]
    H = np.mat(H).T
    observations = np.mat(observations).T

    R_k = R0  # 观测误差协方差不变
    F_k = F0  # 状态转移矩阵不变

    P_K_1 = P0
    X_k_1 = X0.T

    for t in range(1, T):
    
        # 因为此时收敛了，所以协方差矩阵近乎为0，此时要出现系数突变
        # 由于协方差近乎为0，导致加权平均时，导致预测值的权重高，而测量
        # 值的权重低，因此，会出现很长一段时间都找不到正确的突变系数
        if t == 500:  # 在这里重置预测模型的协方差矩阵
            P_K_1 = np.mat([[1.0, 0], [0, 1.0]])
        其余代码同上...

从上图来看，期间只需简单的在迭代到500处时，重置一下协方差矩阵，kalman滤波就能快速的找到真实的第二阶段的斜率。针对何时需要调整协方差矩阵，是一个可以继续研究的点，本文的重点是将协方差矩阵随机化，因此，这里就不在继续讨论这个问题了。

贝叶斯推断定位突变点

        本节将描述如何利用贝叶斯方法，来推断数据的突变点个数和突变点位置。这里的突变点可以理解为分段函数；仍然以下面的分段函数为例，来验证我们检测突变点的方法；
y t = { 12 x t + 5 + ϵ t , t ∈ { 0 , 1 , . . . , 499 } , 20 x t + 5 + ϵ t , t ∈ { 500 , 501 , . . . , 1999 } (5) y_t=\begin{cases} 12x_t+5+\epsilon_t,t\in \{0,1,...,499\}, \\ 20x_t+5+\epsilon_t,t \in \{500,501,...,1999\} \end{cases} \tag{5} yt​={12xt​+5+ϵt​,t∈{0,1,...,499},20xt​+5+ϵt​,t∈{500,501,...,1999}​(5)
   突变中，最重要的一种情形就是数据的均值从一个平台转换到另一个平台，从时间序列的角度来看，也就是说以往历史的样本，和最近的样本均值发生了变化，这种也是现实中常见的情形，为了模拟这种场景，在实验时，要求$x_t$是平稳的，且均值非0，否则的话，突变应该是属于波动率突变的情形了，那种情形应该对波动率进行结构突变的检验，总之，这里将演示如何构造模型，检验均值突变。
当含有一个结构突变点时，根据实y取值的分布，常常可以构造如下两种类型的贝叶斯模型：

• 当观测值y取值可正可负，则可以假设观测值的分布是正态分布：
  $$\{\begin{array}{l}{y}_t\sim N(\lambda ,{\sigma }_2)\\ \lambda =\{\begin{array}{l}{\lambda }_1\sim exp(1/\bar{y}),t<\tau \\ {\lambda }_2\sim exp(1/\bar{y}),t\ge \tau \end{array}\\ \tau \sim DiscreteUniform(1,T)\\ \sigma \sim U[0,2{\sigma }_y]\end{array}\text{\hspace{1em}(3)}$$
  其中，${\sigma }_y$表示观测值y的标准差；$\bar{y}$表示观测值y的均值。T表示观测值的数量。

• 当观测值是非负的情况，我们既可以假设其实正态分布，也可以假设其实泊松分布，或者指数分布，这里给出泊松分布的例子：
  $$\{\begin{array}{l}{y}_t\sim Possion(\lambda )\\ \lambda =\{\begin{array}{l}{\lambda }_1\sim exp(1/\bar{y}),t<\tau \\ {\lambda }_2\sim exp(1/\bar{y}),t\ge \tau \end{array}\\ \tau \sim U[1,T]\end{array}\text{\hspace{1em}(3)}$$

上述两种思路的python实现和结果如下：

状态空间后采样

待续

精度参数gibbs抽样

待续