文章目录
 
线性动态系统
模型定义
LDS 滤波
递推公式
Prediction
Update
   
   
参考文献
 

 
线性动态系统
 
动态系统两大问题：
 $$\{\begin{array}{rl}& Learning\\ & Inference:p(Z|X)\{\begin{array}{rl}& decoding\to HMM{}_{(\text{维特比})}\\ & probofevidence\to p(O|\lambda ){}_{(\text{前向/后向})}\\ & filtering:p(z_t|x_1,x_2,\cdots ,x_t)\\ & smoothing:p(z_t|x_1,x_2,\cdots ,x_T){}_{(\text{给定所有数据})}\\ & prediction:\{\begin{array}{r}p(z_t|x_1,x_2,\cdots ,x_{t-1})\\ p(x_t|x_1,x_2,\cdots ,x_{t-1})\end{array}\end{array}\end{array}$$
 
模型定义
 
HMM 模型适用于隐变量是离散的值的时候，对于连续隐变量的 HMM，常用线性动态系统描述线性高斯模型的状态变量，使用粒子滤波来表述非高斯非线性变量。
 
 
LDS 又叫卡尔曼滤波，其线性体现在上一时刻与当前时刻隐变量以及隐变量与观测变量之间：
 $$\begin{array}{rl}& \begin{array}{rl}{z}_t& =A\cdot z_{t-1}+B+\epsilon \\ x_t& =C\cdot z_t+D+\delta \end{array}\}{}_{\text{线性+噪声}}\\ & \begin{array}{rl}\epsilon & \sim \mathcal{N}(0,Q)\\ \delta & \sim \mathcal{N}(0,R)\end{array}\}{}_{\text{高斯分布噪声}}\end{array}$$
 
由线性高斯模型的性质可知：
 $$\begin{array}{rl}p(z_t|z_{t-1})& \sim \mathcal{N}(A\cdot z_{t-1}+B,Q)\\ p(x_t|z_t)& \sim \mathcal{N}(C\cdot z_t+D,R)\\ z_1& \sim \mathcal{N}({\mu }_1,{\Sigma }_1)\end{array}$$
 
LDS 滤波
 
递推公式
 
LDS 更关心滤波问题：边缘后验概率分布 $p(z_t|x_1,x_2,\cdots ,x_t)$。类似 HMM 中的前向算法，可以找一个递推关系
 
 
令 $x_1,x_2,\cdot ,x_t$ 的简化表示为 $x_{1:t}$ ,有：
 $$\begin{array}{rl}p(z_t|x_{1:t})& =\frac{p(x_{1:t},z_t)}{\sout{p(x_{1:t})}}\\ & \propto p(x_{1:t},z_t)=p(x_t|x_{1:t-1},z_t)p(x_{1:t-1},z_t)\\ & =p(x_t|z_t)p(x_{1:t-1},z_t)\\ & =p(x_t|z_t)\underset{Prediction}{\underbrace{p(z_t|x_{1:t-1})}}\underset{\text{常数}}{\underbrace{p(x_{1:t-1})}}\\ & \propto p(x_t|z_t)p(z_t|x_{1:t-1})\end{array}$$
 
对于上式中的 $prediction$ 项：
 $$\begin{array}{rl}p(z_t|x_{1:t-1})& ={\int }_{z_{t-1}}p(z_t,z_{t-1}|x_{1:t-1})d{z}_{t-1}\\ & ={\int }_{z_{t-1}}p(z_t|z_{t-1},x_{1:t-1})p(z_{t-1}|x_{1:t-1})d{z}_{t-1}\\ & ={\int }_{z_{t-1}}p(z_t|z_{t-1})p(z_{t-1}|x_{1:t-1})d{z}_{t-1}\end{array}$$
 
至此，得到了关于 $p(z_t|x_{1:t})$ 的递推公式，总体上可分为两个部分
 
$p(z_{t-1}|x_{1:t-1})$ 可看作 $update$
$p(z_t|x_{1:t-1})$ 可看作 $pridiction$
 
这其实是一个 Online 的过程：
 $$\begin{array}{rl}& t=1\{\begin{array}{rl}& p(z_1|x_1)⟶update\\ & p(z_2|x_1)⟶pridiction\end{array}\\ & t=2\{\begin{array}{rl}& p(z_2|x_1,x_2)⟶update\\ & p(z_3|x_1,x_2)⟶pridiction\end{array}\\ & ⋮\\ & t=T\{\begin{array}{rl}& p(z_T|x_1,\cdots ,x_T)⟶update\\ & p(z_{T+1}|x_1,\cdots ,x_T)⟶pridiction\end{array}\end{array}$$
 
Prediction
 
对于线性高斯假设，可以得到解析解：
 
$$\begin{array}{rl}p(z_t|x_{1:t-1})& ={\int }_{z_{t-1}}p(z_t|z_{t-1})p(z_{t-1}|x_{1:t-1})d{z}_{t-1}\\ & ={\int }_{z_{t-1}}\mathcal{N}(A{z}_{t-1}+B,Q)\mathcal{N}({\mu }_{t-1},{\Sigma }_{t-1})d{z}_{t-1}\end{array}$$
 
其中第二个高斯分布是上一步的 Update 过程，所以根据线性高斯模型性质：
 $$\boxed{\begin{array}{rl}p(x)& =\mathcal{N}(x|\mu ,{\Lambda }^{-1})\\ p(y|x)& =\mathcal{N}(y|Ax+b,L^{-1})\\ p(y)& =\mathcal{N}(y|A\mu +b,L^{-1}+A{\Lambda }^{-1}{A}^T)\end{array}}$$
 
可直接写出该积分： $$p(z_t|x_{1:t-1})=\mathcal{N}(A{\mu }_{t-1}+B,Q+A{\Sigma }_{t-1}{A}^T)$$
 
Update
 
$$p(z_t|x_{1:t})\propto p(x_t|z_t)p(z_t|x_{1:t-1})$$ 同样可以利用线性高斯模型的性质直接写出这个高斯分布。
 
参考文献
 
【1】线性动态系统
 【2】详解卡尔曼滤波原理
 【3】How a Kalman filter works, in pictures
 【4】线性动态系统与卡尔曼滤波
 【5】隐马尔科夫模型（HMM）与线性动态系统（LDS）
 【6】贝叶斯视角下的卡尔曼滤波
 【7】机器学习中的高斯分布
 【8】[易懂]如何理解那个把嫦娥送上天的卡尔曼滤波算法Kalman filter?

一、概述
 
HMM 模型适用于隐变量是离散的值的时候，对于连续隐变量的 HMM，常用卡尔曼滤波（Kalman Filtering）描述线性高斯模型的态变量，使用粒子滤波（Particle Filter）来表述非高斯非线性的态变量。
 
线性体现在上一时刻和这一时刻的隐变量以及隐变量和观测变量之间，它们的关系可以表示为：
 
 
 
$$\begin{gathered} z_t=A\cdot z_{t-1}+B+\epsilon \\ {x}_t=C\cdot z_t+D+\delta \\ \epsilon \sim N(0,Q) \\ \delta \sim N(0,R) \end{gathered}$$
 
 
类比HMM中几个参数，我们也可以写出类似初始概率、转移概率或发射概率的形式：
 
 
 
$$\begin{gathered} P(z_t|z_{t-1})\sim N(A\cdot z_{t-1}+B,Q) \\ P(x_t|z_t)\sim N(C\cdot z_t+D,R) \\ {z}_1\sim N({\mu }_1,{\Sigma }_1) \end{gathered}$$
 
 
所有的参数为：
 
 
 
$$\theta =(A,B,C,D,Q,R,{\mu }_1,{\Sigma }_1)$$
 
 
二、Filtering问题
 
在多个inference问题中，卡尔曼滤波更关心Filtering问题，即求边缘概率：
 
 
 
$$P(z_t|x_1,x_2,\cdots ,x_t)$$
 
 
类似HMM的前向算法，我们需要找到一个递推关系：
 
 
 
$$\begin{gathered} P(z_t|x_1,x_2,\cdots ,x_t) \\ =\frac{P(x_1,x_2,\cdots ,x_t,z_t)}{P(x_1,x_2,\cdots ,x_t)} \\ \propto P(x_1,x_2,\cdots ,x_t,z_t) \\ =\underset{P(x_t|z_t)}{\underbrace{P(x_t|x_1,x_2,\cdots ,x_{t-1},z_t)}}\cdot P(x_1,x_2,\cdots ,x_{t-1},z_t) \\ =P(x_t|z_t)\cdot P(x_1,x_2,\cdots ,x_{t-1},z_t) \\ =P(x_t|z_t)\cdot \underset{prediction问题}{\underbrace{P(z_t|x_1,x_2,\cdots ,x_{t-1})}}\cdot P(x_1,x_2,\cdots ,x_{t-1}) \\ \propto P(x_t|z_t)\cdot P(z_t|x_1,x_2,\cdots ,x_{t-1}) \end{gathered}$$
 
 
上式结果中，$P(x_t|z_t)$已知，而另一项可做以下转化：
 
 
 
$$\begin{gathered} P(z_t|x_1,x_2,\cdots ,x_{t-1}) \\ ={\int }_{z_{t-1}}P(z_t,z_{t-1}|x_1,x_2,\cdots ,x_{t-1})\mathrm{d}{z}_{t-1} \\ ={\int }_{z_{t-1}}\underset{P(z_t|z_{t-1})}{\underbrace{P(z_t|z_{t-1},x_1,x_2,\cdots ,x_{t-1})}}\cdot \underset{Filtering问题}{\underbrace{P(z_{t-1}|x_1,x_2,\cdots ,x_{t-1})}}\mathrm{d}{z}_{t-1} \\ ={\int }_{z_{t-1}}P(z_t|z_{t-1})\cdot P(z_{t-1}|x_1,x_2,\cdots ,x_{t-1})\mathrm{d}{z}_{t-1} \end{gathered}$$
 
 
因此，我们找到了Filtering问题的递推式：
 
 
 
$$P(z_t|x_1,x_2,\cdots ,x_t)=C\cdot P(x_t|z_t)\cdot {\int }_{z_{t-1}}P(z_t|z_{t-1})\cdot P(z_{t-1}|x_1,x_2,\cdots ,x_{t-1})\mathrm{d}{z}_{t-1}$$
 
 
因此，我们可以确定求解Filtering问题的步骤如下：
 
 
 
$$\begin{gathered} t=1\{\begin{array}{c}P(z_1|x_1)\to update\\ P(z_2|x_1)\to prediction\end{array} \\ t=2\{\begin{array}{c}P(z_2|x_1,x_2)\to update\\ P(z_3|x_1,x_2)\to prediction\end{array} \\ ⋮ \\ t\{\begin{array}{c}P(z_t|x_1,x_2,\cdots ,x_t)\to update\\ P(z_{t+1}|x_1,x_2,\cdots ,x_t)\to prediction\end{array} \end{gathered}$$
 
 
很明显这是一个online的过程。
 
三、Filtering问题求解
 
通过上述转化我们可以确定Filtering问题的计算是通过以下两步迭代计算进行的：
 
 
 
Step1 Prediction:
 $$P(z_t|x_1,x_2,\cdots ,x_{t-1})={\int }_{z_{t-1}}P(z_t|z_{t-1})\cdot P(z_{t-1}|x_1,x_2,\cdots ,x_{t-1})\mathrm{d}{z}_{t-1}$$
 Step2 Update:
 $$P(z_t|x_1,x_2,\cdots ,x_t)=C\cdot P(x_t|z_t)\cdot P(z_t|x_1,x_2,\cdots ,x_{t-1})$$
 
 
我们可以确定的是几个高斯分布经过相乘或者积分运算后仍然是高斯分布，所以我们假设：
 
 
 
$$\begin{gathered} Prediction:P(z_t|x_1,x_2,\cdots ,x_{t-1})=N(z_t|{\mu }_t^{\ast },{\Sigma }_t^{\ast }) \\ Update:P(z_t|x_1,x_2,\cdots ,x_t)=N(z_t|{\mu }_t,{\Sigma }_t) \end{gathered}$$
 
 
代入高斯分布的形式可以得到：
 
 
 
$$\begin{gathered} Prediction:N(z_t|{\mu }_t^{\ast },{\Sigma }_t^{\ast })={\int }_{z_{t-1}}N(z_t|A\cdot z_{t-1}+B,Q)\cdot N(z_{t-1}|{\mu }_{t-1},{\Sigma }_{t-1})\cdot \mathrm{d}{z}_{t-1} \\ Update:N(z_t|{\mu }_t,{\Sigma }_t)=C\cdot N(x_t|C\cdot z_t+D,R)\cdot N(z_t|{\mu }_t^{\ast },{\Sigma }_t^{\ast }) \end{gathered}$$
 
 
接下来的求解需要用到高斯分布|机器学习推导系列（二）第六部分内容中我们得到的结论，即已知$P(x)$和$P(y|x)$来求$P(y)$和$P(x|y)$，这里我们直接套用公式即可。
 
首先，在Prediction过程中：
 
 
 
$$\underset{P(y)}{\underbrace{P(z_t|x_1,x_2,\cdots ,x_{t-1})}}={\int }_{z_{t-1}}\underset{P(y|x)}{\underbrace{P(z_t|z_{t-1})}}\cdot \underset{P(x)}{\underbrace{P(z_{t-1}|x_1,x_2,\cdots ,x_{t-1})}}\mathrm{d}{z}_{t-1}$$
 
 
代入计算$P(y)$的公式可得：
 
 
 
$$\begin{gathered} {\mu }_t^{\ast }=A{\mu }_{t-1}+B \\ {\Sigma }_t^{\ast }=Q+A{\Sigma }_{t-1}{A}^T \end{gathered}$$
 
 
在update过程中：
 
 
 
$$\underset{P(x|y)}{\underbrace{P(z_t|x_1,x_2,\cdots ,x_t)}}=C\cdot \underset{P(y|x)}{\underbrace{P(x_t|z_t)}}\cdot \underset{P(x)}{\underbrace{P(z_t|x_1,x_2,\cdots ,x_{t-1})}}$$
 
 
代入计算$P(x|y)$的的公式也可以得出结果，过程比较复杂，所以省略。
 
注意这里将$x_1,x_2,\cdots ,x_{t-1}$看做已知即可，然后再套用$p(x)$、$p(y|x)$等形式。
 
另外这里只需要根据公式直接套用得出结果即可，具体的得出的Prediction和Update计算公式并没有用到。

 

  

    
协方差矩阵
状态协方差矩阵传递
状态协方差的更新
Matlab 实现
Matlab效果
测试代码
测试效果
 




 

 【原创】Liu_LongPo 转载请注明出处 
【CSDN】http://blog.csdn.net/llp1992


卡尔曼滤波器是一种利用线性系统状态方程，通过系统输入输出观测数据，对系统状态进行最优估计的算法。而且由于观测包含系统的噪声和干扰的影响，所以最优估计也可看做是滤波过程。
卡尔曼滤波器的核心内容就是5条公式，计算简单快速，适合用于少量数据的预测和估计。
下面我们用一个例子来说明一下卡尔曼算法的应用。
假设我们想在有一辆小车，在 t 时刻其速度为 Vt ，位置坐标为 Pt，ut 表示 t 时刻的加速度，那么我们可以用Xt表示 t 时刻的状态，如下：


 

 


则我们可以得到，由t-1 时刻到 t 时刻，位置以及速度的转换如下：


 

 


用向量表示上述转换过程，如下：


 

 


如下图：


 

 


那么我们可以得到如下的状态转移公式：


 

  (1)


其中矩阵 F 为状态转移矩阵，表示如何从上一状态来推测当前时刻的状态，B 为控制矩阵，表示控制量u如何作用于当前矩阵，上面的公式 x 有顶帽子，表示只是估计值，并不是最优的。
有了状态转移公式就可以用来推测当前的状态，但是所有的推测都是包含噪声的，噪声越大，不确定越大，协方差矩阵用来表示这次推测带来的不确定性
协方差矩阵
假设我们有一个一维的数据，这个数据每次测量都不同，我们假设服从高斯分布，那么我们可以用均值和方差来表示该数据集，我们将该一维数据集投影到坐标轴上，如下图：


 

 


可以看到，服从高斯分布的一维数据大部分分布在均值附近。
现在我们来看看服从高斯分布的二维数据投影到坐标轴的情况，如下图：


 

 


二维数据比一维数据稍微复杂一点，投影后有3种情况，分别是： 
左图：两个维的数据互不相关； 
中图：两个维的数据正相关，也就是 y 随着 x 的增大而增大（假设两个维分别为 x 和 y） 
右图：两个维的数据负相关，也就是 y 随着 x 的增大而减小。
那怎么来表示两个维的数据的相关性呢？答案就是协方差矩阵。


 

 



状态协方差矩阵传递
在公式（1）之中，我们已经得到了状态的转移公式，但是由上面可知，二维数据的协方差矩阵对于描述数据的特征是很重要的，那么我们应该如何更新或者说传递我们的二维数据的协方差矩阵呢？假如我们用 P 来表示状态协方差，即


 

 


那么加入状态转换矩阵 F ，得到


 

  (2)


也即：


 

 


因此我们便得到了协方差的转换公式。
现在我们得到了两个公式，运用这两个公式能够对现在状态进行预测。按照我们的正常思路来理解，预测结果不一定会对嘛，肯定有误差。而且在我们大多数回归算法或者是拟合算法中，一般思路都是先预测，然后看看这个预测结果跟实际结果的误差有多大，再根据这个误差来调整预测函数的参数，不断迭代调整参数直到预测误差小于一定的阈值。
卡尔曼算法的迭代思想也类似，不过这里根据误差调整的是状态 X 。
在这里，我们的实际数据就是 Z， 如下图：


 

 


其中，矩阵 H 为测量系统的参数，即观察矩阵，v 为观测噪声， 其协方差矩阵为R
那么我们的状态更新公式如下：


 

 


其中K 为卡尔曼系数， Z-Hx 则为残差，也就是我们说的，预测值与实际值的误差。
K的作用：

 
1.K 权衡预测协方差P和观察协方差矩阵R那个更加重要，相信预测，残差的权重小，相信观察，残差权重大，由 K 的表达是可以退出这个结论 
2,将残差的表现形式从观察域转换到状态域（残差与一个标量，通过K 转换为向量），由 状态 X 的更新公式可得到该结论。

至此，我们已经得到了 t 状态下的最优估计值 Xt。但为了能让我们的迭代算法持续下去，我们还必须更新状态协方差的值。
状态协方差的更新


 

 


以上就是卡尔曼滤波算法的思想，只有简单的 5 条公式，总结如下：


 

 


Matlab 实现
function kalmanFiltering
%%
clc
close all

%%
%
%   Description : kalmanFiltering
%   Author : Liulongpo
%   Time：2015-4-29 16:42:34
%

%%
Z=(1:2:200); %观测值  汽车的位置  也就是我们要修改的量
noise=randn(1,100); %方差为1的高斯噪声
Z=Z+noise;
X=[0 ; 0 ]; %初始状态
P=[1 0;0 1]; %状态协方差矩阵
F=[1 1;0 1]; %状态转移矩阵
Q=[0.0001,0;0 , 0.0001]; %状态转移协方差矩阵
H=[1,0]; %观测矩阵
R=1; %观测噪声协方差矩阵
figure;
hold on;
for i = 1:100
%基于上一状态预测当前状态  
X_ = F*X;
% 更新协方差  Q系统过程的协方差  这两个公式是对系统的预测
P_ = F*P*F'+Q;
% 计算卡尔曼增益
K = P_*H'/(H*P_*H'+R);
% 得到当前状态的最优化估算值  增益乘以残差
X = X_+K*(Z(i)-H*X_);
%更新K状态的协方差
P = (eye(2)-K*H)*P_;
scatter(X(1), X(2),4); %画点，横轴表示位置，纵轴表示速度
end

end

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
Matlab效果


 

 


其中 x 轴为位置，y轴为速度。 
在代码中，我们设定x的变化是 1:2:200，则速度就是2，可以由上图看到，值经过几次迭代，速度就基本上在 2 附近摆动，摆动的原因是我们加入了噪声。
接下来来看一个实际例子。
我们的数据为 data = [149.360 , 150.06, 151.44, 152.81,154.19 ,157.72]; 
这是运用光流法从视频中获取角点的实际x轴坐标，总共有6个数据，也就是代表了一个点的连续6帧的x轴坐标。接下来这个例子，我们将实现用5帧的数据进行训练，然后预测出第6帧的x轴坐标。
在上一个matlab例子中，我们的训练数据比较多，因此我们的初始状态设置为[0,0]，也就是位置为0，速度为0，在训练数据比较多的情况下，初始化数据为0并没有关系，因为我们在上面的效果图中可以看到，算法的经过短暂的迭代就能够发挥作用。
但在这里，我们的训练数据只有5帧，所以为了加快训练，我们将位置状态初始化为第一帧的位置，速度初始化为第二帧与第一帧之差。
测试代码：
KF.m
function [predData,dataX] = KF(dataZ)

%%
%
%   Description : kalmanFiltering
%   Author : Liulongpo
%   Time：2015-4-29 16:42:34
%
%%
Z = dataZ';
len = length(Z);
%Z=(1:2:200); %观测值  汽车的位置  也就是我们要修改的量
noise=randn(1,len); %方差为1的高斯噪声
dataX = zeros(2,len);
Z=Z+noise;
X=[Z(1) ; Z(2)-Z(1) ]; %初始状态  分别为 位置 和速度
P=[1 0;0 1]; %状态协方差矩阵
F=[1 1;0 1]; %状态转移矩阵
Q=[0.0001,0;0 , 0.0001]; %状态转移协方差矩阵
H=[1,0]; %观测矩阵
R=1; %观测噪声协方差矩阵
%figure;
%hold on;
for i = 1:len
%基于上一状态预测当前状态  
% 2x1  2x1
X_ = F*X;
% 更新协方差  Q系统过程的协方差  这两个公式是对系统的预测
%   2x1  2x1  1x2  2x2
P_ = F*P*F'+Q;
% 计算卡尔曼增益
K = P_*H'/(H*P_*H'+R);
% 得到当前状态的最优化估算值  增益乘以残差
X = X_+K*(Z(i)-H*X_);
%更新K状态的协方差
P = (eye(2)-K*H)*P_;
dataX(:,i) = [X(1);X(2)];
%scatter(X(1), X(2),4); %画点，横轴表示位置，纵轴表示速度
end
predData = F*X;
end
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
testKF.m
function testKF
%% 
clc
close all
%%
%data = load('D:\\a.txt');
%data = [149.360 , 150.06, 151.44, 152.81,154.19,157.72,157.47,159.33,153.66];
data = [149.360 , 150.06, 151.44, 152.81,154.19 ,157.72];
[predData , DataX] = KF(data');
error = DataX(1,:) - data;

i = 1:length(data);
figure
subplot 311
scatter(i,data,3),title('原始数据')
subplot 312
scatter(i,DataX(1,:),3),title('预测数据')
subplot 313
scatter(i,error,3),title('预测误差')
predData(1)
%{
scatter(i,error,3);
figure
scatter(i,data,3)
figure
scatter(i,predData(1,:),3)
%}
end
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
测试效果：


 

 


预测结果为： 155.7493 ，跟实际结果 157.72 仅有1.9 的误差，可以看到，卡尔曼滤波器算法对于少量数据的预测效果还是挺不错的。当然，预测位置的同时，我们也得到了预测速度。
参考文献： 视频: 卡尔曼滤波的原理以及在MATLAB中的实现

 
文章目录
 
Particle Filter
SIS
参考文献
 

 
Particle Filter
 
Kalman 滤波根据线性高斯模型可求得解析解，但是在非线性、非高斯的情况下无法得到解析解，对这类一般的情况，则要采用粒子滤波的方法
 
对于复杂的概率分布，可以通过一个简单的概率分布 $q(z)$ 作为桥梁进行采样（重要性采样），求得概率分布函数积分值（期望）:
 $$\begin{array}{rl}\mathbb{E}[f(z)]& ={\int }_{z}f(z)p(z)dz={\int }_{z}f(z)\frac{p(z)}{q(z)}q(z)dz\\ & =\frac{1}{N}\sum _{i=1}^{N}f(z_i)\boxed{\frac{p(z_i)}{q(z_i)}}\to weight\end{array}$$
 
通过对 $q(z)$ 采样，然后对每一个采样的样本应用权重就得到了期望的近似，为了概率分布的特性，还需要对权重进行归一化。
 
在滤波问题中，需要求解 $p(z_t|x_{1:t})$，其权重为：
 $$w_t^i=\frac{p(z_t^i|x_{1:t})}{q(z_t^i|x_{1:t})},i=1,2,\cdots ,N$$
 
若每一个时刻 $t$ 都需要采样 $N$ 个点，则每个时刻都要计算$N$个权重，使得计算较为繁杂，希望找到一个关于权重的递推公式，于是就引入了序列重要性采样（SIS）。
 
SIS
 
在 SIS (sequential importance sampling filter)中，解决的问题是 $p(z_{1:t}|x_{1:t})$。 $$w_t^i\propto \frac{p(z_{1:t}|x_{1:t})}{q(z_{1:t}|x_{1:t})}$$
 
根据 LDS 中的推导可以得到上式分子的递推式：
 $$\begin{array}{rl}p(z_{1:t}|x_{1:t})\propto p(x_{1:t},z_{1:t})& =p(x_t|z_{1:t},x_{1:t-1})p(z_{1:t},x_{1:t-1})\\ & =p(x_t|z_t)p(z_t|x_{1:t-1},z_{1:t-1})p(x_{1:t-1},z_{1:t-1})\\ & =p(x_t|z_t)p(z_t|z_{t-1})p(x_{1:t-1},z_{1:t-1})\\ & \propto p(x_t|z_t)p(z_t|z_{t-1})p(z_{1:t-1}|x_{1:t-1})\end{array}$$
 
对于分母(指定的提议分布)，直接分解为：
 $$q(z_{1:t}|x_{1:t})=q(z_t|z_{1:t-1},x_{1:t})q(z_{1:t-1}|x_{1:t-1})$$ 所以有： $$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)p(z_t|z_{t-1})p(z_{1:t-1}|x_{1:t-1})}{q(z_t|z_{1:t-1},x_{1:t})q(z_{1:t-1}|x_{1:t-1})}=\frac{p(x_t|z_t)p(z_t|z_{t-1})}{q(z_t|z_{1:t-1},x_{1:t})}{w}_{t-1}^i$$
 
得到计算权重的的步骤：
 
$t-1$ 时刻，采样完成并计算得到权重
$t$ 时刻，根据 $q(z_t|z_{1:t-1},x_{1:t})$ 进行采样得到 $z_t^i$，然后计算得到 $N$ 个权重
最后对权重归一化
 
SIS 算法会出现权值退化的情况，在一定时间后，可能会出现大部分权重很小，极少数较大，这是高维空间容易出现的情况，解决方法：
 
重采样，以权重作为概率分布，重新在已经采样的样本中采样，然所有样本赋予相同权重。这个方法的思路是将权重作为概率分布，然后得到累积密度函数，在累积密度上取点（阶梯函数）
选择合适的提议分布，$q(z_t|z_{1:t-1},x_{1:t})=p(z_t|z_{t-1})$，即：直接选择其状态转移概率分布，粒子$z_t^i\sim p(z_t|z_{t-1}^i)$，这也称为生成与测试方法：$$w_t^i\propto \frac{p(x_t|z_t)\cancel{p(z_t|z_{t-1})}}{\cancel{p(z_t|z_{t-1})}}{w}_{t-1}^i=p(x_t|z_t)w_{t-1}^i$$
 
采用重采样的 SIS 算法就是基本的粒子滤波算法。若选择上述提议分布，则称之为 SIR 算法。
 
参考文献
 
【1】粒子滤波
 【2】粒子滤波算法
 【3】Particle Filter Tutorial 粒子滤波：从推导到应用（四）
 【4】粒子滤波（Particle filter）算法简介及MATLAB实现
 【5】Particle Filter Tutorial 粒子滤波：从推导到应用（一）到（ 四 ）