贝叶斯滤波框架概述

贝叶斯滤波是采用贝叶斯原理来处理滤波问题的滤波方法。为了透彻理解贝叶斯滤波，我们首先从统计学中的估计问题谈起。对于某个总体，若其概率分布已知则称为总体已知，因此总体是定义在概率空间上的概率测度。我们可以对总体做一些独立实验，从而得到一系列的独立同分布（i.i.d）的样本 $X_{i}$ ，它是一个随机变量，其取值为 $x_{i} \in R^{d_{x}}$ （可以拓展到复数域），是一个变量。估计的关键在于通过有限的样本获得我们期望的估计量。估计量是一种特殊的统计量，统计量的定义为样本的一个已知函数 $T (X)$ ，当样本取值确定时统计量的取值也随之确定，因此统计量是一个随机变量，服从一定的分布，虽然这个分布可能并不一定存在解析表达式。常用的统计量有：样本平均值、样本方差和次序统计量。包含一定目的的统计量为估计量，一般假设带估计参数为 $θ$ ，那么估计量为函数 $\hat{θ} (X)$ ，同样估计量也是一个随机变量。

频率学派参数估计

对于待估参数 $θ$ ，经典参数估计理论认为它是一个固定的未知量，是一个变量，不服从任何分布假设。用柏拉图的洞穴比喻来说明，待估参数就是真实的物体，而我们看到的（做的独立实验）是投影在洞穴壁上的影子，我们需要观察影子来估计真实的待估参数。

与估计量相关的几个概念定义如下：

误差： $e (x) = \hat{θ} (x) - θ$ ；
均方误差（MSE）： $M S E (\hat{θ}) = E [(\hat{θ} (X) - θ)^{2}]$ ，MSE描述了统计量和真实值之间的误差；
偏差： $B (\hat{θ}) = E (\hat{θ}) - θ$ ， $B (\hat{θ}) = 0$ 的估计量为无偏估计量
样本偏差： $d (x) = \hat{θ} (x) - E [\hat{θ}]$ ，样本偏差是给定样本的情况下的一个值；
方差： $var (\hat{θ}) = E [(\hat{θ} - E [\hat{θ}])^{2}]$ ，方差描述了估计量的性能；

其中估计量的均方误差、偏差和方差均为描述估计量性能的指标。他们之间的的关系为 $M S E (\hat{θ}) = var (\hat{θ}) + (B (\hat{θ}))^{2}$ 在传统的估计框架下，由于 $θ$ 是固定的未知量，因此MSE往往依赖于 $θ$ ，从而导致无法调整估计量来使得最小进而获得最小MSE估计量。因此往往限定估计量是无偏的，再进一步寻找方差最小的估计量，称为最小方差无偏估计量（MVU），由（1）可知方差最小的无偏估计量在无偏估计量中也是MSE最小的。

尽管我们期望能够得到MVU估计量，但是它往往很难求得。CRLB给出了无偏估计量的下限，达到CRLB的估计量即为MVU估计量。可以藉由充分统计量来确定MVU，具体方法是通过Neyman-Fisher因子分解定理。在无法确定MVU估计量的情况下还存在其他的方法如：最大似然法、矩估计法和BLUE等。

贝叶斯估计

而在贝叶斯统计中，待估参数是一种「信念」。而我们提前拥有关于这个信念的一部分知识，称为先验。而对于事件的观测则会影响我们对该事物的知识（似然），影响后的知识为后验。先验知识和后验结果都以概率的形式体现出来：先验分布 $p (θ)$ ，似然 $p (x | θ)$ 和后验分布 $p (θ | x)$ 。而这三者通过贝叶斯公式关联起来：

\begin{matrix} (2) & p (θ | x) = \frac{p (x | θ) p (θ)}{\int p (x | θ) p (θ) d θ} \end{matrix}

