状态空间模型中实际参数估计

（参数估计所有内容）

状态扩增法

将参数扩充为状态的一部分，例如有一个含有未知参数 $\mathbf{\theta }$ 的非线性模型为：
上述模型可以表述为：
两个参数的动态模型是一致的，重新定义状态为：$\tilde{\mathbf{x}}=({\mathbf{x}}_{\mathbf{k}},{\mathbf{\theta }}_{\mathbf{k}})$，则状态空间模型变为：

此时，模型中不含任何未知参数。状态扩增法在整个系统为线性时，可取得良好的效果。

线性状态空间模型的参数估计

假设一个含有未知参数 $\mathbf{\theta }$ 的线性高斯状态空间模型为：
其中,
线性高斯模型的能量函数：
能量函数的递归表达式：
式中，${\mathbf{v}}_{\mathbf{k}}(\mathbf{\theta }),{\mathbf{S}}_{\mathbf{k}}(\mathbf{\theta })$ 由只含 $\mathbf{\theta }$ 的卡尔曼滤波器给出。
卡尔曼滤波预测过程：
卡尔曼滤波更新过程：

因此，若已知 $\mathbf{\theta }$，先计算 $\mathbf{k}=0$ 时刻的值 ${\mathbf{\phi }}_{\mathbf{T}}=-\mathbf{l}\mathbf{o}\mathbf{g}\mathbf{p}(\mathbf{\theta })$，在依次执行上述算法直到 $\mathbf{k}=\mathbf{T}$，可得全局能量函数 ${\mathbf{\phi }}_{\mathbf{T}}(\mathbf{\theta })$。
已知能量函数后，可用基于Metropolis-Hastings的MCMC采样，产生后验分布的蒙特卡洛逼近等，得到参数等MAP或ML估计。

计算线性高斯模型中的 $\mathcal{Q}$：
如上诉所示的线性高斯模型，$\mathcal{Q}$ 的表达式为：
式中，所设计的参数由以 ${\mathbf{\theta }}^{(\mathbf{n})}$ 为参数的 $\mathbf{R}\mathbf{T}\mathbf{S}$ 平滑器 给出：

EM 算法对线性空间模型有效。

通过分别求：
上式中各个元素的偏导数 $\partial \mathcal{Q}(\mathbf{\theta },{\mathbf{\theta }}^{(\mathbf{n})})/\partial \mathbf{\theta }$，并分别令其为0，可得如下定理：

线性参数模型参数 $\mathcal{Q}$ 的极大化：
当参数选自模型参数：
时， ${\mathbf{\theta }}^{\ast }$ 的极大值为：
可由下式计算得到：

最终，得到关于初始均值 $\mathbf{\theta }={\mathbf{m}}_0$ 的极大值为：

线性状态空间模型的EM算法
令 $\mathbf{\theta }$ 包含模型参数 { A , H , Q , R , P 0 , M 0 } \bm\{{A,H,Q,R,P_0,M_0\}} {A,H,Q,R,P0​,M0​} 的部分子集，利用下述步骤，可获得参数的极大似然估计值。
首先，假设初始参数为值 ${\mathbf{\theta }}^{(0)}$。
对于 $\mathbf{n}=0,1,2,\cdot \cdot \cdot$ 执行如下步骤：
1). $\mathbf{E}$ 步骤：
利用 ${\mathbf{\theta }}^{(\mathbf{n})}$ 中当前的参数值进行 $\mathbf{R}\mathbf{T}\mathbf{S}$ 平滑器解算，并结合平滑结果计算下式：

也就是上上一个定理线性高斯模型中的Q中的式子。
2). $\mathbf{M}$ 步骤：
结合上一个定理：线性模型参数 $\mathcal{Q}$ 的极大化中的式子，计算新的参数，并将其储存在 ${\mathbf{\theta }}^{(\mathbf{n}+1)}$ 中。

$\mathcal{Q}(\theta ,{\theta }^{(n)})$ 的表达式也提供了一种简单的梯度计算方法，利用 $\mathbf{F}\mathbf{i}\mathbf{s}\mathbf{h}\mathbf{e}\mathbf{r}$ 恒定式计算能量函数的梯度。
$\mathbf{F}\mathbf{i}\mathbf{s}\mathbf{h}\mathbf{e}\mathbf{r}$ 恒定式:

利用高斯滤波与平滑的参数估计(非线性模型)

在非线性模型中，可以用相应的非线性卡尔曼滤波器及 $\mathbf{R}\mathbf{T}\mathbf{S}$ 平滑器进行参数估计。
考虑下述模型的参数估计问题：
式中，
此时，能量函数可由下述基于高斯滤波的算法得到。

基于高斯滤波的能量函数：
能量函数逼近的迭代表达式为：
式中，${\mathbf{v}}_{\mathbf{k}}(\mathbf{\theta }),{\mathbf{S}}_{\mathbf{k}}(\mathbf{\theta })$ 由只含 $\mathbf{\theta }$ 的高斯滤波器给出。
高斯滤波预测过程：

高斯滤波更新过程：

上述的能量函数可直接应用于MCMC采样。
当然，还有其他如计算ML或MAP估计中的无梯度最优化。

为计算 $\mathbf{E}\mathbf{M}$ 算法中所需要的期望值，可用高斯假设密度逼近法(如时矩匹配)对下式进行积分逼近：
式中，

积分逼近得到的 $\mathcal{Q}$ 表达式如下：
即：利用高斯平滑计算 $\mathcal{Q}$:
对本节中的非线性状态空间模型， $\mathcal{Q}$ 表达式如下：
式中的期望值由高斯平滑器得到。
实际中，高斯平滑器和高斯积分可由： $\mathbf{E}\mathbf{K}\mathbf{F}/\mathbf{E}\mathbf{R}\mathbf{T}\mathbf{S}\mathbf{S}$，或 $\mathbf{s}\mathbf{i}\mathbf{g}\mathbf{m}\mathbf{a}$ 点方法 进行逼近。$\mathbf{s}\mathbf{i}\mathbf{g}\mathbf{m}\mathbf{a}$ 点方法包括：Gauss-Hermite或球面容积积分，以及无迹变换等计算。
也可用上一节中线性模型极大化的方法：对于噪声协方差 $\mathcal{Q}$ 的极值表示为：

对于其他参数的 $\mathbf{M}$ 步骤实现，需依据实际的 $\mathbf{f},\mathbf{h}$ 而定。若是线性的，可用上述类似方法。
若为非线性， $\mathbf{F}\mathbf{i}\mathbf{s}\mathbf{h}\mathbf{e}\mathbf{r}$ 恒等式给出一种计算能量函数的方法。

基于粒子滤波与平滑的参数估计

粒子滤波也可用于逼近参数估计中的边缘似然函数及能量函数的计算，在粒子滤波中，考虑的一般模型为：

式中，$\mathbf{\theta }\in {\mathbb{R}}^{\mathbf{d}}$ 为未知参数。

算法：基于能量函数逼近的 $\mathbf{S}\mathbf{I}\mathbf{R}$:
参数边缘似然的逼近可以在执行**序列重要性采样(SIR)**算法(粒子滤波)中进行计算，
具体步骤为：
先从先验信息中抽取 $\mathbf{N}$ 个采样点 ${\mathbf{x}}_0^{(\mathbf{i})}$，即有：
并令 ${\mathbf{w}}_0^{(\mathbf{i})}=1/\mathbf{N},andi=1,\cdot \cdot \cdot ,N$。
对于每个 $\mathbf{k}=1,\cdot \cdot \cdot ,\mathbf{T}$，执行如下操作：
1). 从重要性分布中抽取采样 ${\mathbf{x}}_{\mathbf{k}}^{(\mathbf{i})}$，即有：
2). 计算权值：
并计算 $\mathbf{p}({\mathbf{y}}_{\mathbf{k}}|{\mathbf{y}}_{1:\mathbf{k}-1},\mathbf{\theta })$ 的估计值：
3). 计算归一化权值：
4). 若有效粒子数过少，则执行重采样。

参数边缘似然的逼近值为：
相应的能量函数逼近值为：

能量函数逼近可在执行 $\mathbf{S}\mathbf{I}\mathbf{R}$ 算法时进行递归解算。

