16 粒子滤波

16 粒子滤波

16.1 背景介绍

16.1.1 Particle Filter是什么？

Dynamic Model包含：

• HMM——关注Decoding问题
• Linear Dynamic System——关注Filtering问题
• Patricle Filter——Nan-Linear，Nan-Gauss，关注Filtering问题

16.1.2 Patricle Filter的状态如何转移？

在HMM中有$\lambda =(\pi ,A,B)$，用于表示状态转移矩阵和发射矩阵。

由于Linear Dynamic System和Particle Filter中的隐变量与观测变量连续，状态转移矩阵和发射矩阵不用矩阵A、B表示，表示为：
$$\begin{array}{rl}{Z}_t& =g(Z_{t-1},u,\epsilon )\mapsto A\\ X_t& =h(Z_t,u,\delta )\mapsto B\end{array}$$
在Kalman Filter中我们假设以上两个公式均为线性，且噪声为Gauss。表示为：
$$\begin{array}{rlr}{Z}_t& =A\cdot Z_{t-1}+B+\epsilon & \epsilon \backsim N(0,Q)\\ X_t& =C\cdot Z_t+D+\delta & \delta \backsim N(0,R)\end{array}$$

回顾：Kalman Filter通过预测+更新的方式求解Filtering问题

Step1: 求解Prediction问题
$$P(Z_t|x_1,\dots ,x_{t-1})={\int }_{Z_{t-1}}P(Z_t|Z_{t-1})\cdot P(Z_{t-1}|x_1,\dots ,x_{t-1}){\mathrm{d}}_{Z_{t-1}}$$
Step2: 求解update问题
$$P(Z_t|x_1,\dots ,x_t)\propto P(X_t|Z_t)\cdot P(Z_t|x_1,\dots ,x_{t-1})$$
具体可以通过条件概率相互求解的公式求解。

而在Particle Filter中转移方程非线形非高斯，只能通过采样的方式求解。

16.1.3 如何通过采样求解Particle Filter

由于转移方程非线性非高斯，所以只能采取近似方法求解Filtering问题。这里使用Monte Carlo Method，通过采样求取期望：
$$P(Z|X)\to E_{Z|X}[f(z)]={\int }_{Z}f(z)\cdot P(Z|X)\mathrm{d}Z\approx \frac{1}{N}\sum _{i=1}^{N}f(Z^{(i)})$$
其中$Z^{(i)}$为样本，且$Z^{(1)},Z^{(2)},\dots ,Z^{(N)}\backsim P(Z|X)$。

16.2 重要性采样

16.2.1 重要性采样方法

已知问题：
$$E[f(Z)]=\int f(Z)p(Z)\mathrm{d}Z$$
求解方法：

1. 但$p(Z)$的分布复杂，无法直接采样，所以我们引入已知分布$q(Z)$，$q(Z)$也称为提议分布（Proposed dist）：
   $$\begin{array}{rl}E[f(Z)]& =\int f(Z)p(Z)\mathrm{d}Z\\ & =\int f(Z)\cdot \frac{p(Z)}{q(Z)}\cdot q(Z)\mathrm{d}Z\\ & =\frac{1}{N}\sum _{i=1}^{N}f(Z^{(i)})\cdot \frac{p(Z)}{q(Z)}\end{array}$$

2. 其中$\frac{p(Z)}{q(Z)}$被称为weight，表示为$w^{(i)}$，用于表示提议分布与实际分布之间的相似度：
   $$E[f(Z)]=\frac{1}{N}\sum _{i=1}^{N}f(Z^{(i)})\cdot w^{(i)}$$

所以我们通过采样可以求出$f(Z^{(i)})$，然后我们的目标就是求出对应的$w^{(i)}$。

16.2.2 Sequential Importance Sampling

引入SIS的原因：

• 由于Filtering问题在递推过程中求解的是$P(Z_t|X_{1:t})$，所以对应就会有$w_t^{(i)}=\frac{P(Z_t|X_{1:t})}{q(Z_t^{(i)}|X_{1:t})}$，但是随着$t$增加，每次$w$都要求$N$遍，时间开销大。所以引入Sequential Importance Sampling，通过递推的方式求解$w$（$w_{t-1}^{(i)}\to w_t^{(i)}$）。

推导过程：

1. 已知：
   $$w_t^{(i)}\propto \frac{P(Z_{1:t}|X_{1:t})}{q(Z_{1:t}|X_{1:t})}$$

2. 分解$P(Z_{1:t}|X_{1:t})$，其中将已知量（只由观测变量构成的数据）假设为常数：
   $$\begin{array}{rl}P(Z_{1:t}|X_{1:t})& =\frac{P(Z_{1:t},X_{1:t})}{P(X_{1:t})}\\ & =\frac{1}{C}\cdot P(X_t|Z_{1:t},X_{1:t-1})\cdot P(Z_t,X_{1:t-1})\\ & =\frac{1}{C}\cdot P(X_t|Z_t)\cdot P(Z_t|Z_{1:t-1},X_{1:t-1})\cdot P(Z_{1:t-1},X_{1:t-1})\\ & =\frac{1}{C}\cdot P(X_t|Z_t)\cdot P(Z_t|Z_{t-1})\cdot P(Z_{1:t-1},X_{1:t-1})\\ & =\frac{1}{C}\cdot P(X_t|Z_t)\cdot P(Z_t|Z_{t-1})\cdot P(Z_{1:t-1}|X_{1:t-1})\cdot P(X_{1:t-1})\\ & =\frac{D}{C}\cdot P(X_t|Z_t)\cdot P(Z_t|Z_{t-1})\cdot P(Z_{1:t-1}|X_{1:t-1})\end{array}$$
   通过以上推导可将$P(Z_{1:t}|X_{1:t})$分解为由$P(Z_{1:t-1}|X_{1:t-1})$组成的公式