在贝叶斯框架下对估计量性能进行衡量，最先想到的仍然是使用均方误MSE，与传统估计不同的是待估参数也是随机变量，因此有 $B M S E (\hat{θ}) = E [(θ - \hat{θ} (X))^{2}] = \iint (θ - \hat{θ})^{2} p (x, θ) d x d θ$ 应当注意定义中的 $θ$ 和 $\hat{θ}$ 的位置与经典框架下的区别。此时的贝叶斯MSE（BMSE）已经与待估参数无关，因此可以调整 $\hat{θ}$ 来使得BMSE最小从而获得最小均方误差（MMSE）估计量，对（3）关于 $\hat{θ}$ 求导并令其等于零可以得到MMSE估计量为后验期望： ${\hat{θ}}_{M M S E} = E [θ | x] = \int θ p (θ | x) d θ$ 更一般地，令误差 $e = θ - \hat{θ}$ ，定义关于误差的代价函数 $C (e)$ ，代价函数是一个确定性的函数，定义平均代价为贝叶斯风险 $R = E [C (e)]$ 。令不同的代价函数（准则）最小可以得到不同的估计量，总结如下表

代价函数 $C (e)$	估计量 $\hat{θ}$
$e^{2}$ 二次型误差	$E [θ \| x]$ ，后验PDF期望
$\| e \|$ 绝对值误差	后验PDF中值
成功-失败误差	后验PDF众数（最大值）

贝叶斯估计最令人难以接受的地方便是计算后验分布。尤其是在确定MMSE估计量时，首先需要计算（2）中分母上的积分（在实际问题中可能是多维积分）。然后还需要计算复杂的期望中的积分。不过Monte-Carlo方法为积分的数值求解提供了极大的便利，我们将在本文后续介绍Monte-Carlo方法。

小结

经典的参数估计理论把参数作为未知的确定量即变量，这往往导致估计量的MSE自身与待估参数有关，从而无法求出MMSE估计量，取而代之的是我们希望可以求得MVU估计量，它是无偏估计量中最优的。但是MSE准则下的最优估计量并不一定是无偏的。贝叶斯估计量把待估参数看作是随机变量，在求估计量的MSE时待估参数也同时被求期望，因此得到的MSE与参数无关，可以求得MMSE估计量，经过简单的推导可知MMSE估计量就是后验PDF的数学期望。

最优线性滤波器——LMMSE滤波器

在这一部分我们进入贝叶斯估计在信号处理领域的应用，首先考虑最一般性的滤波问题：通过输入数据 $X$ 来推测期望响应 $Y$ 。一般 $X$ 代表了当前时刻之前的 $M$ 个输入信号 $[x_{n - M + 1}, . . ., x_{n}]$ ， $Y$ 为被估计的参数 $y_{n}$ 。我们假定 $y$ 和 $x$ 均是零均值的。应用贝叶斯公式（2）可以求得后验分布为

p (y | x) = \frac{p (x | y) p (y)}{p (x)}

上述公式中的任意一项都很难求得，往往只能使用统计手段估计出来，这需要大量的数据量。一个简单的解决方案是假定输入信号之间不存在相关性，则似然函数可以写作 $p (x | y) = \prod_{i} p (x_{i} | y)$ ，这种方法在在朴素贝叶斯分类器中使用。

LMMSE滤波器推导

回想DSP中的数字滤波器，两种最经典的滤波器都是线性的：IIR和FIR滤波器。现在我们约束对 $y$ 的估计也是线性的： $\hat{y} = c^{H} x$ ，再从一定准则来推导线性系数 $c$ 。此时的估计器的MSE为系数 $c$ 的函数：

\begin{aligned} (1) & E [| y - \hat{y} ‖^{2}] & = E [(y - \hat{y})^{H} (y - \hat{y})] \\ (2) & = E [‖ y ‖^{2}] - c^{H} E [x y^{*}] - E [y x^{H}] c + c^{H} E [x x^{H}] c \end{aligned}

令上式中第一项为 $P_{y}$ 代表了期望响应的功率，第二项中的期望为 $d$ 代表了输入数据和期望响应中间的互相关，最后一项中的期望为输入信号的协方差矩阵。那么上式为系数的确定函数称为误差性能表面

$P (c) = P_{y} - c^{H} d - d^{H} c + c^{H} R c$ 协方差矩阵一定是半正定矩阵，实际数据的协方差矩阵通常都是正定矩阵。约束 $R$ 是正定矩阵，那么（4）可以转换为