$\mathbf{E}\mathbf{M}$ 算法可由粒子平滑实现。$\mathcal{Q}$ 逼近的实际表达式根据具体选用的粒子平滑器而定。若选用后向平滑器，可得以下算法：

算法：利用后向平滑器计算 $\mathcal{Q}$：
设已知模拟轨迹：
并利用后向平滑器，参数为 ${\mathbf{\theta }}^{(\mathbf{n})}$，并利用下式
可得上述 $\mathcal{Q}$ 式为：

若选用重新加权粒子平滑器，可得如下算法：

算法： 利用重新(调整)加权粒子平滑器计算 $\mathcal{Q}$;
设粒子集为 { x k ( i ) , k = 0 , ⋅ ⋅ ⋅ T , i = 1 , ⋅ ⋅ ⋅ N } \bm{\{x_k^{(i)},k=0,\cdot\cdot\cdot T, i=1,\cdot\cdot\cdot N\}} {xk(i)​,k=0,⋅⋅⋅T,i=1,⋅⋅⋅N}，并利用重新加权粒子平滑器计算权值 { w k ∣ T ( i ) , k = 0 , ⋅ ⋅ ⋅ T , i = 1 , ⋅ ⋅ ⋅ N } \bm{\{w_{k|T}^{(i)},k=0,\cdot\cdot\cdot T, i=1,\cdot\cdot\cdot N\}} {wk∣T(i)​,k=0,⋅⋅⋅T,i=1,⋅⋅⋅N}，并利用下式
可得上述 $\mathcal{Q}$ 式为：
可利用 $\mathbf{F}\mathbf{i}\mathbf{s}\mathbf{h}\mathbf{e}\mathbf{r}$ 特性计算能量函数梯度的逼近值。

参数的 Rao-Blackwell 化

利用Rao-Blackwell化方法可以边缘化状态空间模型中的静态参数。有下面的普通摆模型介绍：
设有：
式中，
量测噪声 $\mathbf{R}$ 的方差未知，只给出 $\mathbf{R}$ 的逆 $\mathbf{c}\mathbf{h}\mathbf{i}$ 平方先验分布。
获取 $\mathbf{R}$ 分布参数步骤为：
1). 设此前产生的粒子群为 { w k ( i ) , w 0 : k ( i ) , i = 1 , ⋅ ⋅ ⋅ , N } \bm{\{w_k^{(i)},w_{0:k}^{(i)},i=1,\cdot\cdot\cdot,N\}} {wk(i)​,w0:k(i)​,i=1,⋅⋅⋅,N} 该粒子群逼近的状态完全分布为：
结合给定的量测值以及采样状态信息时，$\mathbf{R}$ 的条件分布为：
由此，关于状态及参数的完全分布为：
2). 对重要性分布进行采样：
3). 结合 ${\mathbf{x}}_{\mathbf{k}}^{(\mathbf{i})}$ 和此前的量测量，可得当前量测量的似然为：
式中，$\mathbf{s}\mathbf{t}\mathbf{u}\mathbf{d}\mathbf{e}\mathbf{n}\mathbf{t}$ 的 $\mathbf{t}$ 分布参数为：
由此，可以计算出 $\mathbf{S}\mathbf{I}\mathbf{R}$ 算法中，下一步的重要性权值，进而可计算有效粒子数，并判断是否需要进行重采样：
4). 结合量测量 ${\mathbf{y}}_{1:\mathbf{k}}$ 及状态量 ${\mathbf{x}}_{0:\mathbf{k}}^{(\mathbf{i})}$，可进一步计算出 $\mathbf{R}$ 的条件分布：
5). 令 $\mathbf{k}\leftarrow \mathbf{k}+1$，跳至步骤 1)。

上述的方法就是 Rao-Blackwell粒子滤波器，其中静态参数 $\mathbf{R}$ 已被边缘化，并通过对 ${\mathbf{v}}_{\mathbf{k}}^{(\mathbf{i})},{\mathbf{R}}_{\mathbf{K}}^{(\mathbf{i})}$ 充分统计，结合过去粒子 ${\mathbf{x}}_{0:\mathbf{k}}^{(\mathbf{i})}$ 和量测量得到。