3. 分解$q(Z_{1:t}|X_{1:t})$：
   $$\begin{array}{rl}q(Z_{1:t}|X_{1:t})& =q(Z_t|Z_{1:t-1},X_{1:t})\cdot q(Z_{1:t-1}|X_{1:t})\\ & =q(Z_t|Z_{1:t-1},X_{1:t})\cdot q(Z_{1:t-1}|X_{1:t-1})\end{array}$$

4. 结合起来就是：
   $$\begin{array}{rl}{w}_t^{(i)}& \propto \frac{P(Z_{1:t}|X_{1:t})}{q(Z_{1:t}|X_{1:t})}\\ & \propto \frac{P(X_t|Z_t)\cdot P(Z_t|Z_{t-1})\cdot P(Z_{1:t-1}|X_{1:t-1})}{q(Z_t|Z_{1:t-1},X_{1:t})\cdot q(Z_{1:t-1}|X_{1:t-1})}\\ & \propto \frac{P(X_t|Z_t)\cdot P(Z_t|Z_{t-1})}{q(Z_t|Z_{1:t-1},X_{1:t})}\cdot w_{t-1}^{(i)}\\ & \propto \frac{P(X_t|Z_t)\cdot P(Z_t|Z_{t-1})}{q(Z_t|Z_{t-1},X_{1:t})}\cdot w_{t-1}^{(i)}\end{array}$$

具体可以表示为一个算法：

条件：t-1时刻的采样已完成$\to w_{t-1}^{(i)}$已知。

t时刻:
for i = 1 to N:
$Z_t^{(i)}\backsim q(Z_t|Z_{t-1},X_{1:t})$ // 采样
$w_t^{(i)}\propto \frac{P(X_t|Z_t)\cdot P(Z_t|Z_{t-1})}{q(Z_t|Z_{t-1},X_{1:t})}\cdot w_{t-1}^{(i)}$ // 计算
end

$w_t^{(i)}$要归一化，${\sum }_{i=1}^N{w}_t^{(i)}$

但是通过SIS直接求解有一个问题：$w_t^{(i)}$的权值会退化——有的大有的小，随着维度上升，可能会出现如：$w_t^{(1)}\to 1$但$w_t^{(N)}\to 0$的情况。解决方案有：

1. Resampling——重采样（通过别的方法重新采样）
2. 选择一个合适的proposed dist q(Z)

16.2.3 Resampling

这里介绍一种最简单的重采样方法。

倘若第一遍的采样结果为第二列：

数据编号	权重（weight）	pdf	cdf
$x^{(1)}$	0.1	0.1	0.1
$x^{(2)}$	0.1	0.1	0.2
$x^{(3)}$	0.8	0.8	1

我们将权重假设为当前数据的概率，通过权重建立概率密度函数，并求出其分布函数，即可通过分段函数进行采样。

这样的优点是可以将数据集中在权重大的地方。

16.2.4 采样总结——Basic Particle Filter

结合：重要性采样方法+SIS+Resampling，就是简单的粒子滤波求解方案：Basic Particle Filter

16.3 具体算法——SIR Filter

Particle Filter整体就是通过每个时刻的采样与迭代地计算权重，通过Monte Carlo方法预测的方法。

根据16.2.2已知迭代公式为：
$$w_t^{(i)}\propto \frac{P(X_t|Z_t)\cdot P(Z_t|Z_{t-1})}{q(Z_t|Z_{t-1},X_{1:t})}\cdot w_{t-1}^{(i)}$$
其中我们令$q(Z_t|Z_{t-1},X_{1:t})$为用于采样的分布，我们假设采样的分布就是状态转移函数：
$$q(Z_t|Z_{t-1},X_{1:t})=p(Z_t|Z_{t-1}^{(i)})$$
可以简化计算，算法可以总结为"generate and test"：

1. generate表示采样：采样的方式变成了：
   $$Z_t^{(i)}\backsim q(Z_t|Z_{t-1},X_{1:t})⟹Z_t^{(i)}\backsim p(Z_t|Z_{t-1}^{(i)})$$

2. test表示通过权重的迭代计算进行预测：变成了：
   $$w_t^{(i)}\propto P(X_t|Z_t^{(i)})\cdot w_{t-1}^{(i)}$$

上面的方法总结下来就是：SIR Filter（Sampling-Importance-Resampling）——SIS + Resampling + （$q(Z_t|Z_{t-1},X_{1:t})=p(Z_t|Z_{t-1}^{(i)})$）