P (c) = P_{y} - d^{H} R^{- 1} d + (R c - d)^{H} R^{- 1} (R c - d)

上式中只有最后一项与 $c$ 有关，因此当系数满足正则方程 ${\begin{cases} R c - d = 0 \\ R ≻ 0 \end{cases}$ 此时的线性滤波器的MSE最小，这种滤波器称为LMMSE滤波器。当输入和期望响应均为平稳过程时，它也被称为维纳滤波器。

小结

LMMSE滤波器对观测值和待估参数的分布没有作出任何假设，它约束估计量为观测的线性组合。所需要的先验知识仅仅为观测和待估参数的前两阶矩，因此是一种非参数方法。在实际应用中设计LMMSE滤波器时，为了能够获得数据的二阶矩，我们可以通过大量的数据进行统计来得到二阶矩的估计量，进而应用LMMSE滤波器；或者我们可以假定系统模型和分布来计算二阶矩，但是这种方法就失去了LMMSE的意义。另外的解决方案是使用LSE滤波器。

序贯贝叶斯滤波

在上一部分的LMMSE滤波中，我们没有对观测和待估参数进行建模，这省去了贝叶斯滤波中繁琐的积分计算而以简洁的正则方程进行代替。这种方法在线性条件下是最优的，但是在非线性情况下则性能很差。考虑极端情况，假定观测和待估参数是线性不相关的，则有 $d = 0$ ，此时的LMMSE滤波器甚至不存在。

而且很多时候观测值和待估参数的关系明显可以使用某种模型表达出来，即可以表达为 $y = g (x, d)$ 。其中 $x$ 即为我们需要估计的参数，而 $y$ 通常为我们需要测量的参数，而式中的 $d$ 代表了这个过程中的噪声。如果假定 $x$ 是固定的未知量，那么这个测量方程定义了一个似然函数 $L (y; x)$ ，可以由最大似然估计得到。或者给出 $g$ 的更加确切的形式来求得含参 $y$ 的pdf从而使用一般的经典参数估计方法。实际上，待估参数 $x$ 往往具有一定的物理意义，比如说我们需要定位目标的位置，可以认为目标在匀速、加速，而 $x$ 代表了某种状态。这种状态可以使用一个函数来描述，在离散情况下，状态可以定义为 $x_{n + 1} = f_{n} (x_{n}, v_{n})$ 。该状态方程为 $x$ 的估计提供了一个先验信息，我们称为状态方程。

通过上述讨论可以看到可以通过一个测量方程和一个状态方程把观测和物理过程结合起来。我们称其为状态空间模型（SSM）（应当注意与上一部分在符号上的差异，在推导LMMSE滤波器时我们使用 $x$ 作为观测 $y$ 作为待估参数，而在这一部分状态向量 $x$ 作为待估参数而 $y$ 作为观测）： ${\begin{cases} x_{n + 1} = f_{n} (x_{n}, v_{n}) \\ y_{n} = g_{n} (x_{n}, d_{n}) \end{cases}$ 结合上面的讨论我们看到状态方程提供了贝叶斯估计中需要的先验知识，它通过上一时刻的状态给出这一时刻的状态的先验分布。状态方程定义了状态转移概率 $p (x_{n + 1} | x_{n})$ 。而测量方程提供了似然 $p (y_{n} | x_{n})$ 。而我们希望通过滤波得到的是 $p (x_{n} | Y_{n})$ （有时是全部状态的后验分布 $p (X_{n} | Y_{n} ）$ )，其中 $Y_{n} = y_{1 : n}$ 。接下来我们将通过状态空间模型把贝叶斯滤波推广为序贯滤波的形式。

假定初始分布 $p (x_{0})$ 、噪声的 $v_{n}$ 和 $d_{n}$ 的概率分布是已知的。我们已经得到了上一时刻的后验分布 $p (x_{n - 1} | Y_{n - 1})$ ，那么有 $\begin{aligned} (3) & p (x_{n} | Y_{n}) & = \frac{p (Y_{n} | x_{n}) p (x_{n})}{p (Y_{n})} \\ (4) & = \frac{p (y_{n} | Y_{n - 1}, x_{n}) p (Y_{n - 1} | x_{n}) p (x_{n})}{p (y_{n} | Y_{n - 1}) p (Y_{)}} \\ (5) & = \frac{p (y_{n} | Y_{n - 1}, x_{n}) p (x_{n} | Y_{n - 1}) p (Y_{n - 1}) p (x_{n})}{p (y_{n} | Y_{n - 1}) p (Y_{n - 1}) p (x_{n})} \\ (6) & = \frac{p (y_{n} | x_{n}) p (x_{n} | Y_{n - 1})}{p (y_{n} | Y_{n - 1})} \end{aligned}$ 其中分子第一项即为似然，第二项为先验，表示了给定之前的观测推理出来的现在的知识： $p (x_{n} | Y_{n - 1}) = \int p (x_{n} | x_{n - 1}) p (x_{n - 1} | Y_{n - 1}) d x_{n - 1}$ 一般文献中称计算先验（7）为预测，表示根据过去的观测预测新的状态；称共识（6）为更新，表示观测数据对预测结果的校正。于是序贯贝叶斯就是通过不断地预测–观测–更新估计来得到待估参数的后验分布的。

对于后验分布 $p (X_{n} | Y_{n})$ 的后验分布有递推公式： $p (X_{n} | Y_{n}) = p (X_{n - 1} | Y_{n - 1}) \frac{p (x_{n} | x_{n - 1}) p (y_{n} | x_{n})}{p (y_{n} | Y_{n - 1})}$ 很容易出可以验证： $p (x_{n} | Y_{n}) = \int p (X_{n} | Y_{n}) d X_{n - 1}$ 。

但是很遗憾的是我们以上过程往往并没有闭式解，因为其中的任何概率分布都可能非常难以计算。但是在一类特殊条件下闭式解是存在的，此时的贝叶斯滤波器称为卡尔曼滤波（KF）。

卡尔曼滤波

卡尔曼滤波考虑最一般的SSM情况：线性模型且噪声均为加性高斯白噪声，此时的SSM为 ${\begin{cases} x_{n + 1} = F_{n} x_{n} + d_{n} \\ y_{n} = G_{n} x_{n} + v_{n} \end{cases}$ 其中， $x \in \R^{d_{x}}$ ， $y \in \R^{d_{y}}$ ， $d_{n} \sim N (0, Σ_{d})$ ， $v_{n} \sim N (0, Σ_{v})$ ，各矩阵的维度均满足矩阵相乘的计算要求，且满足以下假设：

噪声为白噪声： $E [d_{n} d_{m}] = Σ_{d} δ_{m n}$ ， $E [v_{n} v_{m}] = Σ_{v} δ_{m n}$ ；
状态和噪声独立： $E [x_{n} d_{m}^{T}] = 0 ， n \leq m$ ， $E [x_{n} v_{m}^{T}] = 0$ ；
测量噪声和模型噪声独立： $E [d_{n} v_{m}^{T}] = 0$ ；
初始分布为高斯分布： $p (x_{0}) \sim N$ 。

由此可知似然函数满足$$p(y_n

x_n)\sim\mathcal{N}(\bf{G}nx_n,\bf{\Sigma}_d) $。因为状态变量的初始分布为高斯分布，噪声也为高斯分布，故所有的先验分布均为高斯的，又似然函数是高斯分布，高斯分布的共轭先验也是高斯分布，故后验分布$ p(x{n-1}

\bf{Y}{n-1}) $仍为高斯分布，其均值即为估计结果，记为$ \hat x{n-1

n-1}

（ 简 记 为

\hat x_{n-1}

） ， 设 其 协 方 差 矩 阵 为

\bf{P}{n-1} $，那么有$ \bf{P}{n-1}=\mathrm{cov}[e_{n-1}]

， 其 中

e_{n-1}=x_{n-1}-\hat x_{n-1}$$为上一时刻的后验误差向量。

经过预测公式（7）的演化得到后验分布为 $p (x_{n} | Y_{n - 1})$ ，因为前一时刻的后验分布和状态转移概率分布均为高斯分布，故此分布仍为高斯分布。可以证明其期望 ${\hat{x}}_{n | n - 1} = F_{n - 1} {\hat{x}}_{n - 1}$ 其方差为 $P_{n | n - 1} = E [(x_{n} - {\hat{x}}_{n | n - 1})^{2}] = E [e_{n | n - 1}^{2}]$ ，它代表了先验误差的协方差矩阵。根据先验误差和上一时刻的后验误差的关系 $e_{n | n - 1} = F_{n - 1} e_{n - 1} + d_{n - 1}$ ，可以得到先验误差的递推公式： $P_{n | n - 1} = F_{n - 1} P_{n - 1} F_{n - 1}^{T} + Σ_{d}$ 将后验分布带入更新公式（6），其中似然 $p (y_{n} | x_{n}) \sim N (G_{n} x_{n}, Σ_{d})$ 。两个高斯分布相乘再除以归一化因子仍为高斯分布，此时分布的指数内部分为： $- \frac{1}{2} (‖ y_{n} - G_{n} x_{n} ‖_{Σ_{v}^{- 1}}^{2} + ‖ x_{n} - {\hat{x}}_{n | n - 1} ‖_{P_{n | n - 1}}^{2})$ 为二次型的形式，为了求其中心值可以对其求偏导数等于零或者对二次型进行化简，最后得到的结果为： $Missing superscript or subscript argument$ 使用矩阵求逆引理得到： ${\hat{x}}_{n} = {\hat{x}}_{n | n - 1} + K_{n} (y_{n} - G_{n} {\hat{x}}_{n | n - 1})$ 其中 $K_{n}$ 为卡尔曼增益矩阵： $K_{n} = P_{n | n - 1} G_{n}^{T} {(G_{n} P_{n | n - 1} G_{n}^{T} + Σ_{v})}^{- 1}$ 此时（10）的物理意义是明显的，公式右侧第一项为先验知识的预测结果；而第二项括号中的部分表示了测量数据和先验预测之间的差异称为新息（innovation），代表了新的观测数据对模型的修正。其中卡尔曼增益矩阵起到了权衡先验和观测的作用。考虑到得到状态的预测时有观测值的预测： $\hat{y} = G_{n} {\hat{x}}_{n | n - 1} + v_{n}$ ，那么预测测量协方差矩阵为： $R_{\hat{y}} = E [y y^{T}] = G_{n} P_{n | n - 1} G_{n}^{T} + Σ_{v}$ 因此卡尔曼增益可以进一步写作： $K_{n} = P_{n | n - 1} G_{n}^{T} R_{\hat{y}}$ 把（10）带入后验误差公式 $e_{n} = x_{n} - {\hat{x}}_{n}$ 并求其协方差矩阵可以得到后验协方差矩阵的递推公式： $P_{n} = (I - K_{n} G_{n}) P_{n | n - 1}$ 至此到下一时刻滤波所需要的全部信息已经被递推地表达了出来。公式（8）到公式（12）构成了卡尔曼滤波的全部内容，其中（8）和（9）对应贝叶斯滤波的预测阶段，公式（9）到（12）为贝叶斯滤波的更新阶段。

从序贯贝叶斯滤波的推导中我们可以看到序贯贝叶斯滤波是一个不断更新概率分布的过程：由之前的先验分布结合状态模型得到先验分布，再结合测量数据得到新的后验分布。卡尔曼滤波限制状态方程和测量方程都是线性加性高斯白噪声的，这保证了不停演化的分布均为高斯分布，因为高斯分布完全由均值和方差决定，因此在贝叶斯滤波中只需要均值和期望在序贯滤波的过ß程中进行更新即可。卡尔曼滤波的PDF演化过程如下：从前一时刻的后验分布$$p(x_{n-1}

\bf{Y}{n-1})\sim\mathcal{N}(\hat x{n-1},\bf{P}_{n-1})

， 线 性 加 性 高 斯 的 状 态 模 型 保 证 了 预 测 的 先 验 分 布 也 是 高 斯 的

p(x_n

\bf{Y}{n-1})\sim\mathcal{N}(\hat x{n

n-1},\bf{P}_{n

n-1})

， 线 性 加 性 高 斯 的 测 量 模 型 保 证 了 似 然 是 高 斯 的 ， 因 此 后 验 分 布 也 是 高 斯 的

p(x_n

\bf{Y}_n)\sim\mathcal{N}(\hat x_n, \bf{P}_n)$$。

卡尔曼滤波取了后验分布的期望作为估计值，因此它是MSE准则下最优的线性估计器（实际上高斯分布的期望中值和众数是一致的，因此他也是MAP准则下最优的）。高斯分布的优良性质保证了其解析的递推公式存在。但是很多情况下，卡尔曼滤波假设的SSM过于理想，因此很多新的贝叶斯滤波方法在非线性非高斯的假设下被提出。

序贯蒙特卡洛滤波

蒙特卡洛滤波是一种特殊的贝叶斯滤波，他应用蒙特卡洛方法来计算不好计算的积分问题。简单的来说，蒙特卡洛方法就是使用随机采样来模拟某个分布，再通过这个分布来估计与此分布相关的量如均值方差等。

假设存在一个我们可以采样但是不知道其具体分布的概率分布 $p (x)$ ，我们可以独立采样 $N_{s}$ 个随机样本 $X^{i} \sim p (x)$ 。当采样点足够多时，根据大数定律该PDF就可以由这些采样点决定： $$ \hat p(x) = \frac{1}{N_s}\sum_{i=1}^{N_s}\delta(x-X^i)

\tag{11} $虽然在每个采样数据处的值均为$ 1/N_s$$，但是其邻域内的样本的疏密足以反映概率密度的大小。以一维标准正态分布为例，按照分布采样出500个样本点，如图所示。在图中标记出了所有样本点的位置（黄色竖线），标准正态分布的PDF以及以0.25位区间长度的统计直方图。可以看到采样的样本点可以很好的近似PDF，当样本点的数量很大时，近似效果更好。

1Dgaussian

该分布对于任意函数的期望于是可以得到： $I^{M C} (φ) = \frac{1}{N_{s}} \sum_{i = 1}^{N_{s}} φ (X^{i})$ 而且蒙特卡洛估计量是无偏估计量： $E [I^{M C}] = \frac{1}{N_{s}} \sum_{i = 1}^{N_{s}} E [φ (X^{i})] = E [φ]$ 蒙特卡洛估计量的方差为： $var [I^{M C}] = \frac{1}{N_{s}} (\int φ^{2} (x) p (x) d x - I^{2} (φ))$ 可以看到蒙特卡洛估计的方差随着样本数量的增加而降低且与 $x$ 的维度无关。但是一个关键的问题是我们经常并不能直接从某个分布中 $p (x)$ 采样。针对这个问题学者提出了很多蒙特卡洛采样方法，常见的有重要性采样（IS）、拒绝采样（RS）、序贯重要性采样（SIS）、重要性重采样（SIR），此外还有一些基于马尔可夫模特卡洛（MCMC）的采样方法。这里主要说明重要性采样方法。

首先考虑关于分布 $p (x)$ 对某个函数 $φ (x)$ 求期望： $E [φ] = \int φ (x) p (x) d x$ 。我们对这个公式做一定的变换可以得到： $E [φ] = \int φ (x) \frac{p (x)}{q (x)} q (x) d x$ 于是上式又可以看作是函数关于 $q (x)$ 求数学期望，而且在每个概率点上都乘上了一个权重因子： $w (x) = p (x) / q (x)$ 。于是对服从 $q (x)$ 分布的数据进行采样再乘以权重因子就可以求得关于分布 $p (x)$ 对于所有函数的期望。从这里我们可以看到：要求出权重 $w$ 的值需要同时知道 $p$ 和 $q$ 的逐点取值。而且 $q$ 的定义域非零部分（支撑集）要包含 $p$ 的。我们称便于采样的分布 $q$ 为建议分布，对于任何概率分布 $p$ 假设可以表示成未归一化的形式： $p (x) = \frac{γ (x)}{Z} = \frac{w (x) q (x)}{Z}$ 其中 $w (x) = \frac{γ (x)}{q (x)}$ 为（未归一）权重函数。 $Z = \int w (x) q (x) d x$ 为归一化系数。现在从简易分布 $q$ 中采出 $N_{s}$ 个样本 $X^{i}$ ，那么归一化系数 $Z$ 的估计根据之前的讨论即为： $\hat{Z} = \frac{1}{N_{s}} \sum_{i = 1}^{N_{s}} w (X^{i})$ 根据蒙特卡洛PDF估计公式可以得到： $\hat{p} (x) = \frac{1}{N_{s}} \sum_{i = 1}^{N_{s}} \frac{w (X^{i})}{\hat{Z}} δ_{X^{i}} (x) = \sum_{i = 1}^{N_{s}} W^{i} δ_{X^{i}} (x)$ 其中 $W$ 为归一化的权重，公式为： $W^{i} = \frac{w (X^{i})}{\sum_{j = 1}^{N_{s}} w (X^{i})}$ 此时使用样本来表示PDF的估计时，除了采样点的疏密程度以外，还有用来表示PDF大小的权值。这里以二维的高斯分布来逼近均匀分布作为示例。点的大小表示权重的大小。左图为均匀分布的范围和从高斯分布中采样得到的数据，右图为乘以权重后的采样点。可以看到越靠近中心的点的权重越小，而越靠近边缘的点权重越大。

2Dcompare

建议分布 $q$ 的选取对于IS来说至关重要。当 $p$ 和 $q$ 成正比时此方法的方差最小。当 $p$ 和 $q$ 的PDF的模式差别很大时，IS的性能很差，考虑极端情况：使用均值在-100处的标准正态分布区逼近均值在100处的标准正态分布，从第一个PDF中产生的采样点几乎不会在第二个PDF的高概率处出现，此时的IS性能非常差。

粒子滤波

粒子滤波是基于模特卡洛原理的贝叶斯序贯滤波方法，它使用一系列的粒子 ${x_{0 : n}^{i}, w_{n}^{i}}_{i = 1}^{N_{p}}$ 来模拟概率分布的PDF（在粒子滤波部分不再以大写字母标注样本，而是以小写字母带右上角索引角标来标注，仍然有大写加粗的字母为从0时刻累计到当前时刻的所有变量），其中 $n$ 代表了粒子更新的时刻， $i$ 为各个粒子的索引， $N_{s}$ 为所有的粒子数目。它和一般的贝叶斯滤波一样，期望求得的信息是后验概率 $p (X_{n} | Y_{n})$ 。根据上一部分IS方法的描述，后验PDF可以表示为： $\hat{p} (x_{0 : n} | Y_{n}) = \sum_{i = 1}^{N_{p}} w_{n}^{i} δ (x_{0 : n} - x_{0 : n}^{i})$

于是我们可以选取合适的建议分布 $q (x_{0 : n})$ 并从中抽取粒子并计算权重。根据之前的讨论，建议分布应当尽可能的和真实的分布成正比，因此我们应当尽可能的利用所给信息。所以建议分布最好应用上观测数据带来的信息 $q (x_{0 : n} | Y_{n})$ 。此时粒子的权重为： $w_{n}^{i} \propto \frac{p (x_{0 : n}^{i} | Y_{n})}{q (x_{0 : n}^{i} | Y_{n})}$ 注意到计算权重需要我们知道后验概率的逐点取值，但这是我们无法求得的，我们只有上一个阶段的后验概率估计，也就是上一阶段的粒子及其权重。为了计算当前粒子的权重，我们需要利用上一阶段的粒子权重，而上一时刻的粒子是从分布 $q (x_{0 : n - 1} | Y_{n})$ 中采样得到的，通过一个条件概率可以把上一时刻的建议分布和当前时刻的建议分布联系起来： $q (x_{0 : n} | Y_{n}) = q (x_{n} | x_{0 : n - 1}, Y_{n}) q (x_{0 : n - 1} | Y_{n}) = q (x_{n} | x_{0 : n - 1}, Y_{n}) q (x_{0 : n - 1} | Y_{n - 1})$ 同理，后验分布也有递推关系： $\begin{aligned} p (X_{n} | Y_{n}) & = p (X_{n - 1} | Y_{n - 1}) \frac{p (x_{n} | x_{n - 1}) p (y_{n} | x_{n})}{p (y_{n} | Y_{n - 1})} \\ \propto p (y_{n} | x_{n}) p (x_{n} | x_{n - 1}) p (x_{0 : n - 1} | y_{1 : n - 1}) \end{aligned}$ 那么可以得到权重的递推公式： $\begin{aligned} w_{k}^{i} & \propto \frac{p (y_{n} | x_{n}^{i}) p (x_{n}^{i} | x_{n - 1}^{i}) p (x_{0 : n - 1}^{i} | y_{1 : n - 1})}{q (x_{n}^{i} | x_{0 : n - 1}^{i}, Y_{n}) q (x_{0 : n - 1}^{i} | Y_{n - 1})} \\ = w_{k - 1}^{i} \frac{p (y_{n} | x_{n}^{i}) p (x_{n}^{i} | x_{n - 1}^{i})}{q (x_{n}^{i} | x_{0 : n - 1}^{i}, y_{0 : n - 1})} \end{aligned}$ 于是，在序贯蒙特卡洛滤波中我们可以调整的便是建议分布 $q (x_{n} | x_{0 : n - 1}^{i}, y_{0 : n - 1})$ 的选择，不同的建议分布的选择产生了各种粒子滤波器之间的差异。其中最常用的SIR粒子滤波器就采用了最简捷最直观的建议分布： $q (x_{n} | x_{0 : n - 1}^{i}, y_{0 : n - 1}) = p (x_{n} | x_{n - 1}^{i})$ 此时有权重更新公式： $w_{k}^{i} \propto w_{k - 1}^{i} p (y_{n} | x_{n}^{i})$ 通过归一化后得到新的粒子权重 $W_{k}^{i}$ 。到目前为止已经完成了贝叶斯框架下的蒙特卡洛序贯滤波，但是这还不足以完成SIR粒子滤波。因为此时的粒子会出现退化现象：只有少量粒子占有很大的权重而其余的粒子的权重小到几乎可以忽略。为了解决这一问题需要引入重采样（Resampling）：从给出的 $N_{s}$ 个粒子中按照其其权重的大小为概率有放回的抽取 $N_{s}$ 次，得到的新的等权重的粒子。这一过程可以看作有权重的Bootstrap，因此SIR滤波器也被称为Bootstrap滤波器。综上，SIR粒子滤波的过程可以被概括为：

预测：根据上一时刻的粒子采样出新的粒子

x_{n}^{i} \sim p (x_{n} | x_{n - 1}^{i})

这一过程可以由系统的状态方程得到： $x_{n}^{i} = f_{n} (x_{n}^{i}, v_{n})$

更新：根据新的观测 $y_{n}$ 更新粒子权重，其中由于重采样上一时刻的粒子权重均为 $1 / N_{s}$ 。可以看到每个粒子的权重是正比于其似然的，这是一个很直观的结果，似然大的粒子拥有更大的权重。

W_{k}^{i} = \frac{p (y_{n} | x_{n}^{i})}{\sum_{i = 1}^{N_{s}} p (y_{n} | x_{n}^{i})}

重采样：从已经得到的粒子重采样得到新的等权重粒子并通过新的粒子表示后验PDF。

{x_{n}^{i}, W_{k}^{i}} \overset{R e s a m p l i n g}{\to} {x_{n}^{j}, 1 / N_{s}} p (x_{n} | Y_{n}) \approx \frac{1}{N_{s}} \sum_{j = 1}^{N_{s}} δ (x_{n} - x_{n}^{j})

集综卡尔曼EnKF

考虑卡尔曼滤波中的先验和后验协方差矩阵 $P_{k | k - 1}$ 和 $P_{k}$ ，假如我们知道 $k$ 时刻真实的状态值 $x_{k}^{t}$ 。那么就有 $\bf P_{k|k-1}=\mathbb{E} [(x_{k|k-1}-x_k^t)(x_{k|k-1}-x_k^t)^T]\ \bf P_{k}=\mathbb{E} [(x_{k}-x_k^t)(x_{k}-x_k^t)^T]$ 但是真实状态是未知的，因此定义集综误差矩阵

A Brief Overview on the Framework of Bayesian Filtering